/* 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 */
/*lines_helium put energetics, H, and He lines into line intensity stack */
/*TempInterp interpolates on a grid of values to produce predicted value at current Te.*/
#include "cddefines.h"
#include "taulines.h"
#include "dense.h"
#include "phycon.h"
#include "lines_service.h"
#include "iso.h"
#include "trace.h"
#include "lines.h"
#include "path.h"
#include "helike.h"
#include "bevington.h"

#define NUMTEMPS	22

typedef struct 
{
	/* index for upper and lower levels of line */
	long int ipHi;
	long int ipLo;

	char label[5];

} stdLines;

static void GetStandardHeLines(void);

static double TempInterp( double* TempArray , double* ValueArray, long NumElements );

static int lgFirstRun = TRUE;
static double CaABTemps[NUMTEMPS];
static long NumLines;
static double ***CaABIntensity;
static stdLines **CaABLines;

void lines_helium(void)
{
	long int i, j, 
	  ipHi, 
	  ipLo, 
	  nelem;
	char chLabel[5]="    ";
	double 
	  /*cb206, 
	  cb5016, 
	  cb5876 ,
	  em, 
	  rec,*/ 
	  sum,
	  Pop2_3S;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "   lines_helium called\n" );
	}

	/* HeI */
	i = StuffComment( "He-like iso-sequence" );
	linadd( 0., (float)i , "####", 'i');

#if	0
	/* He I 4471, */
	if( phycon.te <= 1e4 )
	{
		em = 2.42e-22/(phycon.te/phycon.te10);
	}
	else
	{
		em = 8.79e-22/phycon.te/phycon.te03/phycon.te01;
	}

	/* >>chng 97 jan 04, label was HeI, changed to Ca B for symmetry with others
	 * He I 4471 recombination only, fit to Brocklehurst '72 */
	rec = dense.eden*dense.xIonDense[ipHELIUM][1]*em;
	linadd( rec   , 4471 , "Ca B",'i');

	/* He I 5876 REC, simple fit to Brocklehurst */
	cb5876 = dense.eden*dense.xIonDense[ipHELIUM][1]*4.22e-21/(phycon.te*phycon.te10);
	linadd( cb5876 , 5876 , "Ca B" , 'i' );

	/* this is collisional enhancement of 5876 from 
	 * >>refer	He1	coll	KingdonJ., & Ferland, G.J. 1995, ApJ, 442, 714-725 */
	if( dense.eden < 1e8 )
	{
		double D = 1. + 3131.*100./phycon.sqrte / dense.eden ;
		em = (6.78*0.52581*phycon.te07*sexp(3.776e4/phycon.te) +
			1.67*3.9811/(phycon.te10*phycon.te05)*sexp(4.545e4/phycon.te ) +
			0.60*22.9087/(phycon.te30*phycon.te04)*sexp(4.901e4/phycon.te) ) / D;
	}
	else
	{
		em = 0.;
	}
	/* simple He1 5876 with collisions */
	linadd( cb5876*(1.+em) , 5876 , "wCol" , 'i' );

	/* He I 6678 REC, simple fit to Brocklehurst */
	rec = dense.eden*dense.xIonDense[ipHELIUM][1]*1.32e-21/(phycon.te*phycon.te10*phycon.te01);
	linadd( rec , 6678 , "Ca B" ,'i' );

	em = dense.eden*dense.xIonDense[ipHELIUM][1];
	cb206 = 1.064e-22/phycon.te70/phycon.te10/phycon.te03/phycon.te03*
	  em;

	/* case B He I 2.06 micron - new or old wavelength scaling?  */
	linadd(cb206,20600.,"Ca B",'i');

	/* Case B He I 5016 */
	cb5016 = 5.43e-23/phycon.te70/phycon.te10*em;
	linadd( cb5016 , 5016 , "Ca B", 'i' );
#endif

	/* read in Case A and B lines from data file	*/
	if( lgFirstRun )
	{
		GetStandardHeLines();
		lgFirstRun = FALSE;
	}

	/* this is the main printout, where line intensities are entered into the stack */
	for( nelem=1; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			if( nelem == ipHELIUM )
			{
				double *qTotEff;

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

				qTotEff[0] = 0.;
				qTotEff[1] = 0.;
					
				for( i=ipHe2s3S+1; i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
				{
					qTotEff[i] = 0.;
					for( j = i; j<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; j++ )
					{
						/*if( iso.quant_desig[ipHE_LIKE][nelem][i].s == 1 )
						{*/
							qTotEff[i] += 
								EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].ColUL*dense.EdenHCorr*
								iso.Boltzmann[ipHE_LIKE][nelem][j][ipHe2s3S] *
								(double)EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].gHi / 
								(double)EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].gLo*
								helike.CascadeProb[j][i];
						/*}*/
					}
				}

				/* print simple 2^3S population and the one cloudy got	*/
				/* helike.qTot2TripS already has dense.EdenHCorr included, but we just want eden.*/
				/* chng 05 dec 14, RadRec_caseB now contains only recombinations into triplets,
				 * so no need for the 0.75, which is inaccurate anyway */
				Pop2_3S = dense.eden*iso.RadRec_caseB[ipHE_LIKE][nelem]/
					( EmisLines[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].Aul+ dense.eden*helike.qTot2TripS);

				for( i=0; i< NumLines; i++ )
				{
					ipHi = CaABLines[nelem][i].ipHi;
					ipLo = CaABLines[nelem][i].ipLo;

					if( ipHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem]*(iso.n_HighestResolved_max[ipHE_LIKE][nelem]+1))
					{
						double intens = TempInterp( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
							dense.xIonDense[nelem][nelem]*dense.eden;

						linadd( intens,
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
							CaABLines[nelem][i].label,
							'i');

						if( lgMatch("Ca B",CaABLines[nelem][i].label) )
						{
							/* chng 05 dec 14, branching ratio was missing here! 
							 * not a big effect because lines with biggest collision
							 * enhancements tend to be dominant decay route from upper level. */
							linadd( intens + 
								Pop2_3S*qTotEff[ipHi]*dense.xIonDense[nelem][nelem]*
								helike.BranchRatio[nelem][ipHi][ipLo]*
								EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
								EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
								"+Col",
								'i');
						}
					}
					/* Make sure to at least do 4471 */
					else if( ipLo==ipHe2p3P1 && ipHi==ipHe4d3D && lgMatch("Ca B",CaABLines[nelem][i].label) )
					{
						double intens = TempInterp( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
							dense.xIonDense[nelem][nelem]*dense.eden;

						linadd( intens, 4471, CaABLines[nelem][i].label, 'i');
					}

				}
				free( qTotEff );
			}

			/* NB NB - low and high must be in this order so that all balmer, paschen,
			 * etc series line up correctly in final printout */
			/* >>chng 01 jun 13, bring 23P lines back together */
			/* two photon is special, not a line and no ipCont array index, add here */

			EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].xIntensity = 
				(double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Aul*
				(double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].PopHi*
				(double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Pesc*
				(double)dense.xIonDense[nelem][nelem]*
				(double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg;
			if( LineSave.ipass == 0 )
			{
				/* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" 
				 * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
				/* total two photon emission */
				chIonLbl(chLabel, &EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S]);
			}
			/* two photon continuum */
			linadd( EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].xIntensity , 0,chLabel,'r');

			/* induced two photon emission */
			linadd(
				iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s1S]*
				dense.xIonDense[nelem][nelem]*
				iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*
				EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg,
				22, chLabel ,'i');

			/* here we will create an entry for the three lines 
			 * coming from 2 3P to 1 1S - the entry called TOTL will
			 * appear before the lines of the multiplet */
			sum = 0.;
			for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
			{
				sum += 
					(double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul*
					(double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].PopHi*
					(double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Pesc*
					(double)dense.xIonDense[nelem][nelem]*
					(double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].EnergyErg;
			}

			/* total emission in He-like forbidden lines from 2p3P to ground */
			linadd(sum,EmisLines[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng,"TOTL",'i' );
			
			/* now do real permitted lines */
			for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
			{
				for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
				{
					/* >>chng 01 may 30, do not add fake he-like lines (majority) to line stack */
					/* >>chng 01 dec 11, use variable for smallest A */
					if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 ) 
						continue;

					/* recombine fine-structure lines since the energies are
					 * not resolved anyway.	*/
					if( helike.lgFSM>=1 && ( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l -
						iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)==1 )
						&& (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l>1) 
						&& (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l>1) 
						&& ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n ==
						iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) )
					{
						if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==0) 
							&& (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==0) )
						{
							continue;
						}
						else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==1) 
							&& (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==1) )
						{

							/* singlet to singlet*/
							EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].phots = 
								( (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].Aul*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].PopHi*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].Pesc +
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Aul*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].PopHi*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Pesc )*
								(double)dense.xIonDense[nelem][nelem];
							
							EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].xIntensity = 
								EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].phots *
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].EnergyErg;

							PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1]);
						
							/* triplet to triplet */
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots = 
								( (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc +
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].Aul*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].PopHi*
								(double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].Pesc )*
								(double)dense.xIonDense[nelem][nelem];
							
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity = 
								EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *
								(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;

							PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
						}
					}
					
					else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
					{
						/* here we will create an entry for the three lines 
						 * coming from 2 3P to 2 3S - the entry called TOTL will
						 * appear before the lines of the multiplet 
						 * for He I this is 10830 */

						sum = 0.;
						for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
						{
							sum += 
								(double)EmisLines[ipHE_LIKE][nelem][i][ipLo].Aul*
								(double)EmisLines[ipHE_LIKE][nelem][i][ipLo].PopHi*
								(double)EmisLines[ipHE_LIKE][nelem][i][ipLo].Pesc*
								(double)dense.xIonDense[nelem][nelem]*
								(double)EmisLines[ipHE_LIKE][nelem][i][ipLo].EnergyErg;
						}

						/* total emission in He-like lines */
						linadd(sum,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,"TOTL",'i' );
#						ifdef PRT_HELIKE
#						endif
						
						/*>>chng 05 sep 8, added the following so that the component
						 * from ipHe2p3P0 is printed, in addition to the total. */
						
						/* find number of photons escaping cloud */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots = 
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
							dense.xIonDense[nelem][nelem];

						/* now find line intensity */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity = 
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
							(double)dense.xIonDense[nelem][nelem]*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
						
						if( helike.lgRandErrGen )
						{
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
								helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *= 
								helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
						}
						PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
					}

					else
					{

						/* find number of photons escaping cloud */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots = 
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
							dense.xIonDense[nelem][nelem];

						/* now find line intensity */
						/* >>chng 01 jan 15, put cast double to force double evaluation */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity = 
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
							(double)dense.xIonDense[nelem][nelem]*
							(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;

						if( helike.lgRandErrGen )
						{
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
								helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *= 
								helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
						}

						/* this will enter .xIntensity into the line stack
						fprintf(ioQQQ,"1 loop %li %li %.1f\n", ipLo,ipHi, 
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng ); */
						PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
						{
							/* option to print particulars of some line when called
							 * a prettier print statement is near where chSpin is defined below*/
							/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
							/*@+redef@*/
						if( DEBUG_LOC )
							{
								if( nelem==1 && ipLo==0 && ipHi==1 )
								{
									fprintf(ioQQQ,"he1 626 %.2e %.2e \n", 
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].TauIn,
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].TauTot
										);
								}
							}
						}
					}
				}
			}

			/* this sum will bring together the three lines going to J levels within 23P */
			for( ipHi=ipHe2p3P2+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
			{
				double sumcool , sumheat ,
					save , savecool , saveheat;

				sum = 0;
				sumcool = 0.;
				sumheat = 0.;
				for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
				{
					/* find number of photons escaping cloud */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots = 
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
						dense.xIonDense[nelem][nelem];

					/* now find line intensity */
					/* >>chng 01 jan 15, put cast double to force double evaluation */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity = 
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
						(double)dense.xIonDense[nelem][nelem]*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;

					if( helike.lgRandErrGen )
					{
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
							helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *= 
							helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
					}

					sumcool += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cool;
					sumheat += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].heat;
					sum += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity;
				}

				/* >>chng 01 dec 11, use variable for smallest A */
				/* >>chng 02 feb 10, replaced check on A with ipCont */
				if( EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].ipCont < 1 ) 
					continue;
				/* this will enter .xIntensity into the line stack */
				save = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity;
				savecool = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool;
				saveheat = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat;

				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity = sum;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool = sumcool;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat = sumheat;

				/*fprintf(ioQQQ,"2 loop %li %li %.1f\n", ipHe2p3P2,ipHi, 
					EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].WLAng );*/
				PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2]);

				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity = save;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool = savecool;
				EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat = saveheat;
			}
			for( ipLo=ipHe2p3P2+1; ipLo < iso.numPrintLevels[ipHE_LIKE][nelem]-1; ipLo++ )
			{
				for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
				{
					/* >>chng 01 may 30, do not add fake he-like lines (majority) to line stack */
					/* >>chng 01 dec 11, use variable for smallest A */
					/* >>chng 02 feb 12, use ipCont to find bogus lines */
					/*if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA )*/
					if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 ) 
						continue;

					/* find number of photons escaping cloud */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots = 
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
						dense.xIonDense[nelem][nelem];

					/* now find line intensity */
					/* >>chng 01 jan 15, put cast double to force double evaluation */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity = 
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
						(double)dense.xIonDense[nelem][nelem]*
						(double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;

					if( helike.lgRandErrGen )
					{
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
							helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *= 
							helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
					}

					/* this will enter .xIntensity into the line stack 
					fprintf(ioQQQ,"2 loop %li %li %.1f\n", ipLo,ipHi, 
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng );*/
					PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
				}
			}
		}
	}

	/* ====================================================
	 * end helium
	 * ====================================================*/

#if	0
	/* Need to free these things, but only one time at the end. */
	/* Free some malloc'd arrays */
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		/* TODO	2	- this structure is currently only used for helium itself...
		 * stuff numbers in for other elements, or drop the [nelem] dimension off
		 * of CaABLines	*/
		if( nelem != ipHELIUM )
			continue;

		/* only grab core for elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
		{
			for( j = 0; j < NumLines ; ++j )
			{
				free( CaABIntensity[nelem][j] );
			}
			free( CaABIntensity[nelem] );
			free( CaABLines[nelem] );
		}
	}
	free( CaABIntensity );
	free( CaABLines );
#endif

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "   lines_helium returns\n" );
	}

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


static void GetStandardHeLines()
{
	FILE *ioDATA;
	int lgEOL, lgHIT;
	long i, i1, i2, j, nelem;

#	define chLine_LENGTH 1000
	char chLine[chLine_LENGTH] , 
		/* this must be longer than chDataPath, set in path.h */
		chFilename[FILENAME_PATH_LENGTH_2];

	/* check on path if file not here and path set */
	/* path was parsed in getset */
	if( lgDataPathSet == TRUE )
	{
		/*path set, so look only there */
		strcpy( chFilename , chDataPath );
		strcat( chFilename , "he1_case_ab.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "he1_case_ab.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," lines_helium opening he1_case_ab.dat:");

	if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " lines_helium could not open he1_case_ab.dat\n" );
		if( lgDataPathSet == TRUE )
			fprintf( ioQQQ, " even tried path\n" );

		if( lgDataPathSet == TRUE )
		{
			fprintf( ioQQQ, " lines_helium could not open he1_case_ab.dat\n");
			fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
			fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
		}

		puts( "[Stop in lines_helium]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check that magic number is ok */
	if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
	{
		fprintf( ioQQQ, " lines_helium could not read first line of he1_case_ab.dat.\n");
		puts( "[Stop in lines_helium]" );
		cdEXIT(EXIT_FAILURE);
	}
	i = 1;
	/* first number is magic number, second is number of lines in file	*/
	i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
	NumLines = i2;

	/* the following is to check the numbers that appear at the start of he1_case_ab.dat */
	if( i1 !=CASEABMAGIC )
	{
		fprintf( ioQQQ, 
			" lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
		fprintf( ioQQQ, 
			" lines_helium: I expected to find the number %i and got %li instead.\n" ,
			CASEABMAGIC, i1 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in lines_helium]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* get the array of temperatures */
	lgHIT = FALSE;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		/* only look at lines without '#' in first col */
		if( chLine[0] != '#')
		{
			lgHIT = TRUE;
			j = 0;
			lgEOL = FALSE;
			i = 1;
			while( !lgEOL && j < NUMTEMPS)
			{
				CaABTemps[j] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
				++j;
			}
			break;
		}
	}

	if( !lgHIT )
	{
		fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
		puts( "[Stop in lines_helium]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* create space for array of values, if not already done */
	if( (CaABIntensity = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	/* create space for array of values, if not already done */
	if( (CaABLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		/* TODO	2	- this structure is currently only used for helium itself...
		 * stuff numbers in for other elements, or drop the [nelem] dimension off
		 * of CaABLines	*/
		if( nelem != ipHELIUM )
			continue;

		/* only grab core for elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
		{
			/* create space for array of values, if not already done */
			if( (CaABIntensity[nelem] = (double **)MALLOC(sizeof(double *)*(unsigned)(i2) ) )==NULL )
				BadMalloc();
			if( (CaABLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) ))==NULL )
				BadMalloc();

			/* avoid allocating 0 bytes, some OS return NULL pointer, PvH 
			CaABIntensity[nelem][0] = NULL;*/
			for( j = 0; j < i2 ; ++j )
			{
				if( (CaABIntensity[nelem][j] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS ))==NULL )
					BadMalloc();

				CaABLines[nelem][j].ipHi = -1;
				CaABLines[nelem][j].ipLo = -1;
				strcpy( CaABLines[nelem][j].label , "    " );

				for( i=0; i<NUMTEMPS; ++i )
				{
					CaABIntensity[nelem][j][i] = 0.;
				}
			}
		}
	}

	/* now read in the case A and B line data */
	lgHIT = FALSE;
	nelem = ipHELIUM;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		static long line = 0;
		char *chTemp;
				
		/* only look at lines without '#' in first col */
		if( (chLine[0] == ' ') || (chLine[0]=='\n') )
			break;
		if( chLine[0] != '#')
		{
			lgHIT = TRUE;

			/* get lower and upper level index */
			j = 1;
			/* the first number is the wavelength, which is not used
			 * for anything, but must skip over it. */
			i1 = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			CaABLines[nelem][line].ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			CaABLines[nelem][line].ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			
			ASSERT( CaABLines[nelem][line].ipLo < CaABLines[nelem][line].ipHi );
			/*if( CaABLines[nelem][line].ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] )
				continue;*/

			chTemp = chLine;
			/* skip over 4 tabs to start of data */
			for( i=0; i<3; ++i )
			{
				if( (chTemp = strchr( chTemp, '\t' )) == NULL )
				{
					fprintf( ioQQQ, " lines_helium no init case A and B\n" );
					puts( "[Stop in lines_helium]" );
					cdEXIT(EXIT_FAILURE);
				}
				++chTemp;
			}

			strncpy( CaABLines[nelem][line].label, chTemp , 4 );
			CaABLines[nelem][line].label[4] = 0;
				
			for( i=0; i<NUMTEMPS; ++i )
			{
				float b;
				if( (chTemp = strchr( chTemp, '\t' )) == NULL )
				{
					fprintf( ioQQQ, " lines_helium could not scan case A and B lines, current indices: %li %li\n",
						CaABLines[nelem][line].ipHi,
						CaABLines[nelem][line].ipLo );
					puts( "[Stop in lines_helium]" );
					cdEXIT(EXIT_FAILURE);
				}
				++chTemp;
				sscanf( chTemp , "%e" , &b );
				CaABIntensity[nelem][line][i] = pow(10.,(double)b);
			}
			line++;
		}
	}

	/* close the data file */
	fclose( ioDATA );

	return;
}

static double TempInterp( double* TempArray , double* ValueArray, long NumElements )
{
	/*static long int ipTe=-1;*/
	long int ipTe=-1;
	double rate = 0.;
	double Xarray[4], Yarray[4];
	long i, index[4], numterms = 3;

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

	ASSERT( fabs( 1. - (double)phycon.alogte/log10(phycon.te) ) < 0.0001 );

	/* te totally unknown */
	if( phycon.alogte <= TempArray[0] )
	{
		return ValueArray[0];
	}
	else if( phycon.alogte >= TempArray[NumElements-1] )
	{
		return ValueArray[NumElements-1];
	}

	/* now search for temperature */
	for( i=0; i<NumElements-1; ++i )
	{
		if( ( phycon.alogte > TempArray[i] ) && 
			( phycon.alogte <= TempArray[i+1] ) )
		{
			/* found the temperature - use it */
			ipTe = i;
			break;
		}
	}

	ASSERT( (ipTe >=0) && (ipTe < NumElements-1)  );
			
#if	0
	rate = ((double)phycon.alogte - TempArray[ipTe]) / (TempArray[ipTe+1] - TempArray[ipTe]) *
		(ValueArray[ipTe+1] - ValueArray[ipTe]) + ValueArray[ipTe] ;
#endif
	
	/* Do a four-point interpolation */
	if( ipTe==0 || ipTe==1 )
	{
		index[0]=0;
		index[1]=1;
		index[2]=2;
		index[3]=3;
	}
	else if( ipTe==NumElements-2 )
	{
		index[0]=NumElements-4;
		index[1]=NumElements-3;
		index[2]=NumElements-2;
		index[3]=NumElements-1;
	}
	else
	{
		index[0]=ipTe-1;
		index[1]=ipTe;
		index[2]=ipTe+1;
		index[3]=ipTe+2;
	}

	for( i=0; i<=3; i++ )
	{
		Xarray[i] = TempArray[index[i]];
		Yarray[i] = ValueArray[index[i]];
	}

	interp(Xarray, Yarray, 4, &numterms, phycon.alogte, &rate);

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

	return rate;
}
