/*HydroCreate create data for hydrogen and helium, 1 per coreload, called by CreatePoint 
 * in turn called after commands parsed */
#include <float.h>
#include <limits.h>
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "nhe1lvl.h"
#include "powi.h"
#include "he3tau.h"/* needed for NHE3TAU */
#include "he1dmp.h"
#include "hlmax.h"
#include "hwidth.h"
#include "taumin.h"
#include "ophf.h"
#include "heopfr.h"
#include "he1as.h"
#include "assav.h"
#include "hlife.h"
#include "hstat.h"
#include "he1stat.h"
#include "hphoto.h"
#include "he1tau.h"
#include "he1nxt.h"
#include "phe1lv.h"
#include "he1ind.h"
#include "hdeltat.h"
#include "he1deltat.h"
#include "he1nionryd.h"
#include "hydroeinsta.h"
#include "elmton.h"
#include "abscf.h"
#include "emlinejunk.h"
#include "refindex.h"
#include "getgf.h"
#include "ph1com.h"

void HydroCreate(void)
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  ipZ, 
	  j, 
	  lo, 
	  nHi, 
	  n2 ,
	  nLo;
	/* will be used to confirm that this is called once per coreload */
	static int nCalled=0;
	double HIonPoten, 
	  Potential, 
	  gHi, 
	  gLo, 
	  *hdamp/*[LMHLVL+1]*/, 
	  scrhe1, 
	  test, 
	  z4;
	double He1EnerRyd[NHE1LVL][NHE1LVL];

#	ifdef DEBUG_FUN
	fputs( "<+>HydroCreate()\n", debug_fp );
#	endif

	/* increment the nCalled counter, and confirm that it is exactly equal to 1 */
	++nCalled;
	assert( nCalled == 1 );

	/* hdamp is freed at end of sub */
	if( (hdamp=(double*)malloc(sizeof(double)*(unsigned)(nhlevel+1)) )==NULL )
	{
		fprintf( ioQQQ, " could not malloc hdamp\n" );
		puts( "[Stop in HydroCreate]" );
		exit(1);
	}

	/* now malloc the space for the hydrogenic lines, 
	 * but this must be done exactly one time per coreload */
	if( !lgHydroMalloc )
	{
		/*if( (HydroLines = (float ****)malloc(sizeof(float ***)*LIMELM )) ==NULL)*/
		if( (HydroLines = (EmLine ***)malloc(sizeof(EmLine **)*LIMELM )) ==NULL)
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc HydroLines 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}

		for( ipZ=0; ipZ < LIMELM; ++ipZ )
		{
			/* only grab core for elements that are turned on */
			if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
			{
				/*if( (HydroLines[ipZ]=(float***)malloc(sizeof(float **)*(unsigned)(nhlevel+1))) ==NULL)*/
				if( (HydroLines[ipZ]=(EmLine **)malloc(sizeof(EmLine *)*(unsigned)(nhlevel+1))) ==NULL)
				{
					fprintf( ioQQQ, 
						" HydroCreate could not malloc HydroLines 2\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}

				for( ipHi=1; ipHi <= nhlevel; ++ipHi )
				{
					/*if( (HydroLines[ipZ][ipHi]=(float**)malloc(sizeof(float *)*(unsigned)ipHi))==NULL)*/
					if( (HydroLines[ipZ][ipHi]=(EmLine*)malloc(sizeof(EmLine )*(unsigned)ipHi))==NULL)
					{
						fprintf( ioQQQ, 
							" HydroCreate could not malloc HydroLines 3\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
#if 0
					/* no longer needed since structure */
					/* arrays are dim'd nhlevel+1 */
					for( ipLo=IP1S; ipLo < ipHi; ++ipLo )
					{
						if( (HydroLines[ipZ][ipHi][ipLo] = (float*)malloc(sizeof(float)*NTA)) ==NULL)
						{
							fprintf( ioQQQ, 
								" HydroCreate could not malloc HydroLines 4\n" );
							puts( "[Stop in HydroCreate]" );
							exit(1);
						}
					}
#endif
				}
			}
		}

		/* This array is [LIMELM][nhlevel+1] */
		/* double **hcbolt == hcbolt[LIMELM][LMHLVL+1]; */
		if( (hydro.hcbolt = (double **)malloc(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hcbolt 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*double **hbn hbn[LIMELM][LMHLVL+1]*/;
		if( (hydro.hbn = (double **)malloc(sizeof(double *)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hbn 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}

		/*double ***hlbolt hlbolt[LIMELM][LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.hlbolt = (double ***)malloc(sizeof(double **)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hlbolt 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float ***esesc esesc[LIMELM][LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.esesc = (float ***)malloc(sizeof(float **)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.esesc 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float ***HRadNet HRadNet[LIMELM][LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.HRadNet = (float ***)malloc(sizeof(float **)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.HRadNet 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float ***hcol hcol[LIMELM][LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.hcol = (float ***)malloc(sizeof(float **)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hcol 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float ***hrec hrec[LIMELM][LMHLVL+1][3]*/;
		if( (hydro.hrec = (float ***)malloc(sizeof(float **)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hrec 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float **hcolnc hcolnc[LIMELM][LMHLVL+1]*/;
		if( (hydro.hcolnc = (float **)malloc(sizeof(float *)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hcolnc 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float **hn hn[LIMELM][LMHLVL+1]*/;
		if( (hydro.hn = (float **)malloc(sizeof(float *)*(unsigned)LIMELM ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.hn 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float **strkar strkar[LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.strkar = (float **)malloc(sizeof(float *)*(unsigned)(nhlevel+1) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.strkar 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		/*float **pestrk pestrk[LMHLVL+1][LMHLVL+1]*/;
		if( (hydro.pestrk = (float **)malloc(sizeof(float *)*(unsigned)(nhlevel+1) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" HydroCreate could not malloc hydro.pestrk 1\n" );
			puts( "[Stop in HydroCreate]" );
			exit(1);
		}
		for( ipZ=0; ipZ < LIMELM; ++ipZ )
		{
			if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
			{
				if( (hydro.hcbolt[ipZ] = 
					(double*)malloc(sizeof(double)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hcbolt 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hbn[ipZ] = 
					(double*)malloc(sizeof(double)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hbn 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hlbolt[ipZ] = 
					(double**)malloc(sizeof(double*)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hlbolt 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.esesc[ipZ] = 
					(float**)malloc(sizeof(float*)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.esesc 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.HRadNet[ipZ] = 
					(float**)malloc(sizeof(float*)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.HRadNet 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hcol[ipZ] = 
					(float**)malloc(sizeof(float*)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hcol 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hrec[ipZ] = 
					(float**)malloc(sizeof(float*)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hrec 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hcolnc[ipZ] = 
					(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hcolnc 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				if( (hydro.hn[ipZ] = 
					(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
				{
					fprintf( ioQQQ, 
						" HydroCreate could not second malloc hydro.hn 1\n" );
					puts( "[Stop in HydroCreate]" );
					exit(1);
				}
				for( ipLo=0; ipLo<= nhlevel ;++ipLo )
				{
					if( (hydro.strkar[ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not second malloc hydro.strkar 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.pestrk[ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not second malloc hydro.pestrk 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.hlbolt[ipZ][ipLo] = 
						(double*)malloc(sizeof(double)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not third malloc hydro.hlbolt 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.esesc[ipZ][ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not third malloc hydro.esesc 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.HRadNet[ipZ][ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not third malloc hydro.HRadNet 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.hcol[ipZ][ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)(nhlevel+1) ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not third malloc hydro.hcol 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
					if( (hydro.hrec[ipZ][ipLo] = 
						(float*)malloc(sizeof(float)*(unsigned)3 ))==NULL )
					{
						fprintf( ioQQQ, 
							" HydroCreate could not third malloc hydro.hrec 1\n" );
						puts( "[Stop in HydroCreate]" );
						exit(1);
					}
				}
			}
		}

		/* we will never do this again, in this coreload,
		 * following says never change number of levels in hydrogen atom again,
		 * future hydrogenic levels commands will be ignored*/
		lgHydroMalloc = TRUE;	

		/* main hydrogenic arrays, initialize so they will blow up if misused
		 * only do this when array created, since pointers to continuum energies
		 * will not need to be recreated */
		for( ipZ=0; ipZ < LIMELM; ipZ++ )
		{
			if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
			{
				/* arrays are dim'd nhlevel+1 */
				for( ipLo=IP1S; ipLo < nhlevel; ipLo++ )
				{
					/* set ENTIRE array to impossible values, in case of bad pointer */
					for( ipHi=ipLo+1; ipHi <= nhlevel; ipHi++ )
					{
						/*for( i=0; i<NTA; ++i )*/
						/*{*/
							/*HydroLines[ipZ][ipHi][ipLo].t[i] = -FLT_MAX;*/
						/*}*/
						EmLineJunk( &HydroLines[ipZ][ipHi][ipLo] );
					}
				}
			}
		}
	}


	/* upper level 4 and above, going to lower level 2 and above,
	 * come from the HyLife and branch product formed in routine pesc
	 * all other values come from here*/

	for( i=IP1S; i <= nhlevel; i++ )
	{
		hlmaxCom.hlmax[i] = 0.;
		hwidth.HLineWidth[i] = 0.;
	}


	/* main hydrogenic arrays, fill with sane values */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		/* must always do helium even if turned off */
		if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
		{
			/* charge to 4th power, needed for scaling laws */
			z4 = POW2(ipZ+1.);
			z4 *= z4;
			for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ )
			{
				if( ipLo == IP1S )
				{
					/* ground */
					nLo = 1;
					gLo = 2.;
				}
				else if( ipLo == IP2S )
				{
					/* 2s level */
					nLo = 2;
					gLo = 2.;
				}
				else if( ipLo == IP2P )
				{
					/* 2p level */
					nLo = 2;
					gLo = 6.;
				}
				else
				{
					nLo = ipLo;
					gLo = 2.*nLo*nLo;
				}
				for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
				{
					if( ipHi == IP2S )
					{
						nHi = 2;
						gHi = 2.;
					}
					else if( ipHi == IP2P )
					{
						nHi = 2;
						gHi = 6.;
					}
					else
					{
						nHi = ipHi;
						gHi = 2.*nHi*nHi;
					}
					/* pointers to transitions, ipLo, ipHi
					 * actual quantum numbers, nLo nHi
					 * atomic number or charge and stage: */
					HydroLines[ipZ][ipHi][ipLo].nelem = ipZ+1;
					HydroLines[ipZ][ipHi][ipLo].IonStg = ipZ+1;

					/* Dima's array has ionization potentials in eV, but not on same
					 * scale as cloudy itself*/
					/* extra factor accounts for this */
					HIonPoten = PH1COM.PH1[0][0][ipZ][0]/EVRYD* 0.9998787;
					assert(HIonPoten > 0.);

					/* transition energy in various units:
					 * 1e-10 is to stop 2s-2p from having zero energy */
					if( ipLo==IP2S && ipHi==IP2P )
					{
						HydroLines[ipZ][ipHi][ipLo].EnergyRyd = 1e-10f;
					}
					else
					{
						HydroLines[ipZ][ipHi][ipLo].EnergyRyd = 
							(float)(HIonPoten*
						  (  1./POW2((double)nLo) -1./POW2((double)nHi) ) );
					}
					assert(HydroLines[ipZ][ipHi][ipLo].EnergyRyd > 0.);

					HydroLines[ipZ][ipHi][ipLo].EnergyErg = 
						(float)(HydroLines[ipZ][ipHi][ipLo].EnergyRyd*
					  EN1RYD);
					assert(HydroLines[ipZ][ipHi][ipLo].EnergyErg > 0.);

					HydroLines[ipZ][ipHi][ipLo].EnergyK = 
						(float)(HydroLines[ipZ][ipHi][ipLo].EnergyRyd*
					  TE1RYD);
					assert(HydroLines[ipZ][ipHi][ipLo].EnergyK > 0.);

					/* the energy of the transition in wavenumbers */
					HydroLines[ipZ][ipHi][ipLo].EnergyWN = 
						(float)(HydroLines[ipZ][ipHi][ipLo].EnergyRyd/
					  WAVNRYD);
					assert(HydroLines[ipZ][ipHi][ipLo].EnergyWN > 0.);

					/* make following an air wavelength */
					HydroLines[ipZ][ipHi][ipLo].WLAng = 
						(float)(1.0e8/
					  HydroLines[ipZ][ipHi][ipLo].EnergyWN/
					  RefIndex(&HydroLines[ipZ][ipHi][ipLo]));
					assert(HydroLines[ipZ][ipHi][ipLo].WLAng > 0.);

					/* stat. weights and Einstein As & gfs: */
					HydroLines[ipZ][ipHi][ipLo].gLo = (float)gLo;

					HydroLines[ipZ][ipHi][ipLo].gHi = (float)gHi;

					/* transition prob, EinstA uses current H atom pointers */
					HydroLines[ipZ][ipHi][ipLo].Aul = 
						(float)(HydroEinstA(ipLo,ipHi)*z4);
					assert(HydroLines[ipZ][ipHi][ipLo].Aul > 0.);

					/* redistribution function = complete */
					HydroLines[ipZ][ipHi][ipLo].iRedisFun = -1;

					HydroLines[ipZ][ipHi][ipLo].gf = 
						(float)(GetGF(HydroLines[ipZ][ipHi][ipLo].Aul,
					  HydroLines[ipZ][ipHi][ipLo].EnergyWN,
					  HydroLines[ipZ][ipHi][ipLo].gHi));
					assert(HydroLines[ipZ][ipHi][ipLo].gf > 0.);

					/* derive the abs coef, call to function is gf, wl (A), g_low */
					HydroLines[ipZ][ipHi][ipLo].opacity = 
						(float)(abscf(HydroLines[ipZ][ipHi][ipLo].gf,
					  HydroLines[ipZ][ipHi][ipLo].EnergyWN,
					  HydroLines[ipZ][ipHi][ipLo].gLo));
					assert(HydroLines[ipZ][ipHi][ipLo].opacity > 0.);

					HydroLines[ipZ][ipHi][ipLo].Pesc = 1.0;

					/* destruction probability*/
					HydroLines[ipZ][ipHi][ipLo].Pdest = 0.0;

					/* upper and lower level populations */
					HydroLines[ipZ][ipHi][ipLo].PopHi = 0.0;
					HydroLines[ipZ][ipHi][ipLo].PopLo = 0.0;

					HydroLines[ipZ][ipHi][ipLo].xIntensity = 0.0;

					/* basic collision Info:
					 * these are not Mewe lines, they are high quality: */
					HydroLines[ipZ][ipHi][ipLo].cs1 = 0.0;

					/* following three will eventually be set by model atoms, but
					 * must have reasonable value at start */
					HydroLines[ipZ][ipHi][ipLo].pump = 0.;
					HydroLines[ipZ][ipHi][ipLo].AovTot = 0.;
					HydroLines[ipZ][ipHi][ipLo].cool = 0.;
					HydroLines[ipZ][ipHi][ipLo].heat = 0.;
					HydroLines[ipZ][ipHi][ipLo].ColOvTot = 0.;
					HydroLines[ipZ][ipHi][ipLo].PopOpc = 0.;
				}
			}

			/* total toptical depth in 1s - 2s */
			HydroLines[ipZ][IP2S][IP1S].TauTot = 0.;

			/* 2s - 2p cannot have zero opacity */
			HydroLines[ipZ][IP2P][IP2S].opacity = 1e-30f;

			/* Lya has special redistribution function */
			HydroLines[ipZ][IP2P][IP1S].iRedisFun = 1;
		}
	}

	/* following comes out very slightly off, correct here */
	HydroLines[ipHE][3][IP2S].WLAng = 1640.f;
	HydroLines[ipHE][3][IP2P].WLAng = 1640.f;

	/**************************************************************************
	 * now set HyLife for hydrogen itself, does not include Z^4 or Lyman lines 
	 **************************************************************************/
	hlife.HyLife[IP2P] = HydroLines[0][IP2P][IP1S].Aul;
	for( ipHi=3; ipHi <= nhlevel; ipHi++ )
	{
		hlife.HyLife[ipHi] = 0.;

		/* HyLife does not include trans to ground since case b branch used later */
		for( ipLo=IP2S; ipLo < ipHi; ipLo++ )
		{
			/* array of damping constants for hydrogen itself */
			hlife.HyLife[ipHi] += HydroLines[0][ipHi][ipLo].Aul;
		}
	}
	/* these are not defined and must never be used */
	hlife.HyLife[IP1S] = -FLT_MAX;
	hlife.HyLife[IP2S] = -FLT_MAX;

	/**************************************************************************
	 *  now fill in full array of damping constants 
	 **************************************************************************/
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
		{
			/* charge to 4th power, needed for scaling laws */
			z4 = POW2(ipZ+1.);
			z4 *= z4;
			/* put something in for 2s - 2p */
			HydroLines[ipZ][IP2P][IP2S].damprel = 1e-30f;

			for( ipHi=IP2P; ipHi <= nhlevel; ipHi++ )
			{
				for( ipLo=IP1S; ipLo < ipHi; ipLo++ )
				{
					/* get A's from H itself, scaling by charge */
					/* ipLnDampRel is related to the damping constant.  The number stored
					 * is the ratio Gamma lambda / 4 pi where lambda is the wavelength in
					 * cm (obtained by dividing by the energy in wavenumbers) and Gamma is
					 * the inverse lifetime, or sum of A's out of that level */
					HydroLines[ipZ][ipHi][ipLo].damprel = 
						(float)(
						/* this is sum of A's to n=2 for H by itself */
						(hlife.HyLife[ipHi]*z4 + 
						/* in following must add on the Lyman lines */
					  HydroLines[ipZ][ipHi][IP1S].Aul)/
					  PI4/HydroLines[ipZ][ipHi][ipLo].EnergyWN);

					assert(HydroLines[ipZ][ipHi][ipLo].damprel> 0.);
				}
			}
		}
	}

	/* 2p 1s is special since do not assume 2s 2p mixed */
	hdamp[2] = HydroLines[0][IP2P][IP1S].Aul*7.958e-6/0.75;

	/* NB hdamp is only dimensioned up through lmhlvl, so 2s-1s has no damp cnst */
	hlife.HyLife[IP2P] = HydroLines[0][IP2P][IP1S].Aul;

	for( ipHi=3; ipHi <= nhlevel; ipHi++ )
	{
		hlife.HyLife[ipHi] = 0.;
		for( ipLo=IP2S; ipLo < ipHi; ipLo++ )
		{
			/* array of lifetimes for hydrogen, scaled values used in HydroPesc to
			 * generate density dependent A's, Lyman lines are not included in this loop */
			hlife.HyLife[ipHi] += HydroLines[0][ipHi][ipLo].Aul;
		}
		/* lyman lines not included in HyLife, since branching ratios only defined
		 * to get case B right.  for damp now add on lyman lines */
		hdamp[ipHi] = (hlife.HyLife[ipHi] + HydroLines[0][ipHi][IP1S].Aul)*7.958e-6/
		  (1. - 1./POW2( (float)ipHi));
	}

	/* levels not possible */
	hdamp[0] = 0.;
	hdamp[1] = 0.;

	/* define excitation and ionization temperatures for hydrogen */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		/* define these for H and He always */
		if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
		{
			/* geneate pointers that will blow if used */
			for( ipLo=IP1S; ipLo<nhlevel; ++ipLo )
			{
				for( ipHi=ipLo+1; ipHi<=nhlevel; ++ipHi )
				{
					HydroLines[ipZ][IP2P][IP1S].ipCont = INT_MIN;
				}
			}

			/* get ground ionization threshold from Dima's fits */
			Potential = PH1COM.PH1[0][0][ipZ][0]/EVRYD* 0.9998787;

			for( ipHi=IP2P; ipHi <= nhlevel; ipHi++ )
			{
				hydro.HLevNIonRyd[ipZ][ipHi] = 
					(float)(Potential/POW2((float)ipHi));
				HdeltaT.HCionTe[ipZ][ipHi] = 
					(float)(TE1RYD*hydro.HLevNIonRyd[ipZ][ipHi]);
			}

			/* 2s done here */
			hydro.HLevNIonRyd[ipZ][IP2S] = hydro.HLevNIonRyd[ipZ][IP2P];
			HdeltaT.HCionTe[ipZ][IP2S] = 
				(float)(TE1RYD*hydro.HLevNIonRyd[ipZ][IP2P]);

			/* 1s done here */
			hydro.HLevNIonRyd[ipZ][IP1S] = (float)Potential;
			HdeltaT.HCionTe[ipZ][IP1S] = 
				(float)(TE1RYD*hydro.HLevNIonRyd[ipZ][IP1S]);
		}
	}

	/* set statistical weights for hydrogen atom */
	hstat.HStatWght[IP1S] = 2.;
	hstat.HStatWght[IP2S] = 2.;
	hstat.HStatWght[IP2P] = 6.;
	for( ipHi=3; ipHi <= nhlevel; ipHi++ )
	{
		/* NB - this is used for He1 by scaling later on -
		 * if this is killed then move code to HeI
		 * fill in stat weights for hydrogen */
		hstat.HStatWght[ipHi] = (float)(2.*ipHi*ipHi);
	}

	/* statistical weights for helium singlet atom and ion */
	for( ipHi=0; ipHi < NHE1LVL; ipHi++ )
	{
		n2 = (ipHi+1)*(ipHi+1);
		He1NIonRyd.He1IonRyd[ipHi] = 1.f/n2;
		He1deltaT.He1CionTe[ipHi] = (float)(TE1RYD*He1NIonRyd.He1IonRyd[ipHi]);
		he1statCOM.he1stat[ipHi] = (float)n2;
	}

	/* ground and n=2 not hydrogenic */
	He1NIonRyd.He1IonRyd[0] = 1.80802f;
	He1NIonRyd.He1IonRyd[1] = 0.2478f;
	He1NIonRyd.He1IonRyd[NHE1LVL] = He1NIonRyd.He1IonRyd[1];
	he1statCOM.he1stat[NHE1LVL] = 1.;
	he1statCOM.he1stat[1] = 3.;

	for( ipHi=1; ipHi < NHE1LVL; ipHi++ )
	{
		for( lo=0; lo <= (ipHi - 1); lo++ )
		{
			He1EnerRyd[lo][ipHi] = He1NIonRyd.He1IonRyd[lo] - 
			  He1NIonRyd.He1IonRyd[ipHi];
			He1deltaT.He1deltaTe[lo][ipHi] = 
				(float)(TE1RYD*He1EnerRyd[lo][ipHi]);
		}
	}

	for( i=0; i < (NHE1LVL + 1); i++ )
	{
		he1ind.he1qin[i] = 0.;
		he1ind.he1rin[i] = 0.;
	}

	/* these are all the helium triplet lines */
	for( i=0; i < NHE3TAU; i++ )
	{
		he3tau[i].TauIn = tauminCom.taumin;
		he3tau[i].TauTot = 1e20f;
		he3tau[i].Pesc = 1.;
		he3tau[i].FracInwd = 1.;
		/* 2 is atomic number of helium */
		he3tau[i].nelem = 2;
		he3tau[i].IonStg = 1;
	}
	/* now set wavelengths */
	he3tau[IPT10830-1].WLAng = 10830.;
	he3tau[IPT3889-1].WLAng = 3889.;
	he3tau[IPT5876-1].WLAng = 5876.;
	he3tau[IPT7065-1].WLAng = 7065.;
	/* now set A's */
	he3tau[IPT10830-1].Aul = 1.02e7f;
	he3tau[IPT3889-1].Aul = 9.48e6f;
	he3tau[IPT5876-1].Aul = 7.06e7f;
	he3tau[IPT7065-1].Aul = 2.78e7f;

	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		if( ipZ < 2 || elmton.lgElmtOn[ipZ] )
		{
			for( i=IP1S; i <= nhlevel; i++ )
			{
				ophf.HOpacRatio[ipZ][i] = 1.;
				hydro.hbn[ipZ][i] = -FLT_MAX;
				/* this is used for pressure before hatom defined*/
				hydro.hn[ipZ][i] = 0.;
				HPhoto.hgamnc[ipZ][i] = -DBL_MAX;
				HPhoto.cinduc[ipZ][i] = -DBL_MAX;
				HPhoto.finduc[ipZ][i] = -DBL_MAX;
				HPhoto.rinduc[ipZ][i] = -DBL_MAX;
				hydro.hcolnc[ipZ][i] = -FLT_MAX;
				hydro.hrec[ipZ][i][ipRecRad] = -FLT_MAX;
				hydro.hrec[ipZ][i][ipRecNetEsc] = -FLT_MAX;
				hydro.hrec[ipZ][i][ipRecEsc] = -FLT_MAX;

				for( j=i + 1; j <= nhlevel; j++ )
				{
					hydro.hcol[ipZ][i][j] = -FLT_MAX;
				}
			}
		}
	}
	/* ground state of H and He is different since totally determine
	 * their own opacities */
	ophf.HOpacRatio[ipH][IP1S] = 1e-5f;
	ophf.HOpacRatio[ipHE][IP1S] = 1e-5f;

	/* this is hydrogenic screening factor for HeI singlet ground state */
	scrhe1 = 1.808;
	for( i=0; i < (NHE1LVL - 1); i++ )
	{
		/* damping constant scales as A*lambda, or z */
		/* for hydrogen, hdamp[1] is 1s, [2] is 2,
		 * so in following changed i_ to i in hdamp */
		he1dmpCOM.he1dmp[i] = (float)hdamp[i+1];
		phe1lv.he1gam[i] = 0.;
		phe1lv.he1n[i] = 1e-35f;
		phe1lv.he1clc[i] = 0.;
		phe1lv.he1rec[0][i] = 0.;
		phe1lv.he1rec[1][i] = 1.;
		phe1lv.he1rec[2][i] = 1.;
		/* NHE1LVL is 9 in nhe1lvl.h */
		for( j=i + 1; j < NHE1LVL; j++ )
		{
			phe1lv.he1cll[i][j] = 0.;

			/* scale HeII A's from hydrogen */
			he1as.He1EinA[i][j] = (float)(scrhe1*AsSav.HeAsSav[i][j]);

			he1tauCOM.he1opc[i][j] = 
				(float)(he1as.He1EinA[i][j]*2.249e-26*
				/* HStatWght[0] is 1s, 1 is 2s, etc,
				 * it was from 0 in the original fortran so there
				 * was no ofset in the c version */
			  hstat.HStatWght[j+1]/hstat.HStatWght[i+1]*
			  /* HEnerRyd has 1 for 1s, so is offset by one in the c version,
			   * I changed the i_ and j_ to i,j */
				/*>>>>chng 99 jan 26, crashed when H atom smaller than He atom*/
			  /*powi(RYDLAM/HEnrRyd.HEnerRyd[i+1][j+1],3));*/
			  powi(RYDLAM/(HIONPOT*(1./((i+1.)*(i+1.)) - 1./((j+1.)*(j+1.)) ) ),3));
			/* destruction probabilities */
			phe1lv.he1dst[i][j] = 0.;
			phe1lv.he1mis[i][j] = 0.;
		}
		scrhe1 = 1.;
	}

	for( i=0; i < NHE1LVL; i++ )
	{
		heopfr.ophe1f[i] = 0.;
	}

	for( lo=0; lo < (NHE1LVL - 1); lo++ )
	{
		if( lo == 1 )
		{
			test = 8.;
		}
		else
		{
			test = he1statCOM.he1stat[lo];
		}
		for( ipHi=lo + 1; ipHi < NHE1LVL; ipHi++ )
		{
			/*lint -e771 He1EnerRyd possibley not init */
			he1tauCOM.he1opc[lo][ipHi] = (float)(he1as.He1EinA[lo][ipHi]*
			  2.24484e-26*he1statCOM.he1stat[ipHi]/test*powi(RYDLAM/
			  He1EnerRyd[lo][ipHi],3));
			/*lint +e771 He1EnerRyd possibley not init */
		}
	}

	/* HeI ground state is special */
	he1tauCOM.he1opc[0][1] = 2.42e-8f;
	he1as.He1EinA[0][1] = 17.99e8f;
	he1dmpCOM.he1dmp[1] *= 1.37f;
	phe1lv.he1n[NHE1LVL-1] = 1e-30f;
	phe1lv.he1gam[NHE1LVL-1] = 0.;
	phe1lv.he1clc[NHE1LVL-1] = 0.;
	phe1lv.he1rec[0][NHE1LVL-1] = 0.;
	phe1lv.he1rec[1][NHE1LVL-1] = 1.;
	phe1lv.he1rec[2][NHE1LVL-1] = 1.;
	he1dmpCOM.he1dmp[NHE1LVL-1] = (float)hdamp[8];
	phe1lv.he12s = 1e-30f;
	phe1lv.he12p = 1e-30f;

	/* zero out inner optical depths */
	for( i=0; i < (NHE1LVL - 1); i++ )
	{
		for( j=0; j < NHE1LVL; j++ )
		{
			he1nxtCOM.he1nxt[i][j] = tauminCom.taumin;
			he1tauCOM.he1tau[i][j] = tauminCom.taumin;
		}
	}
	t206.TauIn = tauminCom.taumin;
	t206.TauTot = tauminCom.taumin;
	t206.Pesc = 1.;
	t206.FracInwd = 0.5;
	t206.Aul = 1.976e6;


	/* now free the local space malloced above */
	free(hdamp);

#	ifdef DEBUG_FUN
	fputs( " <->HydroCreate()\n", debug_fp );
#	endif
	return;
}