/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*iso_create create data for hydrogen and helium, 1 per coreload, called by ContCreatePointers 
 * in turn called after commands parsed */
/*iso_zero zero data for hydrogen and helium */
#include "cddefines.h"
#include "physconst.h"
#include "lines_service.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "helike.h"
#include "opacity.h"
#include "hydroeinsta.h"
#include "dense.h"
#include "atmdat.h"
#include "iso.h"

/*iso_zero zero data for hydrogen and helium */
static void iso_zero(void);

void iso_create(void)
{
	long int
	  ipHi, 
	  ipISO ,
	  ipLo, 
	  nelem, 
	  n ,
	  need,
	  nHi, 
	  nLo,
	  NumHeLinesSav = 0;
	long jj, ii, j;

	static int nCalled = 0;

	double HIonPoten, 
	  Potential, 
	  gHi, 
	  gLo, 
	  *hdamp/*[LMHLVL+1]*/, 
	  z4;

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

	/* > 1 if not first call, then just zero arrays out */
	if( nCalled > 0 )
	{
		iso_zero();
#		ifdef DEBUG_FUN
		fputs( " <->iso_create()\n", debug_fp );
#		endif
		return;
	}

	/* this is first call, increment the nCalled counterso never do this again */
	++nCalled;

	/* >>chng 02 feb 21, the entire He-like iso-sequence will have nCollapsed n-grouped
	 * levels in addition to the number of l-resolved levels 
	 * atom he-like levels adjusts number of l-resolved levels 
	 * atom he-like levels collapsed adjusts nCollapsed */
	if( helike.lgSetBenjamin )
	{
		/* This is an atom with maximum principal quantum number = 5	*/
		iso.numLevels_max[ipHE_LIKE][ipHELIUM] = 31;
		iso.numLevels_local[ipHE_LIKE][ipHELIUM] = iso.numLevels_max[ipHE_LIKE][ipHELIUM];
		iso.nCollapsed_max[ipHE_LIKE][ipHELIUM] = 0;
		iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] = 5;
		/* >>chng 06 jan 2, exclude highest level from print out. */
		iso.numPrintLevels[ipHE_LIKE][ipHELIUM] = iso.numLevels_max[ipHE_LIKE][ipHELIUM];

		for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
		{
			/* turn everything else off.	*/
			dense.lgElmtOn[nelem] = FALSE;
		}
	}
	else if( !helike.lgCompileRecomb )
	{
		for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
		{
			/* only grab core for elements that are turned on */
			if( nelem==ipHELIUM || dense.lgElmtOn[nelem] )
			{
				/* increase numPrintLevels while at it.	*/
				ASSERT( iso.nCollapsed_max[ipHE_LIKE][nelem] > 0);

				/* >>chng 06 jan 2, exclude collapsed levels from print out. */
				iso.numPrintLevels[ipHE_LIKE][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];

				/* now add on collapsed part...has to be done here because we don't know
				 * what order the commands for resolved and collapsed will enter parser. */
				iso.numLevels_max[ipHE_LIKE][nelem] += iso.nCollapsed_max[ipHE_LIKE][nelem];
				iso.numLevels_local[ipHE_LIKE][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
			}
		}
	}

	
	/* hdamp is freed at end of sub */
	if( (hdamp=(double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][ipHYDROGEN])) )==NULL )
		BadMalloc();

	/* these are the statistical weights of the ions */
	iso.stat_ion[ipH_LIKE] = 1.f;
	iso.stat_ion[ipHE_LIKE] = 2.f;

	/* create space for arrays of indices for ionization edges within continuum array */
	/* iso.ipOpac[NISO][LIMELM][NHPLPHOT] */
	if( (iso.ipOpac = (long int ***)MALLOC(sizeof(long **)*NISO ))==NULL )
		BadMalloc();

	/* Master array of levels...iso.quant_desig[ipHE_LIKE][nelem][ipLevel].n, .l, and .s	*/
	if( (iso.quant_desig = (Elevels***)MALLOC(sizeof(Elevels**)*(unsigned)NISO ) )==NULL )
		BadMalloc();

	if( (iso.ipIsoLevNIonCon = (long int ***)MALLOC(sizeof(long **)*NISO ))==NULL )
		BadMalloc();

	if( (iso.xIsoLevNIonRyd = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.stat = (float ***)MALLOC(sizeof(float **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.RadRecomb = (double ****)MALLOC(sizeof(double ***)*NISO )) == NULL )
		BadMalloc();

	if( (iso.ExtraLymanLines = (EmLine ***)MALLOC(sizeof(EmLine **)*NISO )) ==NULL)
		BadMalloc();

	if( (iso.Boltzmann = (double ****)MALLOC(sizeof(double ***)*NISO )) == NULL )
		BadMalloc();

	if( (iso.ConBoltz = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.Pop2Ion = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.RateLevel2Cont = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.RateCont2Level = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.ConOpacRatio = (float ***)MALLOC(sizeof(float **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.DepartCoef = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.gamnc = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.RecomInducRate = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.RecomInducCool_Coef = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.PhotoHeat = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.PopLTE = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	if( (iso.ColIoniz = (double ***)MALLOC(sizeof(double **)*NISO )) == NULL )
		BadMalloc();

	/* the hydrogen and helium like iso-sequences */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		long int nmax;

		/* Master array of levels...iso.quant_desig[ipHE_LIKE][nelem][ipLevel].n, .l, and .s	*/
		if( (iso.quant_desig[ipISO] = (Elevels**)MALLOC(sizeof(Elevels*)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		if( (iso.ExtraLymanLines[ipISO] = (EmLine **)MALLOC(sizeof(EmLine *)*LIMELM )) ==NULL)
			BadMalloc();

		if( (iso.ipOpac[ipISO] = (long int **)MALLOC(sizeof(long *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.ipIsoLevNIonCon[ipISO] = (long int **)MALLOC(sizeof(long *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.xIsoLevNIonRyd[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.stat[ipISO] = (float **)MALLOC(sizeof(float*)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.RadRecomb[ipISO] = (double ***)MALLOC(sizeof(double **)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.Boltzmann[ipISO] = (double ***)MALLOC(sizeof(double **)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.ConBoltz[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.Pop2Ion[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.RateLevel2Cont[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.RateCont2Level[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.ConOpacRatio[ipISO] = (float **)MALLOC(sizeof(float *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.DepartCoef[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM ))==NULL )
			BadMalloc();

		if( (iso.gamnc[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		if( (iso.RecomInducRate[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		if( (iso.RecomInducCool_Coef[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		if( (iso.PhotoHeat[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		if( (iso.PopLTE[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		if( (iso.ColIoniz[ipISO] = (double **)MALLOC(sizeof(double *)*LIMELM )) == NULL )
			BadMalloc();

		/* this will remember the greatest number of levels for each iso sequence,
		 * and is used to allocate the stat weight array after this loop */
		nmax = 0;
		for( nelem=ipISO; nelem < LIMELM; ++nelem )
		{
			/* only grab core for elements that are turned on */
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				nmax = MAX2( nmax , iso.numLevels_max[ipISO][nelem] );

				/*fprintf(ioQQQ,"assert failll %li\t%li\t%li\n", ipISO, nelem , iso.numLevels_max[ipISO][nelem] );*/
				ASSERT( iso.numLevels_max[ipISO][nelem] > 0 );
			
				/* Allocate space so that array iso.quant_desig[ipHE_LIKE] contains proper number of elements for each element. */
				if( (iso.quant_desig[ipISO][nelem] = 
					(Elevels*)MALLOC((unsigned)iso.numLevels_max[ipISO][nelem]*sizeof(Elevels)) ) == NULL )
					BadMalloc();

				if( (iso.ExtraLymanLines[ipISO][nelem] = 
					(EmLine *)MALLOC(sizeof(EmLine )*(unsigned)iso.nLyman[ipISO] )) ==NULL)
					BadMalloc();

				if( (iso.ipOpac[ipISO][nelem] = 
					(long int *)MALLOC(sizeof(long)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.ipIsoLevNIonCon[ipISO][nelem] = 
					(long int *)MALLOC(sizeof(long)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.xIsoLevNIonRyd[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.stat[ipISO][nelem] = (float *)MALLOC(sizeof(float)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
					BadMalloc();

				if( (iso.RadRecomb[ipISO][nelem] = 
					(double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.Boltzmann[ipISO][nelem] = 
					(double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.ConBoltz[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.DepartCoef[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.Pop2Ion[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.RateLevel2Cont[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.RateCont2Level[ipISO][nelem] = 
					(double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.ConOpacRatio[ipISO][nelem] = 
					(float *)MALLOC(sizeof(float)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				if( (iso.gamnc[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				if( (iso.RecomInducRate[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				if( (iso.RecomInducCool_Coef[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				if( (iso.PhotoHeat[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				if( (iso.PopLTE[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				if( (iso.ColIoniz[ipISO][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(iso.numLevels_max[ipISO][nelem]) )) == NULL )
					BadMalloc();

				n = 0;
				/* this needs n = 0, but RadRecomb will crash since malloc(0), not used */
				if( (iso.RadRecomb[ipISO][nelem][n] = (double *)MALLOC(sizeof(double)*(unsigned)3 ))==NULL )
					BadMalloc();
				for( n=1; n<iso.numLevels_max[ipISO][nelem]; ++n )
				{
					if( (iso.RadRecomb[ipISO][nelem][n] = (double *)MALLOC(sizeof(double)*(unsigned)3 ))==NULL )
						BadMalloc();

					/* sec to last dim is upper level n,
					 * last dim of this array is lower level, will go from 0 to n-1 */
					if( (iso.Boltzmann[ipISO][nelem][n] = (double *)MALLOC(sizeof(double)*(unsigned)n ))==NULL )
						BadMalloc();
				}
			}
		}
	}


	/* >>chng 06 mar 03, sizeof(EmLine **) to *** */
	if( (EmisLines = (EmLine ****)MALLOC(sizeof(EmLine ***)*NISO )) ==NULL)
		BadMalloc();

	/* For large helium model calculations, the EmisLines structure gets very large.  
	 * Only allocate a small number of lines if doing a huge case B calculation, 
	 * because most of the details of the structure will not be needed, and the 
	 * calculation will be far less memory requirements. */
	if( helike.lgHugeCaseB == TRUE )
	{
		NumHeLinesSav = iso.numLevels_max[ipHE_LIKE][ipHELIUM];
		iso.numLevels_max[ipHE_LIKE][ipHELIUM] = 32;
	}

	/* the hydrogen and helium like iso-sequences */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		long int nmax;

		if( (EmisLines[ipISO] = (EmLine ***)MALLOC(sizeof(EmLine **)*LIMELM ))==NULL )
			BadMalloc();

		/* this will remember the greatest number of levels for each iso sequence,
		 * and is used to allocate the stat weight array after this loop */
		nmax = 0;
		for( nelem=ipISO; nelem < LIMELM; ++nelem )
		{
			/* only grab core for elements that are turned on */
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				nmax = MAX2( nmax , iso.numLevels_max[ipISO][nelem] );

				/*fprintf(ioQQQ,"assert failll %li\t%li\t%li\n", ipISO, nelem , iso.numLevels_max[ipISO][nelem] );*/
				ASSERT( iso.numLevels_max[ipISO][nelem] > 0 );
			
				if( (EmisLines[ipISO][nelem] = 
					(EmLine **)MALLOC(sizeof(EmLine *)*(unsigned)(iso.numLevels_max[ipISO][nelem]) ))==NULL )
						BadMalloc();

				n = 0;
				for( n=1; n<iso.numLevels_max[ipISO][nelem]; ++n )
				{
					if( (EmisLines[ipISO][nelem][n] = (EmLine *)MALLOC(sizeof(EmLine)*(unsigned)n ))==NULL )
						BadMalloc();
				}
			}
		}
	}

	/* now MALLOC the space for the hydrogenic lines, 
	 * but this must be done exactly one time per coreload */
	if( !lgHydroMalloc )
	{

		if( (He1AutoLines = (EmLine **)MALLOC(sizeof(EmLine *)*LIMELM )) ==NULL)
			BadMalloc();

		/* the helium like iso-sequence */
		for( nelem=1; nelem < LIMELM; ++nelem )
		{
			/* only allocate space for elements that are turned on */
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				if( (He1AutoLines[nelem]=(EmLine *)MALLOC(sizeof(EmLine )*(unsigned)(4))) ==NULL)
					BadMalloc();
				
			}
		}

		/*float **strkar strkar[LMHLVL+1][LMHLVL+1]*/
		if( (hydro.strkar = (float **)MALLOC(sizeof(float *)*(unsigned)(iso.numLevels_max[ipH_LIKE][ipHYDROGEN]) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" iso_create could not MALLOC hydro.strkar 1\n" );
			puts( "[Stop in iso_create]" );
			cdEXIT(EXIT_FAILURE);
		}
		/*float **pestrk pestrk[LMHLVL+1][LMHLVL+1]*/
		if( (hydro.pestrk = (double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels_max[ipH_LIKE][ipHYDROGEN]) ) )==NULL )
		{
			fprintf( ioQQQ, 
				" iso_create could not MALLOC hydro.pestrk 1\n" );
			puts( "[Stop in iso_create]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now do second dimension of previous two */
		for( ipLo=0; ipLo< iso.numLevels_max[ipH_LIKE][ipHYDROGEN] ;++ipLo )
		{
			if( (hydro.strkar[ipLo] = 
				(float*)MALLOC(sizeof(float)*(unsigned)(iso.numLevels_max[ipH_LIKE][ipHYDROGEN]) ))==NULL )
			{
				fprintf( ioQQQ, 
					" iso_create could not second MALLOC hydro.strkar 1\n" );
				puts( "[Stop in iso_create]" );
				cdEXIT(EXIT_FAILURE);
			}
			if( (hydro.pestrk[ipLo] = 
				(double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][ipHYDROGEN]) ))==NULL )
			{
				fprintf( ioQQQ, 
					" iso_create could not second MALLOC hydro.pestrk 1\n" );
				puts( "[Stop in iso_create]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
		/*>>chng 03 oct 08, moved from nsset */
		jj = 0;
		/* define array of quantum numbers needed for Puetter Stark broadening */
		/* what should ipHYDROGEN be */
		for( n=ipH1s; n < (iso.numLevels_max[ipH_LIKE][ipHYDROGEN] - 1); n++ )
		{
			if( n == 0 )
			{
				ii = 1;
			}
			else if( n == 1 )
			{
				ii = 2;
			}
			else
			{
				ii = n;
			}
			for( j=n + 1; j < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; j++ )
			{
				if( j == 1 )
				{
					jj = 2;
				}
				else
				{
					jj = jj;
				}
				hydro.strkar[n][j] = (float)pow(((float)ii*jj),1.2f);
				hydro.pestrk[n][j] = 0.;
			}
		}

		/* 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( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				/* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */
				for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipH_LIKE][nelem]-1; ipLo++ )
				{
					/* set ENTIRE array to impossible values, in case of bad pointer */
					for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
					{
						EmLineJunk( &EmisLines[ipH_LIKE][nelem][ipHi][ipLo] );
					}
				}
				/* these are the extra Lyman lines */
				for( ipHi=2; ipHi < iso.nLyman[ipH_LIKE]; ipHi++ )
				{
					EmLineJunk( &iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi] );
				}
			}
		}
		/* the main he-like arrays */
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				/* set ENTIRE array to impossible values, in case of bad pointer */
				for( ipHi=1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
				{
					/* arrays are dim'd iso.numLevels_max[ipHE_LIKE][nelem] */
					for( ipLo=0; ipLo < ipHi; ipLo++ )
					{
						EmLineJunk( &EmisLines[ipHE_LIKE][nelem][ipHi][ipLo] );
					}
				}
				/* these are the extra Lyman lines */
				for( ipHi=2; ipHi < iso.nLyman[ipHE_LIKE]; ipHi++ )
				{
					EmLineJunk( &iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi] );
				}
			}
		}
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				/* There are four autoionization levels to be considered. */
				for( ipHi = 0; ipHi < 4 ; ipHi++ )
				{
					/* set ENTIRE array to impossible values, in case of bad pointer */
					EmLineJunk( &He1AutoLines[nelem][ipHi] );
				}
			}
		}
	}

	/* 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*/

	/* main hydrogenic arrays, fill with sane values */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* must always do helium even if turned off */
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			/* charge to 4th power, needed for scaling laws */
			z4 = POW2(nelem+1.);
			z4 *= z4;

			/* 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][nelem][0]/EVRYD* 0.9998787;
			ASSERT(HIonPoten > 0.);

			/*TODO	2	this will not work if highest level is resolved */
			iso.quant_desig[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem] - 1].n = iso.numLevels_max[ipH_LIKE][nelem] - 1;
			iso.quant_desig[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem] - 1].l = -1;
			iso.quant_desig[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem] - 1].s = -1;

			for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
			{
				/* test code to see what happens with lyman lines in case a */
				if( ipLo == ipH1s )
				{
					/* ground */
					nLo = 1;
					gLo = 2.;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].n = 1;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].l = 0;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].s = 0;
				}
				else if( ipLo == ipH2s )
				{
					/* 2s level */
					nLo = 2;
					gLo = 2.;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].n = 2;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].l = 0;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].s = 0;
				}
				else if( ipLo == ipH2p )
				{
					/* 2p level */
					nLo = 2;
					gLo = 6.;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].n = 2;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].l = 1;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].s = 0;
				}
				else
				{
					nLo = ipLo;
					gLo = 2.*nLo*nLo;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].n = ipLo;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].l = -1;
					iso.quant_desig[ipH_LIKE][nelem][ipLo].s = -1;
				}
				iso.stat[ipH_LIKE][nelem][ipLo] = (float)gLo;
				for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
				{
					double EnergyRyd;

					if( ipHi == ipH2s )
					{
						nHi = 2;
						gHi = 2.;
					}
					else if( ipHi == ipH2p )
					{
						nHi = 2;
						gHi = 6.;
					}
					else
					{
						nHi = ipHi;
						gHi = 2.*nHi*nHi;
					}
					iso.stat[ipH_LIKE][nelem][ipHi] = (float)gHi;
					/* pointers to transitions, ipLo, ipHi
					 * actual quantum numbers, nLo nHi
					 * atomic number or charge and stage: */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].nelem = (int)(nelem+1);
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].IonStg = (int)(nelem+1);

					/* transition energy in various units:
					 * 1e-10 is to stop 2s-2p from having zero energy */
					/* >>chng02 feb 09, use zero energy as sentinal of non-existant transition */
					if( ipLo==ipH2s && ipHi==ipH2p )
					{
						/*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyRyd = 1e-10f;*/
						EnergyRyd = 0.;
					}
					else
					{
						EnergyRyd = (HIonPoten*
						  (  1./POW2((double)nLo) -1./POW2((double)nHi) ) );
					}
					/* >>chng 02 feb 09, change test to >= 0 since we now use 0 for 2s-2p */
					ASSERT(EnergyRyd >= 0.);

					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg = 
						(float)(EnergyRyd*
					  EN1RYD);
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg >= 0.);

					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyK = 
						(float)(EnergyRyd*
					  TE1RYD);
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyK >= 0.);

					/* the energy of the transition in wavenumbers */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN = 
						(float)(EnergyRyd/ WAVNRYD);
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN >= 0.);

					if( ipLo==ipH2s && ipHi==ipH2p )
					{
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].WLAng = 0.;
					}
					else
					{
						/* make following an air wavelength */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].WLAng = 
							(float)(1.0e8/
						  EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN/
						  RefIndex( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN));
						ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].WLAng > 0.);
					}

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

					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gHi = (float)gHi;

					/* transition prob, EinstA uses current H atom indices */
					if( ipHi==ipH2s && ipLo==ipH1s )
					{
						/* two photon transition scales by Z**6	*/
						/* >>refer	H-like	A	Lipesles, M., Novick, R., & Told, N., 1966,
						 * >>refercon	Phys Rev Let, 15, 690 */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul = 
							(float)(HydroEinstA(ipLo,ipHi)*pow(nelem+1.,6.));
					}
					else
					{
						/* all others by Z**4	*/
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul = 
							(float)(HydroEinstA(ipLo,ipHi)*z4);
					}
					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul > 0.);

					/* >>chng 01 feb 27, had been following, for all lines except Lya */
					/* redistribution function = complete 
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].iRedisFun = -1; - this is same as ipCRD */
					if( ipLo == ipH1s && ipHi == ipH2p )
					{
						/* this is La, special redistribution, default is ipLY_A */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].iRedisFun = iso.ipLyaRedist[ipH_LIKE];
					}
					else if( ipLo == ipH1s && ipHi != ipH2p )
					{
						/* these are rest of Lyman lines, 
						 * complete redistribution, doppler core only, K2 core, default ipCRD */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].iRedisFun = iso.ipResoRedist[ipH_LIKE];
					}
					else
					{
						/* all lines coming from excited states, default is complete redis with wings, ipCRDW*/
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].iRedisFun = iso.ipSubRedist[ipH_LIKE];
					}
					if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN <= 0. )
					{
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gf = 0.;
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity = 0.;
					}
					else
					{
						/* all subordinate lines, CRD with wings */
						/* >>chng 01 feb 27, had been -1, crd with core only,
						* change to crd with wings as per discussion with Ivan Hubeny */
						/*EmisLines[ipH_LIKE][nelem][ipHi][ipLo].iRedisFun = -1; = ipCRD*/

						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gf = 
							(float)(GetGF(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul,
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN,
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gHi));
						ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gf > 0.);

						/* derive the abs coef, call to function is gf, wl (A), g_low */
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity = 
							(float)(abscf(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gf,
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN,
						EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gLo));
						ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity > 0.);
					}

					/* create array index that will blow up */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ipCont = INT_MIN;
				}
			}
			/* these are the extra Lyman lines */
			nLo = 1;
			for( ipHi=2; ipHi < iso.nLyman[ipH_LIKE]; ipHi++ )
			{
				double EnergyRyd;
				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].nelem = (int)(nelem+1);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].IonStg = (int)(nelem+1);

				EnergyRyd = 
							(float)(HIonPoten*
						  (  1./POW2((double)nLo) -1./POW2((double)ipHi) ) );

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyErg = 
						(float)(EnergyRyd*
					  EN1RYD);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyK = 
						(float)(EnergyRyd* TE1RYD);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN = 
						(float)(EnergyRyd/ WAVNRYD);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].WLAng = 
							(float)(1.0e8/
						  iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN/
						  RefIndex( iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN));

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gLo = 2.L;

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gHi = (float)(2.*ipHi*ipHi);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].Aul = 
						(float)(HydroEinstA(0,ipHi)*z4);

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].iRedisFun = iso.ipResoRedist[ipH_LIKE];

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gf = 
						(float)(GetGF(iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].Aul,
					  iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN,
					  iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gHi));

				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].opacity = 
						(float)(abscf(iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gf,
					  iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN,
					  iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].gLo));

				/* create array index that will blow up */
				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].ipCont = INT_MIN;
			}

			/* set in lines_service.h */
#			if LOWDEN_LYMAN
			/* for fix Lyman lines, in terms of alpha transition */
			for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
			{
					float Ratio_lyman_alpha_Z1[25]={
						0.,0.,0.,
5.52E-01f,/* this mult by 1.33 since Ha split into two */
3.30E-01f,
2.96E-01f,
2.80E-01f,
2.74E-01f,
2.72E-01f,
2.72E-01f,
2.74E-01f,
2.76E-01f,
2.78E-01f,
2.81E-01f,
2.84E-01f,
2.86E-01f,
2.89E-01f,
2.92E-01f,
2.95E-01f,
2.98E-01f,
3.01E-01f,
3.05E-01f,
3.09E-01f,
3.13E-01f,
3.18E-01f};
					float Ratio_lyman_alpha_Z2[25]={
						0.,0.,0.,
4.52E-01f,
2.38E-01f,
1.98E-01f,
1.80E-01f,
1.71E-01f,
1.66E-01f,
1.64E-01f,
1.63E-01f,
1.63E-01f,
1.64E-01f,
1.65E-01f,
1.66E-01f,
1.68E-01f,
1.69E-01f,
1.71E-01f,
1.72E-01f,
1.73E-01f,
1.75E-01f,
1.76E-01f,
1.78E-01f,
1.79E-01f,
1.80E-01f};

						if( nelem==0 )
						{
							EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].Aul = EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].Aul *
							Ratio_lyman_alpha_Z1[MIN2(23,ipHi) ];
						}
						else
						{
							EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].Aul = EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].Aul *
							Ratio_lyman_alpha_Z2[MIN2(23,ipHi) ];
						}

						/* derive the abs coef, call to function is gf, wl (A), g_low */
						EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity = 
							(float)(abscf(
						  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].gf,
						  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].EnergyWN,
						  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].gLo));
			}
#			endif
			/* total toptical depth in 1s - 2s */
			EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].TauTot = 0.;

			/* opacity in two-photon transition is incorrect - actually a continuum,
			 * so divide by typical width*/
			/* >>chng 00 nov 27, do this correction, when considered as line sometimes thick */
			EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].opacity /= 1e4f;

			/* wavelength of 0 is sentinal for two-photon emission */
			EmisLines[ipH_LIKE][nelem][ipH2s][ipH1s].WLAng = 0.;

			/* 2s - 2p cannot have zero opacity */
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].opacity = 1e-30f;

			/* Lya has special redistribution function
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].iRedisFun = ipPRD; */
		}
	}

	/* following comes out very slightly off, correct here */
	EmisLines[ipH_LIKE][ipHELIUM][3][ipH2s].WLAng = 1640.f;
	EmisLines[ipH_LIKE][ipHELIUM][3][ipH2p].WLAng = 1640.f;

	/**************************************************************************
	 * now set HyLife for hydrogen itself, does not include Z^4 or Lyman lines 
	 * HyLife is used to get lifetimes of full H iso sequence, so we must find
	 * largest number of levels that will ever be needed 
	 **************************************************************************/
	need = 0;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			need = MAX2( need , iso.numLevels_max[ipH_LIKE][nelem] );
		}
	}
	need = MAX2(need , iso.nLyman[ipH_LIKE]);

	/* >>chng 01 jan 20, had been La transition prob */
	/* does not include Lyman lines so this is zero,
	 * will be added to line itself below */
	hydro.HyLife[ipH2p] = 0.;

	/* >>chng 00 dec 16, from iso.numLevels_max[ipH_LIKE] for hydrogen, to most ever needed */
	for( ipHi=3; ipHi < need; ipHi++ )
	{
		hydro.HyLife[ipHi] = 0.;

		/* HyLife does not include trans to ground since case b branch used later */
		for( ipLo=ipH2s; ipLo < ipHi; ipLo++ )
		{
			/* array of damping constants for hydrogen itself */
			/* >>chng 00 dec 16, had used hydrogen line structure to get lifetime,
			 * but failed in case where H itself had few levels, and other elements
			 * had more */
			double ea = HydroEinstA(ipLo,ipHi);
			ASSERT(ea>0.);
			hydro.HyLife[ipHi] += (float)ea;
		}
	}

	/* now fill in rest of array with just to make sure never used */
	for( ipHi=need; ipHi<NHYDRO_MAX_LEVEL; ++ipHi )
	{
		hydro.HyLife[ipHi] = -FLT_MAX;
	}

	/* this are not defined and must never be used */
	hydro.HyLife[ipH1s] = -FLT_MAX;
	/* >>chng 01 jan 26, from 8.23 to 8.226 */
	hydro.HyLife[ipH2s] = 8.226f;

	/**************************************************************************
	 *  now fill in full array of damping constants 
	 **************************************************************************/
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			/* charge to 4th power, needed for scaling laws */
			z4 = POW2(nelem+1.);
			z4 *= z4;
			/* put something in for 2s - 2p */
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].dampXvel = 1e-30f;

			for( ipHi=ipH2p; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
			{
				for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
				{
					if( EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN <= 0. )
						continue;
					/* 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 */
					EmisLines[ipH_LIKE][nelem][ipHi][ipLo].dampXvel = 
						(float)(
						/* this is sum of A's to n=2 for H by itself */
						(hydro.HyLife[ipHi]*z4 + 
						/* in following must add on the Lyman lines */
					  EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].Aul)/
					  PI4/EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN);

					ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].dampXvel> 0.);
				}
			}

			/* now do extra Lyman lines */
			ipLo = 0;
			for( ipHi=2; ipHi < iso.nLyman[ipH_LIKE]; ipHi++ )
			{
				iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].dampXvel = 
					(float)(
					/* this is sum of A's to n=2 for H by itself */
					(hydro.HyLife[ipHi]*z4 + 
					/* in following must add on the Lyman lines */
					iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].Aul)/
					PI4/iso.ExtraLymanLines[ipH_LIKE][nelem][ipHi].EnergyWN);
			}
		}
	}

	/* 2p 1s is special since do not assume 2s 2p mixed */
	/* the constant 7.958e-6 in the damping constant is 1/4pi * conversion fractor for microns */
	hdamp[ipH2p] = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul*7.958e-6/0.75;

	/* NB hdamp is only dimensioned up through lmhlvl, so 2s-1s has no damp cnst */
	hydro.HyLife[ipH2p] = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul;

	for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][ipHYDROGEN]; ipHi++ )
	{
		hydro.HyLife[ipHi] = 0.;
		for( ipLo=ipH2s; 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 */
			hydro.HyLife[ipHi] += EmisLines[ipH_LIKE][ipHYDROGEN][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] = (hydro.HyLife[ipHi] + EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipH1s].Aul)*7.958e-6/
		  (1. - 1./POW2( (float)ipHi));
	}

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

	/* define excitation and ionization temperatures for hydrogen only,
	 * helium is done helike where non-hydrogenic energies are assigned */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* define these for H and He always */
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{

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

			for( ipHi=ipH2p; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
			{
				iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipHi] = 
					Potential/POW2((double)ipHi);
				/*HdeltaT.HCionTe[nelem][ipHi] = 
					(float)(TE1RYD*iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipHi]);*/
			}

			/* 2s done here */
			iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH2s] = iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH2p];
			/*HdeltaT.HCionTe[nelem][ipH2s] = 
				(float)(TE1RYD*iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH2p]);*/

			/* 1s done here */
			iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s] = Potential;
			/*HdeltaT.HCionTe[nelem][ipH1s] = 
				(float)(TE1RYD*iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s]);*/
		}
	}

	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
				{
					iso.PopLTE[ipISO][nelem][n] = 0.;
					iso.ColIoniz[ipISO][nelem][n] = 0.;
					iso.gamnc[ipISO][nelem][n] = -DBL_MAX;
					iso.RecomInducRate[ipISO][nelem][n] = -DBL_MAX;
					iso.DepartCoef[ipISO][nelem][n] = -DBL_MAX;
					iso.Pop2Ion[ipISO][nelem][n] = 0.;
					iso.RateLevel2Cont[ipISO][nelem][n] = 0.;
					iso.RateCont2Level[ipISO][nelem][n] = 0.;
					iso.ConOpacRatio[ipISO][nelem][n] = 1.f;
					iso.RadRecomb[ipISO][nelem][n][ipRecRad] = 0.;
					iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] = 1.;
					iso.RadRecomb[ipISO][nelem][n][ipRecEsc] = 1.;
				}
			}
		}
	}

	/* ground state of H and He is different since totally determine
	 * their own opacities */
	iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][0] = 1e-5f;
	iso.ConOpacRatio[ipH_LIKE][ipHELIUM][0] = 1e-5f;
	iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][0] = 1e-5f;


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

	/* zero out some line information */
	iso_zero();

	if( helike.lgHugeCaseB == TRUE )
	{
		iso.numLevels_max[ipHE_LIKE][ipHELIUM] = NumHeLinesSav;
	}

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

/* ============================================================================== */
static void iso_zero(void)
{
	long int i, 
	  ipHi, 
	  ipISO,
	  ipLo, 
	  nelem;

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

	hydro.HLineWidth = 0.;

	/* main isoelectronic arrays, fill with sane values */
	for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			/* >>chng 05 dec 30, move nelem.numLevel to numLevels_local and int here */
			iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];

			/* must always do helium even if turned off */
			if( nelem < 2 || dense.lgElmtOn[nelem] )
			{
				for( ipLo=0; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
				{
					for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
					{
						EmisLines[ipISO][nelem][ipHi][ipLo].Pesc = 1.0;

						/* destruction probability*/
						EmisLines[ipISO][nelem][ipHi][ipLo].Pdest = 0.0;
						EmisLines[ipISO][nelem][ipHi][ipLo].Pelec_esc = 0.0;

						/* upper and lower level populations */
						EmisLines[ipISO][nelem][ipHi][ipLo].PopHi = 0.0;
						EmisLines[ipISO][nelem][ipHi][ipLo].PopLo = 0.0;

						EmisLines[ipISO][nelem][ipHi][ipLo].xIntensity = 0.0;

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

	/* heliumlike autoionization array, fill with sane values */
	for( nelem=1; nelem < LIMELM; nelem++ )
	{
		/* must always do helium even if turned off */
		if( nelem < 2 || dense.lgElmtOn[nelem] )
		{
			for( ipHi = 0; ipHi < 4 ; ipHi++ )
			{
				He1AutoLines[nelem][ipHi].Pesc = 1.0;

				/* destruction probability*/
				He1AutoLines[nelem][ipHi].Pdest = 0.0;
				He1AutoLines[nelem][ipHi].Pelec_esc = 0.0;

				/* upper and lower level populations */
				He1AutoLines[nelem][ipHi].PopHi = 0.0;
				He1AutoLines[nelem][ipHi].PopLo = 0.0;

				He1AutoLines[nelem][ipHi].xIntensity = 0.0;

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


	/* initialize heavy element line array */
	for( i=0; i <= nLevel1; i++ )
	{
		TauLines[i].TauIn = opac.taumin;
		TauLines[i].TauCon = opac.taumin;
		TauLines[i].ColOvTot = 0.;
		/* outward optical depth */
		TauLines[i].TauTot = 1e30f;
		/* escape probability */
		TauLines[i].Pesc = 1.;
		/* inward part of line */
		TauLines[i].FracInwd = 1.;
		/* destruction probability */
		TauLines[i].Pdest = 0.;
		TauLines[i].Pelec_esc = 0.;
		/* line pumping rate */
		TauLines[i].pump = 0.;
		/* population of lower level */
		TauLines[i].PopLo = 0.;
		/* >>chng 97 jul 21, added following zero
		 * population of upper level */
		TauLines[i].PopHi = 0.;
		/* population of lower level with correction for stim emission */
		TauLines[i].PopOpc = 0.;
		/* following two heat exchange excitation, deexcitation */
		TauLines[i].cool = 0.;
		TauLines[i].heat = 0.;
		/* intensity of line */
		TauLines[i].xIntensity = 0.;
		/* number of photons emitted in line */
		TauLines[i].phots = 0.;
	}

	for( i=0; i < nUTA; i++ )
	{
		UTALines[i].PopLo = 0.;
		UTALines[i].PopHi = 0.;

		/* these are line optical depth arrays
		 * inward optical depth */
		UTALines[i].TauIn = opac.taumin;
		UTALines[i].TauCon = opac.taumin;
		UTALines[i].ColOvTot = 0.;
		/* outward optical depth */
		UTALines[i].TauTot = 1e20f;
		/* escape probability */
		UTALines[i].Pesc = 1.;
		/* inward part of line */
		UTALines[i].FracInwd = 1.;
		/* destruction probability */
		UTALines[i].Pdest = 0.;
		UTALines[i].Pelec_esc = 0.;
		/* line pumping rate */
		UTALines[i].pump = 0.;
		/* population of lower level */
		UTALines[i].PopLo = 0.;
		/* population of upper level */
		UTALines[i].PopHi = 0.;
		/* population of lower level with correction for stim emission */
		UTALines[i].PopOpc = 0.;
		/* following two heat exchange excitation, deexcitation */
		UTALines[i].cool = 0.;
		UTALines[i].heat = 0.;
		/* intensity of line */
		UTALines[i].xIntensity = 0.;
		/* number of photons emitted in line */
		UTALines[i].phots = 0.;
	}

	for( i=0; i < nWindLine; i++ )
	{
		TauLine2[i].PopLo = 0.;
		TauLine2[i].PopHi = 0.;

		/* these are line optical depth arrays
		 * inward optical depth */
		TauLine2[i].TauIn = opac.taumin;
		TauLine2[i].TauCon = opac.taumin;
		TauLine2[i].ColOvTot = 0.;
		/* outward optical depth */
		TauLine2[i].TauTot = 1e20f;
		/* escape probability */
		TauLine2[i].Pesc = 1.;
		/* inward part of line */
		TauLine2[i].FracInwd = 1.;
		/* destruction probability */
		TauLine2[i].Pdest = 0.;
		TauLine2[i].Pelec_esc = 0.;
		/* line pumping rate */
		TauLine2[i].pump = 0.;
		/* population of lower level */
		TauLine2[i].PopLo = 0.;
		/* population of upper level */
		TauLine2[i].PopHi = 0.;
		/* population of lower level with correction for stim emission */
		TauLine2[i].PopOpc = 0.;
		/* following two heat exchange excitation, deexcitation */
		TauLine2[i].cool = 0.;
		TauLine2[i].heat = 0.;
		/* intensity of line */
		TauLine2[i].xIntensity = 0.;
		/* number of photons emitted in line */
		TauLine2[i].phots = 0.;
	}

	for( i=0; i < nHFLines; i++ )
	{
		HFLines[i].PopLo = 0.;
		HFLines[i].PopHi = 0.;

		/* these are line optical depth arrays
		 * inward optical depth */
		HFLines[i].TauIn = opac.taumin;
		HFLines[i].TauCon = opac.taumin;
		HFLines[i].ColOvTot = 0.;
		/* outward optical depth */
		HFLines[i].TauTot = 1e20f;
		/* escape probability */
		HFLines[i].Pesc = 1.;
		/* inward part of line */
		HFLines[i].FracInwd = 1.;
		/* destruction probability */
		HFLines[i].Pdest = 0.;
		HFLines[i].Pelec_esc = 0.;
		/* line pumping rate */
		HFLines[i].pump = 0.;
		/* population of lower level */
		HFLines[i].PopLo = 0.;
		/* population of upper level */
		HFLines[i].PopHi = 0.;
		/* population of lower level with correction for stim emission */
		HFLines[i].PopOpc = 0.;
		/* following two heat exchange excitation, deexcitation */
		HFLines[i].cool = 0.;
		HFLines[i].heat = 0.;
		/* intensity of line */
		HFLines[i].xIntensity = 0.;
		/* number of photons emitted in line */
		HFLines[i].phots = 0.;
	}

	for( i=0; i < nCORotate; i++ )
	{
		C12O16Rotate[i].PopLo = 0.;
		C13O16Rotate[i].PopLo = 0.;
		C12O16Rotate[i].PopHi = 0.;
		C13O16Rotate[i].PopHi = 0.;

		/* these are line optical depth arrays
		 * inward optical depth */
		C12O16Rotate[i].TauIn = opac.taumin;
		C13O16Rotate[i].TauIn = opac.taumin;
		C12O16Rotate[i].TauCon = opac.taumin;
		C13O16Rotate[i].TauCon = opac.taumin;
		C12O16Rotate[i].ColOvTot = 0.;
		C13O16Rotate[i].ColOvTot = 0.;
		/* outward optical depth */
		C12O16Rotate[i].TauTot = 1e20f;
		C13O16Rotate[i].TauTot = 1e20f;
		/* escape probability */
		C12O16Rotate[i].Pesc = 1.;
		C13O16Rotate[i].Pesc = 1.;
		/* inward part of line */
		C12O16Rotate[i].FracInwd = 1.;
		C13O16Rotate[i].FracInwd = 1.;
		/* destruction probability */
		C12O16Rotate[i].Pdest = 0.;
		C13O16Rotate[i].Pdest = 0.;
		C12O16Rotate[i].Pelec_esc = 0.;
		C13O16Rotate[i].Pelec_esc = 0.;
		/* line pumping rate */
		C12O16Rotate[i].pump = 0.;
		C13O16Rotate[i].pump = 0.;
		/* population of lower level */
		C12O16Rotate[i].PopLo = 0.;
		C13O16Rotate[i].PopLo = 0.;
		/* population of upper level */
		C12O16Rotate[i].PopHi = 0.;
		C13O16Rotate[i].PopHi = 0.;
		/* population of lower level with correction for stim emission */
		C12O16Rotate[i].PopOpc = 0.;
		C13O16Rotate[i].PopOpc = 0.;
		/* following two heat exchange excitation, deexcitation */
		C12O16Rotate[i].cool = 0.;
		C13O16Rotate[i].cool = 0.;
		C12O16Rotate[i].heat = 0.;
		C13O16Rotate[i].heat = 0.;
		/* intensity of line */
		C12O16Rotate[i].xIntensity = 0.;
		C13O16Rotate[i].xIntensity = 0.;
		/* number of photons emitted in line */
		C12O16Rotate[i].phots = 0.;
		C13O16Rotate[i].phots = 0.;
	}

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

