/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*prt_He_like_DeparCoef routine to print departure coefficients for he-like species */
/*prt_He_like_Pops routine to print level populations for he-like species */
/*HeCreate create he-like series, called by ContCreatePointers */
/*AGN_He1_CS routine to punch table needed for AGN3 - collision strengths of HeI */
/*defect - calculate quantum defect. */
/*he_energy - calculates energy of a given level.	*/
/*printstuff - in helike.h, printflag determines what will be printed.*/
/*ContinuumLowering - limit max prin. quan. no. due to continuum lowering processes	*/
/*he_assign - assigns quantum numbers and indices to iso.quant_desig and 
 * QuantumNumbers2Index, respectively*/

#include "cddefines.h" 

/* various useful print statements occur throughout this file - search on PRINTIT */
/* 
  
  Energy order within 2 3P

  The order of the levels within the 2 3P level of atomic helium is opposite
  from the order within astrophysically abundant ions.  The indices below
  consistently point to the correct level, and the energies are correct,
  so the J levels within 2 3P are not in increasing energy order for helium itself. 
  This is ok since the atomic data is correct, and the difference in energies is so small.

*/

/* 
the following block of defines occurs in cddefines.c, and echoed in cddefines.h 
const int ipHe1s1S = 0;
const int ipHe2s3S = 1;
const int ipHe2s1S = 2;
const int ipHe2p3P0 = 3;
const int ipHe2p3P1 = 4;
const int ipHe2p3P2 = 5;
const int ipHe2p1P = 6;

level 3
const int ipHe3s3S = 7;
const int ipHe3s1S = 8;
const int ipHe3p3P = 9;
const int ipHe3d3D = 10;
const int ipHe3d1D = 11;
const int ipHe3p1P = 12;

level 4
const int ipHe4s3S = 13;
const int ipHe4s1S = 14;
const int ipHe4p3P = 15;
const int ipHe4d3D = 16;
const int ipHe4d1D = 17;
const int ipHe4f3F = 18;
const int ipHe4f1F = 19;
const int ipHe4p1P = 20;
*/

const int ipRad = 0;
const int ipCollis = 1;
const int ipEnergy = 2;

#include "physconst.h" 
#include "lines_service.h"
#include "tfidle.h"
#include "elementnames.h"
#include "taulines.h"
#include "path.h"
#include "trace.h"
#include "phycon.h"
#include "dense.h"
#include "iso.h"
#include "helike.h"
#include "helike_recom.h"
#include "helike_cs.h"
#include "helike_einsta.h"
#include "hydroeinsta.h"

/*lint -e662 creation of  out of bound pointer */
/*lint -e661 creation of  out of bound pointer */

/*	The Elevels data type, and then the iso.quant_desig[ipHE_LIKE] structure,
 *  which contain the quantum numbers n,l, and s for a given
 *  energy level, are defined in iso.h, for use in multiple subroutines. */

/*typedef struct { 
	long n;
	long s;
	long l;
} Elevels; iso.h */ 

/* An array of structures each containing for a given element, n,l, and s
 * s=0 for singlets, s=1 for triplets
 * in helike.h
static Elevels **iso.quant_desig[ipHE_LIKE]; */

/* the designation of the levels, chLevel[n][string] */
static char **chLevel;

/* compute quantum defect */
static double defect( long nelem, long ipLo );
 
static double he_energy( double Eff_n, long nelem, long ipLo );

/*
static void printstuff(long flag);
*/

static void printCustomAs( void );

/* Function he_assign assigns, in order of increasing energy, 	*/
/* quantum numbers n,l, and s to each element of vector iso.quant_desig[ipHE_LIKE].	*/ 
static void he_assign( long nelem );

static double EthRyd; /*kTRyd,,Xn_S59*/
/*static long ipLev,globalZ;*/
/*ContinuumLowering - limit max prin. quan. no. due to continuum lowering processes	*/
static void ContinuumLowering( long nelem );

static FILE *ofp;  

/* experimental energies, in wavenumbers, for atomic helium */
#define NHE1LEVELS 111
/* These energies from Drake 1996, except these divided by c */
static double He1Energies[NHE1LEVELS] =
{0.0       , 159855.9734, 166277.4390, 169087.8298, 169086.8417, 169086.7652, 171134.8957,
183236.7908, 184864.8281, 185564.6657, 186101.5615, 186104.9656, 186209.3638, 190298.6619,
190940.6075, 191217.0826, 191444.4868, 191446.4547, 191451.8805, 191451.8964, 191492.7108,
193346.9900, 193663.5106, 193800.7280, 193917.1538, 193918.2888, 193921.1207, 193921.1298,
193921.6166, 193921.6209, 193942.4612, 194936.1184, 195114.8672, 195192.7542, 195260.0724,
195260.7694, 195262.4251, 195262.4307, 195262.7236, 195262.7261, 195262.7930, 195262.7947, 
195274.9074, 195868.2357, 195978.8938, 196027.3216, 196069.6730, 196070.1273, 196071.1763,
196071.1800, 196071.3686, 196071.3702, 196071.4141, 196071.4151, 196071.4283, 196071.4290,
196079.0865, 196461.3605, 196534.5628, 196566.7159, 196595.0620, 196595.3730, 196596.0785,
196596.0810, 196596.2092, 196596.2103, 196596.2404, 196596.2411, 196596.2503, 196596.2508,
196596.2541, 196596.2544, 196601.3992, 196861.9861, 196912.9014, 196935.3339, 196955.2261,
196955.4477, 196955.9445, 196955.9463, 196956.0373, 196956.0380, 196956.0595, 196956.0600,
196956.0666, 196956.0670, 196956.0693, 196956.0696, 196956.0705, 196956.0707, 196959.6917,
197145.2320, 197182.0643, 197198.3343, 197212.8252, 197212.9885, 197213.3513, 197213.3527,
197213.4194, 197213.4200, 197213.4358, 197213.4362, 197213.4411, 197213.4414, 197213.4431,
197213.4433, 197213.4440, 197213.4442, 197213.4445, 197213.4446, 197216.0885};
/* Last energy is 10^1P.	*/


/* Ionization potentials (in wavenumber) for each ion in the iso-seq up to Z=30.		*/
/* These are exactly what you get if you take EionRYD below and multiply by RYD_INF.	*/
static double EionWN[29] =       
	/*198307.729760353*/
	{198310.6679     ,610003.839889137,1241136.72201499,2091948.45665631,3162116.52584231,
	 4452446.95015668,5962133.81875305,7692790.05069734,9645221.44709864,11814589.7994457,
	 14209766.0528639,16822685.5022862,19661412.9625169,22717883.6187518,26000162.0663204,
	 29508248.5246975,33234078.1790787,37185715.7345311,41363161.0813172,45766414.4389118,
	 50395475.4781030,55258409.0136949,60339085.8550283,65653635.1927626,71202056.8074231,
	 76976286.4328920,82984388.3352872,89194104.5722390,95726403.3055320};
	
/* Ionization potentials (in Rydbergs) for each ion in the iso-seq up to Z=30.	*/
/* These are exactly what you get if you take EionEV below and divide by EVRYD.	*/
static double EionRYD[29] =       
	/*1.807113*/
	{1.807387521,5.558764,11.310070,19.063237,28.815326,40.573682,54.330961,70.101861,
	87.893725,107.662464,129.488916,153.299590,179.167978,207.020588,236.930910,
	268.898946,302.851204,338.861175,376.928858,417.054255,459.237363,503.551674,
	549.850208,598.279945,648.840883,701.459535,756.209388,812.796486,872.323172};

/* Ionization potentials (in eV) for each ion in the iso-seq up to Z=30. These are 
 * exactly what you get if you take Verner's numbers in PH1COM.PH1[0][1][nelem][0] 
 * and multiply by 0.9998787, exactly as is done elsewhere in Cloudy.	*/
double EionEV[29] =       
	/*24.587017*/
	{24.59077,75.630824,153.881326,259.368529,392.052444,552.033006,739.210311,
	953.784316,1195.854925,1464.822296,1761.786269,2085.746968,2437.704271,
	2816.658298,3223.608929,3658.556163,4120.500123,4610.440686,5128.377852,
	5674.311623,6248.241996,6851.168852,7481.092433,8140.012497,8827.929042,
	9543.842191,10288.751823,11058.658422,11868.560169};

void putError(long nelem, long ipHi, long ipLo, long whichData, float errorToPut)
{
	if( helike.lgRandErrGen )
	{
		ASSERT( whichData <= 2 );
		ASSERT( nelem < LIMELM );
		ASSERT( ipHi <= iso.numLevels[ipHE_LIKE][nelem] + 1 );
		ASSERT( ipLo <= iso.numLevels[ipHE_LIKE][nelem] );

		helike.Error[nelem][ipHi][ipLo][whichData] = errorToPut;
	}

	return;
}

/*He-like main routine to call HeLevel and determine model he-like species atom level balance,
 * called by ionize */
void HeLike(void)
{
	int lowsav , ihisav;
	
	long int ipLo, ipHi, nelem;
	static int lgFinitePop[LIMELM];
	static int lgMustInit=TRUE;

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

	if( lgMustInit )
	{
		for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
			lgFinitePop[nelem] = TRUE;
		lgMustInit = FALSE;
	}

	/* save state of he-like low and high ionziation, since must not change here,
	 * since there is a parallel helium atom that must decrement he */
	/*TODO	2	remove this when this routine really controls helium itself */
	lowsav = dense.IonLow[ipHELIUM];
	ihisav = dense.IonHigh[ipHELIUM];
	for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
	{
		/* do not consider elements that have been turned off */
		if( dense.lgElmtOn[nelem] )
		{
			/* note that nelem scale is totally on c not physical scale, so 0 is h */
			/* evaluate helium-like balance if ionization reaches this high */
			if( (dense.IonHigh[nelem] >= nelem)  )
			{
				ContinuumLowering( nelem );
			
				/* generate gaussian errors */
				if( helike.lgRandErrGen && nzone==0 )
				{
					HeLikeError(nelem);
				}

				/* evaluate recombination rates */
				HeRecom(nelem);

				/* evaluate collisional rates */
				HeCollid(nelem);

				/* >>chng 02 jan 18, move to unified photo rate routine */
				/* evaluate photoionization rates */
				iso_photo(ipHE_LIKE , nelem );

				/* evaluate state specific creation and destruction processes,
				 * also define iso.xIonSimple */
				iso_ionize_recombine( ipHE_LIKE , nelem );

				/* solve for the level populations */
 				HeLikeLevel(nelem);

				/* say that we have set the populations */
				lgFinitePop[nelem] = TRUE;
			}
			else if( lgFinitePop[nelem] )
			{
				/* this branch, pops were set previously, but now are zero,
				 * so must zero them out this time */
				lgFinitePop[nelem] = FALSE;

				iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
				iso.xIonSimple[ipHE_LIKE][nelem] = 0.;

				/* zero it out since no population*/
				iso.Pop2Ion[ipHE_LIKE][nelem][0] = 0.;
				for( ipHi=ipHe2s3S; ipHi < iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
				{
					iso.Pop2Ion[ipHE_LIKE][nelem][ipHi] = 0.;
					for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
					{
						/* population of lower and upper levels rel to ion */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopLo = 0.;
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi =  0.;

						/* population of lower level rel to ion, corrected for stim em */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc =  0.;
					}
				}
			}
		}
	}
	dense.IonLow[ipHELIUM] = lowsav;
	dense.IonHigh[ipHELIUM] = ihisav;

	return;
}

/*defect - calculate quantum defect for a given level and nuclear charge. */
static double defect(long nelem, long ipLo)  
{  
	/* The quantum defect, and parameters a,b, and c	*/
	double qd,a,b,c;	

	double HeDefectAsymptotes[2][10] = {
		{1.40005E-01,-1.20673E-02,2.08056E-03,4.21484E-04,1.14868E-04,
			4.08648E-05,1.73548E-05,8.33891E-06,4.39680E-06,2.42075E-06},
		{2.97063E-01,6.81567E-02,2.82381E-03,4.27703E-04,1.17319E-04,
			4.25254E-05,1.85549E-05,9.24641E-06,5.30882E-06,3.02877E-06}
	};

	/* Parameters for fits to quantum defects for	*/
	/* P triplet and S orbitals. The dimensions are	*/
	/*	first: l									*/
	/*	second: n									*/
	/* 	third: parameters a,b,and c.				*/
	double param[3][4][3]=  
	{ 
		{{0.6451941,0.3119437,-1.2722842},		/* ^3S	*/
			{0.7664874,0.3455675,-1.3976462}, 
			{0.8247101,0.3603131,-1.4520500}, 
			{0.8878402,0.3714450,-1.4995732}}, 

		{{1.4203514,0.5311096,-2.6728087},		/* ^1S	*/
			{1.5733513,0.5997339,-2.9253834}, 
			{1.4531025,0.5924751,-2.8662756}, 
			{1.6038999,0.6342552,-3.0298071}}, 

		{{-2.2323488,0.0890840,-0.5166053},		/* ^3P	*/
			{-2.0463691,0.1222081,-0.6672983}, 
			{-1.9904104,0.1328918,-0.7150879}, 
			{-1.9500974,0.1452111,-0.7649031}} 
	};   

	/* Because they cannot be fit to a funtion of the same form as the other orbitals,
	 * the P singlets are fit to a different function, with these parameters and dimensions */
	/*	first: n									*/
	/*	second: parameters a and b.					*/
	double P1[4][2]=
	{
		{-56.65245,-3.661923},
		{-52.03411,-4.941075},
		{-50.43744,-5.525750},
		{-49.45137,-5.908615}
	};

	long n = iso.quant_desig[ipHE_LIKE][nelem][ipLo].n;
	long lqn = iso.quant_desig[ipHE_LIKE][nelem][ipLo].l;
	long s = iso.quant_desig[ipHE_LIKE][nelem][ipLo].s;
	
	ASSERT(n >= 1L);
	ASSERT(n > lqn);
	/* Only Helium and up, and only those turned on.	*/
	ASSERT((nelem >= ipHELIUM) && (nelem <= LIMELM));       

	if ( n > iso.n_HighestResolved[ipHE_LIKE][nelem] )
		return 0.;

	if ( nelem == ipHELIUM )
	{
		if( ipLo<NHE1LEVELS && n<=iso.n_HighestResolved[ipHE_LIKE][nelem] )
		{
			qd = n-sqrt(1./WAVNRYD/(EionWN[ipHELIUM-1] - He1Energies[ipLo]));
		} 
		else if ( lqn<=9 )
		{
			qd = HeDefectAsymptotes[s][lqn];
		}
		else if ( s == 0 )
		{
			qd = 0.0497*pow((double)lqn, -4.4303);
		}
		else
		{
			qd = 0.0656*pow((double)lqn, -4.5606);
		}

		return qd;
	}
	else if ( ipLo == ipHe1s1S )
	{
		/* Quantum defects for ground state are found from the rydberg
		 * equation, and the ionization potential of the ion.	*/
		ASSERT(nelem>ipHYDROGEN );
		qd = 1.0 - nelem * sqrt(1/EionRYD[nelem-1]);
		return qd;  
	}
	else
	{
		/* For levels with n > 5, the quantum defect	*/
		/* is approximately the same as if n equaled 5.	*/
		if ( n > 5L )
		{
			n = 5L;
		}
		/* For P singlets	*/
		if ( lqn==1L && s==0L )
		{
			qd = 1./(P1[n-2][0] + P1[n-2][1] * (nelem+1) * log((double)nelem+1.) );
			return qd;  
		}
		/* Defects for orbitals with l>2 are approximately equal to zero.	*/
		else if ( lqn < 2L )
		{
			a = param[2*(lqn+1)-s-1][n-2][0];  
			b = param[2*(lqn+1)-s-1][n-2][1];  
			c = param[2*(lqn+1)-s-1][n-2][2];  
			qd = exp((a+c*(nelem+1))/(1.0+b*(nelem+1)));  
			return qd;  
		}
		/* This fit is a simplification of table 11.9 from 
		 * >>refer Helike	defects	Drake, G.W.F., editor.  Atomic, Molecular & Optical Physics Handbook.
		 * >>refercon	Chapter 11, "High Precision Calculations for Helium", G.W.F. Drake. 
		 * >>refercon	AIP Press: Woodbury, New York, 1996
		 * This will give quasi-real energies for all transitions, allowing a reasonable
		 * determination of which decays are zeroed due to being below the plasma frequency.
		 * The 1/nelem dependence is arbitray.  	*/
		else
		{
			ASSERT( lqn >= 2L );
			return ( ( 0.0612/(double)nelem ) / pow((double)lqn, 4.44) );
		}
	}   
}

/*he_energy calculates energy of a given level.	*/
static double he_energy(double Eff_n,long nelem, long ipLo )
{
	double Ef,gap2tripP;
	
	/* for atomic helium use experimental energies */
	if( nelem == ipHELIUM && ipLo < NHE1LEVELS )
	{
		Ef = EionWN[ipHELIUM-1] - He1Energies[ipLo];
	} 
	else
	/* Calculate energy difference for all other trans.	*/
	{
		Ef = RYD_INF * (nelem/Eff_n) * (nelem/Eff_n);  
	}
	
	if (nelem > ipHELIUM)
	{
		/* gap2tripP represents the amount to add or subtract from the median 
		 * 2p3P term in order to find energies of individual levels of the term.	*/
		/* >>refer	Helike	triplets	Porquet and Dubau (2000)	*/
		/* TODO	2	This is not the best source of data.  Replace it.	*/
		gap2tripP = 0.007*pow((double)nelem+1., 4.5);
		if ( ipLo == ipHe2p3P0 )
		{
			Ef += (3*gap2tripP);
		}
		else if ( ipLo == ipHe2p3P1 )
		{
			Ef += (2*gap2tripP);
		}
		else if ( ipLo == ipHe2p3P2 )
		{
			Ef -= (5*gap2tripP);
		}
	}

	ASSERT(Ef > 0.);
	return Ef;
}

/* Function he_assign assigns, in order of increasing energy, 				*/
/* quantum numbers n,l, and s to each element of vector iso.quant_desig[ipHE_LIKE].			*/
static void he_assign( long nelem)
{
	/*	in, il, and is are dummies for n, l, and s.		*/
	/*	i is the index in the vector iso.quant_desig[ipHE_LIKE].			*/
	long in, il, is, ij, i = 0 , level;
	
	/* this loop is over quantum number n */
	for ( in = 1L; in <= iso.n_HighestResolved[ipHE_LIKE][nelem] ; ++in )
	{
		for ( il = 0L; il < in; ++il )
		{
			for ( is = 1L; is >= 0 ; --is )
			{
				/* All levels except singlet P follow the ordering scheme:	*/
				/*	lower l's have lower energy	*/
				/* 	triplets have lower energy	*/
				if ( (il == 1L) && (is == 0L) )
					continue;
				/* n = 1 has no triplet, of course.	*/
				if ( (in == 1L) && (is == 1L) )
					continue;
				ij = 0;
				do 
				{
					iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
					iso.quant_desig[ipHE_LIKE][nelem][i].s = is;
					iso.quant_desig[ipHE_LIKE][nelem][i].l = il;
					QuantumNumbers2Index[nelem][in][il][is] = i;
					++i;
					++ij;
				}	while ( (in == 2) && (il == 1) && (is == 1) && (ij < 3) );
			}
		}
		/*	Insert singlet P at the end of every sequence for a given n.	*/
		if ( in > 1L )
		{
			iso.quant_desig[ipHE_LIKE][nelem][i].n = in;
			iso.quant_desig[ipHE_LIKE][nelem][i].s = 0L;
			iso.quant_desig[ipHE_LIKE][nelem][i].l = 1L;
			QuantumNumbers2Index[nelem][in][1][0] = i;
			++i;
		}
	}
	/* now do the collapsed levels */
	in = iso.n_HighestResolved[ipHE_LIKE][nelem] + 1;
	for( level = i; level< iso.numLevels[ipHE_LIKE][nelem]; ++level)
	{
		iso.quant_desig[ipHE_LIKE][nelem][level].n = in;
		iso.quant_desig[ipHE_LIKE][nelem][level].s = -LONG_MAX;
		iso.quant_desig[ipHE_LIKE][nelem][level].l = -LONG_MAX;
		/* Point every l and s to same index for collapsed levels.	*/
		for( il = 0; il < in; ++il )
		{
			for( is = 0; is <= 1; ++is )
			{
				QuantumNumbers2Index[nelem][in][il][is] = level;
			}
		}
		++in;
	}
	--in;

	/* confirm that we did not overrun the array */
	ASSERT( i <= iso.numLevels[ipHE_LIKE][nelem] );

	ASSERT( (in > 0) && (in < (iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] + 1) ) );


	/* Verify iso.quant_desig[ipHE_LIKE] and QuantumNumbers2Index agree in all cases	*/
	for ( in = 2L; in <= iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem]  ; ++in )
	{
		for ( il = 0L; il < in; ++il )
		{
			for ( is = 1L; is >= 0 ; --is )
			{
				/* Ground state is not triplicate.	*/
				if ( (in == 1L) && (is == 1L) )
					continue;
				
				ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].n == in );
				if( in <= iso.n_HighestResolved[ipHE_LIKE][nelem] )
				{
					/* Must only check these for resolved levels...
					 * collapsed levels in iso.quant_desig[ipHE_LIKE] have pointers for l and s that will blow if used.	*/
					ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].l == il );
					ASSERT( iso.quant_desig[ipHE_LIKE][nelem][ QuantumNumbers2Index[nelem][in][il][is] ].s == is );
				}
			}
		}
	}
	
	return;
}

/*printstuff printflag determines what will be printed.*/
/*
static void printstuff(long flag)
{
	long int i,ic,nelem,ipLo,ipHi;
	int lgHeader = 0;
	double thresh_energy;

	switch(flag)
	{
	case 1:
		if ( (ofp = fopen("HelikeAs.txt","w")) == NULL )
			return;
		fprintf(ofp," Z\tni^JL\tnf^JL\tipLo\tipHi\tA\thydro\tscqdri\n");
		break;
	case 2:
		if ( (ofp = fopen("HelikefAs.txt","w")) ==NULL )
			return;
		fprintf(ofp,"FORBIDDEN	TRANSITIONS	to 1^1S\n");
		fprintf(ofp,"Z\t2^3S\t2^1S\t2^3P1\n");
		break;
	case 3:
		if ( (ofp = fopen("HelikeEs.txt","w")) == NULL )
			return;
		break;
	case 4:
		if ( ( ofp = fopen("HelikeThPCS.txt","w")) == NULL )
			return;
		break;
	case 5:
		if ( ( ofp = fopen("HelikeRR.txt","w")) ==NULL )
			return;
		break;
	case 6:
		if ( ( ofp = fopen("HelikePCSgrid.txt","w")) == NULL )
			return;
		break;
	default:
		return;
	}
	
	for (i=1; i<N_HE_TE_RECOMB; i++)
	{
		if (flag == 5)
		{
			fprintf(ofp,"temp= %2.1f \n",0.5*i);
		}
		else if (lgHeader)
		{
			return;
		}
		if ( (flag > 2) && (flag < 6) && (!lgHeader) )
		{
			fprintf(ofp,"WARNING: These levels are as resolved for He...not all elements necessarily resolved in same way!\n");
			fprintf(ofp,"Z   ");
			for( ipLo=ipHe1s1S; ipLo < iso.numLevels[ipHE_LIKE][ipHELIUM] - iso.nCollapsed[ipHE_LIKE][ipHELIUM]; ipLo++ )
			{
				fprintf(ofp,"n%1lds%1ldl%1ld      ",
					iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].n,
					iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].s,
					iso.quant_desig[ipHE_LIKE][ipHELIUM][ipLo].l);
			}
			fprintf(ofp,"\n");
		}
		lgHeader = TRUE;

		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 
			{
				if ( (flag > 1) && (flag < 6) )
				{
					fprintf(ofp,"%2ld ",nelem+1L);
				}
				
				globalZ = nelem;
				
				for( ipLo=ipHe1s1S; ipLo < MIN2(iso.numLevels[ipHE_LIKE][ipHELIUM] - iso.nCollapsed[ipHE_LIKE][ipHELIUM],
					iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem]); ipLo++ )
				{
					if (ipLo == ipHe1s1S)
					{
						thresh_energy = EionRYD[nelem-1L];
					}
					else
					{
						thresh_energy = EionRYD[nelem-1L] - EmisLines[ipHE_LIKE][nelem][ipLo][ipHe1s1S].EnergyWN * WAVNRYD;
					}
					
					thresh_energy = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo];
					ipLev = ipLo;
					EthRyd = thresh_energy;
					
					switch(flag)
					{
					case 1:
						for( ipHi=ipLo+1L; ipHi < iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem]; ipHi++)
						{
							double RI2, Aul;
							if ( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul>iso.SmallA &&
								abs(L_(ipHi)-L_(ipLo))==1 )
							{
								fprintf(ofp,"%2ld\t%2ld^%1ld",nelem+1L, N_(ipLo), 2*(S_(ipLo))+1 );
								print spectroscopic label here
								fprintf(ofp,"\t%2ld^%1ld", N_(ipHi), 2*(S_(ipHi))+1 );
								print spectroscopic label here
								fprintf(ofp,"\t%3ld\t%3ld\t%2.5e",ipLo,ipHi,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul);
								if( N_(ipHi)>N_(ipLo)  )
								{
									fprintf(ofp,"\t%e", H_Einstein_A(N_(ipHi) ,L_(ipHi) ,N_(ipLo) ,L_(ipLo) , nelem) );
								}
								else
								{
									fprintf(ofp,"\t%e", -1.0);
								}
								RI2 = scqdri(N_(ipHi) - defect(nelem,ipHi),
									L_(ipHi),
									N_(ipLo) - defect(nelem,ipLo),
									L_(ipLo),(double)(nelem)
									);
								Aul = ritoa(L_(ipHi),L_(ipLo),nelem,
									(iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLo]-iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi])*RYD_INF,RI2);
								if ( (ipLo == ipHe2p3P0) || (ipLo == ipHe2p3P1) || (ipLo == ipHe2p3P2) )
								{
									Aul *= (2.*(ipLo-3)+1.0)/9.0;
								}

								fprintf(ofp,"\t%e", Aul );
								fprintf(ofp,"\t%e", EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng );
								if( N_(ipHi)==N_(ipLo) )
									fprintf( ofp,"\t%e\n",
										EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*helike.Lifetime[ipHELIUM][ipHi]);
								else
									fprintf( ofp,"\n");

									
							}
						}
						break;
					case 2:
						if( (ipLo == ipHe2s3S) || (ipLo == ipHe2s1S) || (ipLo == ipHe2p3P1) )
						{
							fprintf(ofp,"\t%2.2e",EmisLines[ipHE_LIKE][nelem][ipLo][0].Aul);
						}
						break;
					case 3:
						fprintf(ofp,"%2.4e ",thresh_energy);
						break;
					case 4:
						fprintf(ofp,"%2.4e ",He_cross_section(thresh_energy, ipLo, nelem)); 
						break;
					case 5:
						fprintf(ofp,"%2.4e ",pow(10.,RRCoef[nelem][ipLo][i]));
						break;
					case 6:
						fprintf(ofp,"Z=%ld  %2ld^%1ld",nelem+1,	N_(ipLo), 2*(S_(ipLo))+1 );
						print spectroscopic label here
						fprintf(ofp,"\n");
						for (ic=0; ic<100; ic++)
						{
							fprintf(ofp,"%2.2e %2.2e\n",thresh_energy*(1+0.1*pow((double)ic,1.5)),
									He_cross_section(thresh_energy*(1+0.1*pow((double)ic,1.5)), ipLo, nelem));
						}
						break;
					default:
						fprintf(ofp,"\n");
					}
				}
				fprintf(ofp,"\n");
			}
		}
	}
	fclose(ofp);
	return;
}
*/

void He1Autoionize(void)
{

/* Autoionizing levels are as follows:	*/
	const int ipHe2s2p1P = 0;
/*	const int ipHe2s2p3P = 1;	*/
/*	const int ipHe2p2_1D = 2;	*/
 	const int ipHe2p2_3P = 3;		

	double a[4] = {75827., 72116., 74853., 74010.};
	double b[4] = {1.8956, 1.9084,    1.9, 1.9021};

	long nelem, ipHi;

	for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
	{
		/* only set elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			for ( ipHi = ipHe2s2p1P; ipHi <= ipHe2p2_3P; ipHi++ )
			{
				He1AutoLines[nelem][ipHi].EnergyWN = (float)(a[ipHi]*pow((double)(nelem+1),b[ipHi]));
			}
		}
	}

	return;
	/* TODO	2	say where these come from...and do something with them!	*/	
	/* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002	*/
	/* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001.	*/
}

#undef DEBUG_LOC 

/* create he-like series, called by ContCreatePointers */
void HeCreate(void)
{
	double **energies, *n_effs, **SumAPerN, **RobbinsP, **RobbinsC;

	long int i, i1;
	long int j, ipLo, ipHi, ipFirstCollapsed, nelem;

	static int nCalled = 0;

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

	FILE *ioDATA;
	int lgEOL, lgHIT;

	char chSpin[2]={'1','3'};
	char chL[6]={'S','P','D','F','G','H'};

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

	/* now MALLOC the space for the helium-like lines, 
	 * but this must be done exactly one time per coreload */
	if( nCalled > 0 )
	{
#		ifdef DEBUG_FUN
		fputs( " <->HeCreate()\n", debug_fp );
#		endif
		return;
	}
	++nCalled;

	if( helike.lgSetBenjamin )
	{
		max_num_levels = iso.numLevels[ipHE_LIKE][ipHELIUM];
		max_n = iso.n_HighestResolved[ipHE_LIKE][ipHELIUM];
		iso.lgInd2nu_On = FALSE;
		
		if( helike.lgCompileRecomb )
		{
			fprintf( ioQQQ, " Not compiling He recomb file...can not do with Benjamin command set.");
			helike.lgCompileRecomb = FALSE;
		}
		if( helike.lgFSM )
		{
			fprintf( ioQQQ, " Can not include f-s mixing with the Benjamin command!  I have turned it off.");
			helike.lgFSM = FALSE;
		}
	}
	else
	{
		/* first find the largest number of levels in any element in the he-like sequence */
		max_num_levels = 0;
		max_n = 0;
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			/* only check elements that are turned on */
			if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
			{
				max_num_levels = MAX2( max_num_levels, iso.numLevels[ipHE_LIKE][nelem] );
				max_n = MAX2( max_n, iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] );
			}
		}
	}

	if( helike.lgRandErrGen )
	{
		if( (helike.Error = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();

		if( (helike.ErrorFactor = (float****)MALLOC(sizeof(float***)*(unsigned)LIMELM ) )==NULL )
			BadMalloc();
	}

	if( (helike.n_HighestResolved = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	
	if( (helike.nCollapsed = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();
	
	if( (helike.numLevels = (long*)MALLOC(sizeof(long)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/*char ***helike chLevel[lev][lev][str]*/
	if( (chLevel = (char **)MALLOC(sizeof(char *)*(unsigned)(max_num_levels) ) )==NULL )
		BadMalloc();

	/* In this array are stored the P values described in Robbins 68.	*/
	if( (RobbinsP = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();

	/* In this array are stored the C values described in Robbins 68.	*/
	if( (RobbinsC = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();

	/* In this array are stored the C values described in Robbins 68.	
	 * this is a permanent one for helium only.	*/
	if( (helike.TotBranchRatio = (double**)MALLOC(sizeof(double*)*(unsigned)(max_num_levels) ))==NULL )
		BadMalloc();

	/* Wavenumber level energies for each element are stored here.	*/
	if( (energies = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* The lifetime of individual levels.  helike.Lifetime[LIMELM][iso.numLevels[ipHE_LIKE][nelem]]	*/
	if( (helike.Lifetime = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* The sum of all A's coming out of a given n...used to make sure trend is monotonic	*/
	if( (SumAPerN = (double**)MALLOC(sizeof(double*)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	/* the reverse of iso.quant_desig[ipHE_LIKE]...QuantumNumbers2Index[nelem][n][l][s] 
	 * gives corresponding index in iso.quant_desig[ipHE_LIKE].	*/
	if( (QuantumNumbers2Index = (long****)MALLOC(sizeof(double***)*(unsigned)LIMELM ) )==NULL )
		BadMalloc();

	for( ipLo=ipHe1s1S; ipLo < iso.numLevels[ipHE_LIKE][ipHELIUM] ;++ipLo )
	{
		if( (chLevel[ipLo] = (char*)MALLOC(sizeof(char)*(unsigned)10 ))==NULL )
			BadMalloc();
	}

	/* >>refer He triplets	Robbins, R.R. 1968, ApJ 151, 497R	*/
	/* >>refer He triplets	Robbins, R.R. 1968a, ApJ 151, 511R	*/
	for( ipLo=ipHe1s1S; ipLo < max_num_levels ;++ipLo )
	{
		if( (RobbinsP[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
			BadMalloc();

		if( (RobbinsC[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
    		BadMalloc();

		if( (helike.TotBranchRatio[ipLo] = (double*)MALLOC(sizeof(double)*(unsigned)(ipLo+1) ))==NULL )
			BadMalloc();
	}
	
	for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
	{
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			if( (energies[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]) ))==NULL )
				BadMalloc();
			
			if( helike.lgRandErrGen )
			{
				if( (helike.Error[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]+1) ) )==NULL )
					BadMalloc();

				if( (helike.ErrorFactor[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]+1) ) )==NULL )
					BadMalloc();
			}

			for( i = 0; i<= iso.numLevels[ipHE_LIKE][nelem]+1; i++ )
			{
				if( helike.lgRandErrGen )
				{
					if( (helike.Error[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]+1) ) )==NULL )
						BadMalloc();

					if( (helike.ErrorFactor[nelem][i] = (float**)MALLOC(sizeof(float*)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]+1) ) )==NULL )
						BadMalloc();

					for( j = 0; j<= iso.numLevels[ipHE_LIKE][nelem]; j++ )
					{
						if( (helike.Error[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
							BadMalloc();

						if( (helike.ErrorFactor[nelem][i][j] = (float*)MALLOC(sizeof(float)*(unsigned)(3) ) )==NULL )
							BadMalloc();

						/* set each of these to negative one */
						for( i1=0; i1<3; i1++ )
						{
							helike.Error[nelem][i][j][i1] = -1.f;
							helike.ErrorFactor[nelem][i][j][i1] = -1.f;
						}
					}
				}
			}

			if( (helike.Lifetime[nelem] = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]) ))==NULL )
				BadMalloc();

			/* The sum of all A's coming out of a given n...used to make sure trend is monotonic	*/
			if( (SumAPerN[nelem] = (double*)MALLOC(sizeof(double)*
				(unsigned)(iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] + 1) ) )==NULL )
				BadMalloc();
			
			/* Allocate proper number of principal quantum number.	*/
			if( (QuantumNumbers2Index[nelem] = (long***)MALLOC(sizeof(double**)*
				(unsigned)(iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] + 1) ) )==NULL )
				BadMalloc();

			for( i = 0; i <= (iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem]); ++i)
			{
				/* Allocate proper number of angular momentum quantum number.	*/
				if( (QuantumNumbers2Index[nelem][i] = (long**)MALLOC(sizeof(double*)*(unsigned)MAX2(1,i) ) )==NULL )
					BadMalloc();

				for( i1 = 0; i1 < i; ++i1 )
				{
					/* Allocate 2 spaces for spin.	*/
					if( (QuantumNumbers2Index[nelem][i][i1] = (long*)MALLOC(sizeof(double)*(unsigned)2 ) )==NULL )
						BadMalloc();
				}
			}
		}
	}

	for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
	{
		/* only check elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			/* Function he_assign assigns, in order of increasing energy, 
			 * quantum numbers n,l,s, and j to each element of vector iso.quant_desig[ipHE_LIKE]. */
			he_assign(nelem);
	    }
	}

	/* >>allocate one more than needed to allow space for end of vector sentinal
	 * in case where no collapsed levels */
	if ( (n_effs= (double *)MALLOC((unsigned)(max_num_levels+1)*sizeof(double)) ) == NULL )
		BadMalloc();
	{
		/* 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 )
		{
			for( nelem = ipHELIUM; nelem <LIMELM; ++nelem )
			{
				/* first print internal representation */
				for( i=0; i < iso.numLevels[ipHE_LIKE][nelem]; ++i )
				{
					fprintf(ioQQQ,"nelem %li i %li n %li s %li l %li\n", 
						nelem, i, N_(i), S_(i), L_(i) );
				}
				if( DEBUG_LOC )
					cdEXIT(EXIT_SUCCESS);
			}
		}
	}

	/* This routine reads in transition probabilities from a file. */ 
	HelikeTransProbSetup();

	/* main helium-like arrays, fill with sane values */
	for( nelem = ipHELIUM; nelem < LIMELM; nelem++ )
	{
		double Atot;
		/* charge to 4th power, needed for scaling laws for As*/
		double z4 = POW2((double)nelem);
		z4 *= z4;

		/* must always do helium even if turned off */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] ) 
		{
			for( i=0; i <= ( iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] ); ++i)
			{
				SumAPerN[nelem][i] = 0.;
			}

			/* this is the number of resolved levels, so first collapsed level is [ipFirstCollapsed] */
			ipFirstCollapsed = iso.numLevels[ipHE_LIKE][nelem]-iso.nCollapsed[ipHE_LIKE][nelem]; 
			
			/* Initialize some ground state stuff, easier here than in loops.	*/
			RobbinsC[ipHe1s1S][ipHe1s1S] = 1.;
			RobbinsP[ipHe1s1S][ipHe1s1S] = 0.;
			helike.Lifetime[nelem][ipHe1s1S] = 0.;
			helike.TotBranchRatio[ipHe1s1S][ipHe1s1S] = 1.;
			iso.stat[ipHE_LIKE][nelem][ipHe1s1S] = 1.f;
			n_effs[ipHe1s1S] = N_(ipHe1s1S) - defect(nelem,ipHe1s1S);
			energies[nelem][ipHe1s1S] = EionWN[nelem-1];
			iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] = energies[nelem][ipHe1s1S]*WAVNRYD;
			ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHe1s1S] > 0. );
			
			/*************************************************************/
			/**  Statistical weights, Line and level energies, and		**/
			/**  transition probabilities are calculated in this loop.	**/
			/*************************************************************/
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels[ipHE_LIKE][nelem] ; ipHi++ )
			{
				if (ipHi >= ipFirstCollapsed )
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = 4.f*POW2((float)N_(ipHi));
				}
				else if ( ipHi>=ipHe2p3P0 && ipHi<=ipHe2p3P2 )
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*J_(ipHi) + 1.f);
				}
				else
				{
					iso.stat[ipHE_LIKE][nelem][ipHi] = (2.f*L_(ipHi) + 1.f) * (2.f*S_(ipHi) + 1.f);
				}

				for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
				{
					double Enerwn, Aul, Aul1;
					
					/* atomic number or charge and stage:			  */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].nelem = (int)(nelem+1);
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].IonStg = (int)(nelem);

					/* store statistical weights in EmisLines */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo = iso.stat[ipHE_LIKE][nelem][ipLo];
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi = iso.stat[ipHE_LIKE][nelem][ipHi];
					ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi> 0.);
					ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo> 0.);
					
					/************************************************/
					/**********  Level and line energies ************/
					/************************************************/
					{
						/* the effective quantum number */
						n_effs[ipHi] = N_(ipHi) - defect(nelem,ipHi);
						
						/* energies (in wavenumbers) that correspond to quantum defect */
						/* for atomic helium use experimental energies */
						if( nelem==ipHELIUM && ipHi<NHE1LEVELS &&
							N_(ipHi)<=iso.n_HighestResolved[ipHE_LIKE][nelem] )
						{
							energies[nelem][ipHi] = EionWN[ipHELIUM-1] - He1Energies[ipHi];
							n_effs[ipHi] = sqrt(1./energies[nelem][ipHi]/WAVNRYD);
						} 
						else if( N_(ipHi) > iso.n_HighestResolved[ipHE_LIKE][nelem] )
							energies[nelem][ipHi] = RYD_INF * POW2((double)nelem/(double)N_(ipHi));  
						else
							energies[nelem][ipHi] = he_energy(n_effs[ipHi],nelem,ipHi);
						
						/* convert wavenumbers to Ryd.	*/
						EthRyd = energies[nelem][ipHi] * WAVNRYD;
						
						/* Store threshold energy in larger cloudy array	*/
						iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] = EthRyd;
						ASSERT( iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipHi] > 0. );

						Enerwn = energies[nelem][ipLo] - energies[nelem][ipHi];

						if ( Enerwn == 0 )
							Enerwn = 0.0001;

						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN = (float)Enerwn;

						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg = (float)(Enerwn*WAVNRYD*EN1RYD);

						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK = (float)(Enerwn*WAVNRYD*TE1RYD );
					}
					
					/************************************************/
					/**********  Transition probabilities ***********/
					/************************************************/
					{
						if ( ipHi >= ipFirstCollapsed )
						{
							if ( ipLo >= ipFirstCollapsed )
							{
								/* Neither upper nor lower is resolved...Aul's are easy.	*/
								Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
								putError(nelem,ipHi,ipLo,ipRad,0.001f);

								ASSERT( Aul > 0.);
							}
							else 
							{
								/* Lower level resolved, upper not. First calculate Aul
								 * from upper level with ang mom one higher.	*/
								Aul = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)+1,
									S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
								Aul *= (2.*L_(ipLo)+3.) * (2.*S_(ipLo)+1.) / 
									(4.*(double)N_(ipHi)*(double)N_(ipHi));
								if ( L_(ipLo) != 0 )
								{
									/* For all l>0, add in transitions from upper level with ang mom one lower.	*/
									Aul1 = he_1trans( nelem, Enerwn, n_effs[ipHi], L_(ipLo)-1,
										S_(ipLo), -1, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
									/* now add in this part after multiplying by stat weight for lHi = lLo-1.	*/
									Aul += Aul1*(2.*L_(ipLo)-1.) * (2.*S_(ipLo)+1.) /
										(4.*(double)N_(ipHi)*(double)N_(ipHi));
								}
								putError(nelem,ipHi,ipLo,ipRad,0.01f);
								ASSERT( Aul > 0.);
							}
							if( N_(ipHi) > N_(ipLo) )
								SumAPerN[nelem][N_(ipHi)] += Aul;
						}
						else
						{
							/* Both levels are resolved...do the normal bit.	*/
							Aul = he_1trans( nelem, Enerwn,  n_effs[ipHi],
								L_(ipHi), S_(ipHi), ipHi-3, n_effs[ipLo], L_(ipLo), S_(ipLo), ipLo-3);
							if( N_(ipHi) > N_(ipLo) )
							{
								SumAPerN[nelem][N_(ipHi)] += Aul*(2.*L_(ipHi) + 1.)*(2.*S_(ipHi) + 1.)/
									(4. * (double)N_(ipHi) * (double)N_(ipHi)) ;
							}
						}
						/* set the transition probability */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul = (float)(Aul);

						/*if( N_(ipHi)==N_(ipLo) && N_(ipHi)>2 )
						{
							EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul = iso.SmallA;
						}*/
							
						ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul > 0.);
					}
				}
			}

			/************************************************/
			/**********  Fine Structure Mixing - FSM ********/
			/************************************************/
			if( helike.lgFSM == TRUE )
			{
				for( ipHi=ipHe2s3S; ipHi<iso.numLevels[ipHE_LIKE][nelem] ; ipHi++ )
				{
					for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
					{
						DoFSMixing( nelem, ipLo, ipHi );
					}
				}
			}

			/************************************************/
			/**********  stuff gf, opacity, etc.   **********/
			/************************************************/
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels[ipHE_LIKE][nelem] ; ipHi++ )
			{
				Atot = 0;
				RobbinsC[ipHi][ipHi] = 1.;
				RobbinsP[ipHi][ipHi] = 0.;
				helike.Lifetime[nelem][ipHi] = 0.;
				helike.TotBranchRatio[ipHi][ipHi] = 1.;

				for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
				{
					RobbinsC[ipHi][ipLo] = 0.;
					RobbinsP[ipHi][ipLo] = 0.;
					helike.TotBranchRatio[ipHi][ipLo] = 0.;

					/* there are some negative energy transitions, where the order
					 * has changed, but these are not optically allowed, these are
					 * same n, different L, forbidden transitions */
					ASSERT( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN > 0. ||
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA );
					
					/* very forbidden transitions may have negative energies */
					if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul <= iso.SmallA ||
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN <= 0.)
					{
						/* this branch - bogus line */
						/* transition energy in various units:*/

						/* take abs value since some have order switched */
						/* make following an air wavelength */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng = 1e6f;

						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf = 1e-20f;

						/* derive the abs coef, call to function is gf, wl (A), g_low */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 0.;
					}
					else
					{
						/* take abs value since some have order switched */
						/* make following an air wavelength */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng = 
							(float)fabs(1.0e8/
						  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/
						  RefIndex(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]));
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf = 
							(float)(GetGF(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul,
						  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
						  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi));
						ASSERT(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf > 0.);

						/* derive the abs coef, call to function is gf, wl (A), g_low */
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 
							(float)(abscf(EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gf,
						  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
						  EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo));
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity = 
							MAX2(0.f, EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].opacity );
					}
					{
						/* option to print particulars of some line when called
						 * a prettier print statement is near where chSpin is defined below
						 * search for "pretty print" */
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC  )
						{
							if( (nelem == ipHELIUM) /*&& (EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul> 1e-19)*/ )
							{
								fprintf(ioQQQ,"Z %li Lo %li n %li s %li l %li \t ", 
									nelem+1, ipLo, N_(ipLo), S_(ipLo), L_(ipLo) );
								fprintf(ioQQQ," Hi %li n %li s %li l %li \t", 
									ipHi, N_(ipHi), S_(ipHi), L_(ipHi) );
								fprintf(ioQQQ,"%.4e\t%.4e\tcs\t%.4e\n",
									/*EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,*/
									1.e8/EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN,
									EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul ,
									EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cs);
							}
						}
					}

					/* redistribution function = partial or complete? */
					if( ipLo == ipHe1s1S )
					{
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipPRD;
					}
					else
					{
						/* >>chng 01 feb 27, had been -1, crd with core only,
						 * change to crd with wings as per discussion with Ivan Hubeny */
						/*EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = -1;*/
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].iRedisFun = ipCRDW;
					}
					/* create array index that will blow up */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont = INT_MIN;
					iso.Pop2Ion[ipHE_LIKE][nelem][ipLo] = 0.;

					/* This is not actually the lifetime...must invert when all A's added.	*/
					/* No stat. weight is needed because this is the lifetime of an individual level,
					 * Collapsed or resolved is irrelevant.	*/
					helike.Lifetime[nelem][ipHi] += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul;
				}

				/* convert sum of As into proper lifetime in seconds */
				helike.Lifetime[nelem][ipHi] = 1./helike.Lifetime[nelem][ipHi];

				if ( ipHi < ipFirstCollapsed )
				{
					for( ipLo=ipHe1s1S; ipLo<ipHi ; ipLo++ )
					{
						if ( (Atot > 0.) && (S_(ipLo) == S_(ipHi)) && 
							( abs((int)(L_(ipLo) - L_(ipHi))) == 1 ) )
							RobbinsP[ipHi][ipLo] = (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul/Atot;
						else 
							RobbinsP[ipHi][ipLo] = 0.;
					}
					
					for( ipLo=ipHe1s1S; ipLo<=ipHi ; ipLo++ )
					{
						for ( i=ipLo /*+ 1*/; i<=ipHi; i++ )
						{
							RobbinsC[ipHi][ipLo] += RobbinsP[ipHi][i] * RobbinsC[i][ipLo];
						}
						if( nelem == ipHELIUM )
							helike.TotBranchRatio[ipHi][ipLo] = RobbinsC[ipHi][ipLo];
					}
				}
			}

			/**********************************************************************/
			/** Allowed decay conversion probabilities. See Robbins68b, Table 1. **/
			/** These are not used by the code.  They're simply calculated for   **/
			/** output, (so Ryan can make unbelieveably pretty tables.)			 **/
			/**********************************************************************/
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
					
				if( DEBUG_LOC && (nelem == ipHELIUM ) )
				{
					/* To output Bm(n,l; ipLo), set ipLo, lo_l, and lo_s accordingly.	*/
					long hi_l,hi_s;
					double Bm;
					
					/* tripS to 2^3P	
					hi_l = 0;
					hi_s = 1;
					ipLo = ipHe2p3P0;*/

					/* tripD to 2^3P	
					hi_l = 2;
					hi_s = 1;
					ipLo = ipHe2p3P0;*/

					/* tripP to 2^3S	*/
					hi_l = 1;
					hi_s = 1;
					ipLo = ipHe2s3S;	
					
					ASSERT( hi_l != iso.quant_desig[ipHE_LIKE][nelem][ipLo].l );
					
					fprintf(ioQQQ,"Bm(n,%ld,%ld;%ld)\n",hi_l,hi_s,ipLo);
					fprintf(ioQQQ,"m\t2\t\t3\t\t4\t\t5\t\t6\n");

					for( ipHi=ipHe2p3P2; ipHi<iso.numLevels[ipHE_LIKE][nelem]-iso.nCollapsed[ipHE_LIKE][nelem] ; ipHi++ )
					{
						/* Pick out excitations from metastable 2tripS to ntripP.	*/
						if ( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l == 1) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s == 1) )
						{
							fprintf(ioQQQ,"\n%ld\t",iso.quant_desig[ipHE_LIKE][nelem][ipHi].n);
							j = 0;
							Bm = 0;
							for ( i = ipLo; i<=ipHi; i++)
							{
								if ( (iso.quant_desig[ipHE_LIKE][nelem][i].l == hi_l) && (iso.quant_desig[ipHE_LIKE][nelem][i].s == hi_s)  )
								{
									if ( (ipLo == ipHe2p3P0) && (i > ipHe2p3P2) )
									{
										Bm += RobbinsC[ipHi][i] * ( RobbinsP[i][ipHe2p3P0] + 
											RobbinsP[i][ipHe2p3P1] + RobbinsP[i][ipHe2p3P2] );
									}
									else
										Bm += RobbinsC[ipHi][i] * RobbinsP[i][ipLo];
										
									if ( (i == ipHe2p3P0) || (i == ipHe2p3P1) || (i == ipHe2p3P2) )
									{
										j++;
										if (j == 3)
										{
											fprintf(ioQQQ,"%2.4e\t",Bm);
											Bm = 0;
										}
									}
									else
									{
										fprintf(ioQQQ,"%2.4e\t",Bm);
										Bm = 0;
									}
								}
							}
						}
					}
					fprintf(ioQQQ,"\n\n");
				}
			}

			/**************************************/
			/*********  Extra Lyman lines    ******/
			/**************************************/
			ipLo = ipHe1s1S;
			/* This ipHi does NOT represent ipHe2s1S, but the n in n^1P. */
			for( ipHi=2; ipHi <iso.nLyman[ipHE_LIKE]; ipHi++ )
			{
				double Enerwn , Aul /*,RI2,npstar */ ;
				float gHi , gLo;

				/* pointers to transitions, ipLo, ipHi
					* actual quantum numbers, nLo nHi
					* atomic number or charge and stage: */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].nelem = (int)(nelem+1);
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].IonStg = (int)(nelem);

				/* stat weight of all n ^1P is 3 */
				gHi = 3.f;

				gLo = 1.f;

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo = gLo;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gHi = gHi;

				/* hydrogenic energy, first ryd, convert to wn */
				Enerwn = EionWN[nelem-1] - RYD_INF*POW2((double)nelem)/POW2((double)ipHi);
				
				if ( nelem == ipHELIUM )
				{
					/* A simple fit for the calculation of Helium lyman Aul's.	*/
					Aul = (1.508e10) / pow((double)ipHi,2.975);
				}
				else 
				{
					/* Fit to values given in 
					 * >>refer	He-like	As	Johnson, W.R., Savukov, I.M., Safronova, U.I., & 
					 * >>refercon	Dalgarno, A., 2002, ApJS 141, 543J	*/
					/* originally astro.ph. 0201454  */
					Aul = 1.375E10 * pow((double)nelem, 3.9) / pow((double)ipHi,3.1);

				}

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul = (float)(Aul);

				/* the energy of the transition in wavenumbers, this can be negative
					* for highly forbidden upper levels where energy switches */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN = (float)Enerwn;

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyErg = 
					(float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
					EN1RYD);

				/* the damping constant times the velocity */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].dampXvel = (float)(Aul/
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/PI4);

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyK = 
					(float)(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN * WAVNRYD*
					TE1RYD );

				/* transition energy in various units:*/

				/* take abs value since some have order switched */
				/* make following an air wavelength */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].WLAng = 
					(float)fabs(1.0e8/
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN/
					RefIndex(&iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi]));

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

				/* derive the abs coef, call to function is gf, wl (A), g_low */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity = 
					(float)(abscf(iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gf,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].EnergyWN,
					iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].gLo));

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

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pesc = 1.0;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pelec_esc = 0.0;

				/* destruction probability*/
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Pdest = 0.0;

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

				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].xIntensity = 0.0;

				/* following three will eventually be set by model atoms, but
					* must have reasonable value at start */
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].pump = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].cool = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].heat = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].ColOvTot = 0.;
				iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].PopOpc = 0.;
				{
					/* option to print particulars of some line when called
						* a prettier print statement is near where chSpin is defined below
						* search for "pretty print" */
					/*@-redef@*/
					enum {DEBUG_LOC=FALSE};
					/*@+redef@*/
					if( DEBUG_LOC )
					{
						fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",
							nelem+1,
							ipHi,
							iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].Aul , 	
							iso.ExtraLymanLines[ipHE_LIKE][nelem][ipHi].opacity
							);
					}
				}
			}

			/* This is above the test below so that, in the case that one of those asserts blows,
			 * we will still have this complete set of A's for diagnostics. */
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC && (nelem == ipHELIUM) )
				{
					printCustomAs();
				}
			}

			/******************************************************/
			/***  Lifetimes should increase monotonically with  ***/
			/***  increasing n...Make sure the A's decrease.	***/
			/******************************************************/
			{
				for( i=2; i < iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] - 1; ++i)
				{
					ASSERT( SumAPerN[nelem][i] > SumAPerN[nelem][i+1] );
				}
			}
			
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC /* && (nelem == ipHELIUM) */)
				{
					for( i = 2; i<= (iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem]); ++i)
					{
						fprintf(ioQQQ,"n %ld\t lifetime %.4e\n", i, 1./SumAPerN[nelem][i]);
					}
				}
			}

			/* this is the two photon transition, designated by wl of 0 */
			EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].WLAng = 0.;

			/* opacity in two-photon continuum is not correct since treated as line,
			 * TODO this is part of the induce 2nu problem, which must be fixed for both
			 * he and h sequences. */
			EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].opacity /= 1e4f;
			
			/* we now have lifetimes for everything - set up true damping constants */
			for( ipHi=ipHe2s3S; ipHi<iso.numLevels[ipHE_LIKE][nelem]; ipHi++ )
			{
				for( ipLo=ipHe1s1S; ipLo<ipHi; ipLo++ )
				{
					/* the damping constant times the velocity */
					EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].dampXvel = (float)(
						1./helike.Lifetime[nelem][ipHi]/
						EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN/PI4);
				}
			}
		}
	}

	/* This routine reads in precompiled recombination coefficients for an
	 * array of temperatures.  Interpolation is later necessary to arrive at 
	 * coefficients for a given temperature. */
	HelikeRecombSetup();
	
/*	printstuff(printflag); */

	/* free the memory	*/
	free(n_effs );

	/* create spectroscopic designation of labels */
	for( ipLo=ipHe1s1S; ipLo < iso.numLevels[ipHE_LIKE][ipHELIUM]-iso.nCollapsed[ipHE_LIKE][ipHELIUM]; ++ipLo )
	{
		nelem = ipHELIUM;
		sprintf( chLevel[ipLo] , "%li %c%c", N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))] );
	}

	/* option to create trace printout */
	if( (trace.lgTrace && trace.lgHeBug) )
	{

		/* "pretty print" only helium, O, and Fe 
		 * comments elsewhere say to search for "pretty print" */
		/* next print designation of levels */
		fprintf(ioQQQ,"These indices may not represent the same levels for each element!\n");
		fprintf(ioQQQ,"index\tdesig\tE(wn,He)\tE(wn,O)\tE(wn,Fe)\n");
		for( ipHi=1; ipHi<( MAX3(iso.numLevels[ipHE_LIKE][ipHELIUM],
			iso.numLevels[ipHE_LIKE][ipIRON],iso.numLevels[ipHE_LIKE][ipOXYGEN]) - iso.nCollapsed[ipHE_LIKE][ipOXYGEN] ); ++ipHi )
		{
			float enHe=0. , enO=0., enFe=0.;
			nelem = -1;
			if( dense.lgElmtOn[ipHELIUM] &&(ipHi < iso.numLevels[ipHE_LIKE][ipHELIUM]) )
			{
				enHe = EmisLines[ipHE_LIKE][ipHELIUM][ipHi][0].EnergyWN;
				nelem = ipHELIUM;
			}
			else
			{
				enHe = 0.;
			}
			if( dense.lgElmtOn[ipOXYGEN] &&(ipHi < iso.numLevels[ipHE_LIKE][ipOXYGEN]) )
			{
				enO = EmisLines[ipHE_LIKE][ipOXYGEN][ipHi][0].EnergyWN;
				nelem = ipOXYGEN;
			}
			else
			{
				enO = 0.;
			}
			if( dense.lgElmtOn[ipIRON] &&(ipHi < iso.numLevels[ipHE_LIKE][ipIRON]) )
			{
				enFe = EmisLines[ipHE_LIKE][ipIRON][ipHi][0].EnergyWN;
				nelem = ipIRON;
			}
			else
			{
				enFe = 0.;
			}
			ASSERT( nelem >0 && nelem < LIMELM );
			if( ipHi < iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] )
			{
				fprintf(ioQQQ,"%li\t%li %c%c\t%.5e\t%.5e\t%.5e\n", 
					ipHi, N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))], 
					enHe, enO, enFe  );
			}
		
		}

		fprintf(ioQQQ,"lo\thi\t He states \n");
		nelem = ipHELIUM;
		for( ipHi=ipHe2s3S; ipHi<iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] ; ++ipHi )
		{
			for( ipLo=ipHe1s1S; ipLo<ipHi ; ++ipLo )
			{
				fprintf(ioQQQ," %li\t%li\t%li %c%c\t%li %c%c\n",
					ipLo, ipHi, 
					N_(ipLo), chSpin[S_(ipLo)], chL[MIN2(5,L_(ipLo))], 
					N_(ipHi), chSpin[S_(ipHi)], chL[MIN2(5,L_(ipHi))] );
			}
		}
	}

	if( iso.numLevels[ipHE_LIKE][ipHELIUM] > ipHe2p1P )
	{
		/* some wavelengths of HeI are expected to be exact */
		EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P0][ipHe2s3S].WLAng = 10830.;
		EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P1][ipHe2s3S].WLAng = 10830.;
		EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p3P2][ipHe2s3S].WLAng = 10830.;
		EmisLines[ipHE_LIKE][ipHELIUM][ipHe2p1P][ipHe2s1S].WLAng = 20600.;

		if(  iso.numLevels[ipHE_LIKE][ipHELIUM] > ipHe3d1D )
		{
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].WLAng = 7065.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].WLAng = 7065.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].WLAng = 7065.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].WLAng = 3889.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].WLAng = 5876.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].WLAng = 5876.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].WLAng = 5876.;
			EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d1D][ipHe2p1P].WLAng = 6678.;

			if(  iso.numLevels[ipHE_LIKE][ipHELIUM] > ipHe4d3D )
			{
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P0].WLAng = 4471.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P1].WLAng = 4471.;
				EmisLines[ipHE_LIKE][ipHELIUM][ipHe4d3D][ipHe2p3P2].WLAng = 4471.;
			}
		}
	}

	{
		/* prints table of wavelengths of lines to ground */
		/*@-redef@*/
		enum {PRINTIT=FALSE};
		/*@+redef@*/
		if( PRINTIT )
		{
			for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
			{
				if( dense.lgElmtOn[nelem] )
				{
					fprintf(ioQQQ,"%li\t%s\t",nelem+1,elementnames.chElementSym[nelem] );

					ipLo = ipHe1s1S;
					for( ipHi=ipLo+1; ipHi<=ipHe2p1P; ++ipHi )
					{
						if( ipHi==ipHe2s1S || ipHi==ipHe2p3P0 ) continue;
						fprintf(ioQQQ,"%.4f\t", EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng);
					}
					fprintf(ioQQQ,"\n");
				}
			}
			cdEXIT(EXIT_SUCCESS);
		}
	}

	/* get the collision strength data for the He 1 lines */

	/* 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_cs.dat" );
	}
	else
	{
		/* path not set, check local space only */
		strcpy( chFilename , "he1_cs.dat" );
	}

	if( trace.lgTrace )
		fprintf( ioQQQ," HeCreate opening he1_cs.dat:");

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

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

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

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

	/* the following is to check the numbers that appear at the start of he1_cs.dat */
	if( i1 !=COLLISMAGIC )
	{
		fprintf( ioQQQ, 
			" HeCreate: the version of he1_cs.dat is not the current version.\n" );
		fprintf( ioQQQ, 
			" HeCreate: I expected to find the number %i and got %li instead.\n" ,
			COLLISMAGIC, i1 );
		fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
		puts( "[Stop in HeCreate]" );
		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;
			helike.nCS = 0;
			lgEOL = FALSE;
			i = 1;
			while( !lgEOL && helike.nCS < HE1CSARRAY)
			{
				helike.CSTemp[helike.nCS] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
				++helike.nCS;
			}
			break;
		}
	}
	--helike.nCS;
	if( !lgHIT )
	{
		fprintf( ioQQQ, " HeCreate could not find line in CS temperatures in he1_cs.dat.\n");
		puts( "[Stop in HeCreate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* create space for array of CS values, if not already done */
	if( (helike.HeCS = (float ****)MALLOC(sizeof(float ***)*(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 helike.HeCS	*/
		if( nelem != ipHELIUM )
			continue;

		/* only grab core for elements that are turned on */
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
		{
			if( (helike.HeCS[nelem] = (float***)MALLOC(sizeof(float**)*(unsigned)(iso.numLevels[ipHE_LIKE][nelem]- iso.nCollapsed[ipHE_LIKE][nelem]) ))==NULL )
				BadMalloc();

			/* avoid allocating 0 bytes, some OS return NULL pointer, PvH */
			helike.HeCS[nelem][0] = NULL;
			for( ipHi=ipHe2s3S; ipHi < iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem];++ipHi )
			{
				if( (helike.HeCS[nelem][ipHi] = (float**)MALLOC(sizeof(float*)*(unsigned)ipHi ))==NULL )
					BadMalloc();

				for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo )
				{
					if( (helike.HeCS[nelem][ipHi][ipLo] = (float*)MALLOC(sizeof(float)*(unsigned)helike.nCS ))==NULL )
						BadMalloc();

					for( i=0; i<helike.nCS; ++i )
					{
						helike.HeCS[nelem][ipHi][ipLo][i] = 0.f;
					}
				}
			}
		}
	}

	/* now read in the CS data */
	lgHIT = FALSE;
	nelem = ipHELIUM;
	while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
	{
		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;
			ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
			ASSERT( ipLo < ipHi );
			if( ipHi >= iso.numLevels[ipHE_LIKE][nelem] - iso.nCollapsed[ipHE_LIKE][nelem] )
				continue;
			else
			{
				chTemp = chLine;
				/* skip over 4 tabs to start of cs data */
				for( i=0; i<3; ++i )
				{
					if( (chTemp = strchr( chTemp, '\t' )) == NULL )
					{
						fprintf( ioQQQ, " HeCreate no init cs\n" );
						puts( "[Stop in HeCreate]" );
						cdEXIT(EXIT_FAILURE);
					}
					++chTemp;
				}

				for( i=0; i<helike.nCS; ++i )
				{
					float a;
					if( (chTemp = strchr( chTemp, '\t' )) == NULL )
					{
						fprintf( ioQQQ, " HeCreate not scan cs, current indices: %li %li\n", ipHi, ipLo );
						puts( "[Stop in HeCreate]" );
						cdEXIT(EXIT_FAILURE);
					}
					++chTemp;
					sscanf( chTemp , "%e" , &a );
					helike.HeCS[nelem][ipHi][ipLo][i] = a;
				}
			}
		}
	}

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

	/* option to print cs data for AGN */
	/* create spectroscopic designation of labels */
	{
		/* option to print particulars of some line when called
		 * a prettier print statement is near where chSpin is defined below*/
		/*@-redef@*/
		enum {AGN=FALSE};
		/*@+redef@*/
		if( AGN )
		{
#			define NTEMP 6
			double te[NTEMP]={6000.,8000.,10000.,15000.,20000.,25000. };
			double telog[NTEMP] ,
				cs ,
				ratecoef;
			nelem = ipHELIUM;
			fprintf( ioQQQ,"trans");
			for(i=0; i<NTEMP; ++i)
			{
				telog[i] = log10( te[i] );
				fprintf( ioQQQ,"\t%.3e",te[i]);
			}
			for(i=0; i<NTEMP; ++i)
			{
				fprintf( ioQQQ,"\t%.3e",te[i]);
			}
			fprintf(ioQQQ,"\n");

			for( ipHi=ipHe2s3S; ipHi< iso.numLevels[ipHE_LIKE][ipHELIUM]; ++ipHi )
			{
				for( ipLo=ipHe1s1S; ipLo < ipHi; ++ipLo )
				{
					/* print the designations of the lower and upper levels */
					fprintf( ioQQQ,"%s - %s",
						chLevel[ipLo] , chLevel[ipHi] );

					/* print the interpolated collision strengths */
					for(i=0; i<NTEMP; ++i)
					{
						phycon.alogte = (float)telog[i];
						/* print cs */
						cs = HeCSInterp( nelem , ipHi , ipLo );
						fprintf(ioQQQ,"\t%.2e", cs );
					}

					/* print the rate coefficients */
					for(i=0; i<NTEMP; ++i)
					{
						phycon.alogte = (float)telog[i];
						phycon.te = (float)pow(10.,telog[i] );
						tfidle(FALSE);
						cs = HeCSInterp( nelem , ipHi , ipLo );
						/* collisional deexcitation rate */
						ratecoef = cs/sqrt(phycon.te)*COLL_CONST/iso.stat[ipHE_LIKE][nelem][ipLo] *
							sexp( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyK / phycon.te );
						fprintf(ioQQQ,"\t%.2e", ratecoef );
					}
					fprintf(ioQQQ,"\n");
				}
			}
			cdEXIT(EXIT_FAILURE);
		}
	}

	for( ipLo=ipHe1s1S; ipLo < max_num_levels ;++ipLo )
	{
		free( RobbinsC[ipLo] );
		free( RobbinsP[ipLo] );
	}
	for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
	{
		if( nelem == ipHELIUM || dense.lgElmtOn[nelem] )
		{
			free( energies[nelem] );
			free( SumAPerN[nelem] );
		}
	}
	free( energies );
	free( SumAPerN );
	
	free( RobbinsC );
	free( RobbinsP );

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

/* routine to punch table needed for AGN3 - collision strengths of HeI */
void AGN_He1_CS( FILE *ioPun )
{

	long int i;

	/* list of temperatures where cs will be printed */
#	define NTE 5 
	double TeList[NTE] = {6000.,10000.,15000.,20000.,25000.};
	float TempSave;

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

	/* put on a header */
	fprintf(ioPun, "Te\t2 3s 33s\n");

	/* Restore the original temp when this routine done.	*/
	TempSave = phycon.te;

	for( i=0; i<NTE; ++i )
	{
		phycon.te = (float)TeList[i];
		tfidle(FALSE);

		fprintf(ioPun , "%.0f\t", 
			TeList[i] );
		fprintf(ioPun , "%.2f\t", 
			HeCSInterp( 1 , ipHe3s3S , ipHe2s3S ) );
		fprintf(ioPun , "%.2f\t", 
			HeCSInterp( 1 , ipHe3p3P , ipHe2s3S ) );
		fprintf(ioPun , "%.2f\t", 
			HeCSInterp( 1 , ipHe3d3D , ipHe2s3S ) );
		fprintf(ioPun , "%.3f\t", 
			HeCSInterp( 1 , ipHe3d1D , ipHe2s3S ) );
		/*fprintf(ioPun , "%.1f\t%.1f\t%.1f\n", */
		fprintf(ioPun , "%.1f\n", 
			HeCSInterp( 1 , ipHe2p3P0 , ipHe2s3S ) +
			HeCSInterp( 1 , ipHe2p3P1 , ipHe2s3S ) +
			HeCSInterp( 1 , ipHe2p3P2 , ipHe2s3S ));
	}

	phycon.te = TempSave;
	/* no need to force update since didn't do above	*/
	tfidle(FALSE);

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

	return;
#	undef NTE
}

static void ContinuumLowering( long nelem )
{
	double a;
	long np, nd, ns, nc;
	
	/* size of rate matrices will be defined according to the n calculated here	*/
	
	ASSERT( nelem < LIMELM );

	/* Particle packing - the density here should be density of all nuclei in the plasma */
	a = sqrt( 1.8887E8 * (nelem+1.) / pow(dense.xNucleiTotal, 0.333) );
	if( a > (double)iso.n_HighestResolved[ipHE_LIKE][nelem]+(double)iso.nCollapsed[ipHE_LIKE][nelem] )
	{
		np = iso.n_HighestResolved[ipHE_LIKE][nelem]+iso.nCollapsed[ipHE_LIKE][nelem] + 1;
	}
	else
		np = (long)a;

	/* Debye shielding - the density here is electron density	*/
	a = 2.6E7 * nelem * nelem * pow( phycon.te/dense.eden, 0.25);
	if( a > (double)iso.n_HighestResolved[ipHE_LIKE][nelem]+(double)iso.nCollapsed[ipHE_LIKE][nelem] )
	{
		nd = iso.n_HighestResolved[ipHE_LIKE][nelem]+iso.nCollapsed[ipHE_LIKE][nelem] + 1;
	}
	else
		nd = (long)a;
	
	/* Stark broadening - this should be the density of singly charged ions, 
	 * both positive and negative.  The sum of protons, electrons, and HeII should be
	 * good enough.	*/
	a = 3171. * pow( (double)nelem, 0.8 ) * pow( dense.eden + (double)dense.xIonDense[ipHYDROGEN][1]
		+ (double)dense.xIonDense[ipHELIUM][1], -0.1333);
	if( a > (double)iso.n_HighestResolved[ipHE_LIKE][nelem]+(double)iso.nCollapsed[ipHE_LIKE][nelem] )
	{
		ns = iso.n_HighestResolved[ipHE_LIKE][nelem]+iso.nCollapsed[ipHE_LIKE][nelem] + 1;
	}
	else
		ns = (long)a;
	
	ASSERT( np > 3 );
	ASSERT( nd > 3 );
	ASSERT( ns > 3 );

	nc = MIN3(np, nd, ns);
	
	/* I assert greater than three because the code depends upon having at least up to n=3, and
	 * because it would take EXTREMELY dense conditions to lower the continuum that much, and
	 * the code would probably get a wrong answer then anyway.  */
	ASSERT( nc > 3 );

	if( nc < iso.n_HighestResolved[ipHE_LIKE][nelem])
	{
		helike.n_HighestResolved[nelem] = nc;
		helike.nCollapsed[nelem] = 0;
		helike.numLevels[nelem] = nc*nc + nc + 1;
		helike.lgFullSize[nelem] = TRUE;
		iso.numPrintLevels[ipHE_LIKE][nelem] = helike.numLevels[nelem];
	}
	else if( nc <= iso.n_HighestResolved[ipHE_LIKE][nelem] + iso.nCollapsed[ipHE_LIKE][nelem] )
	{
		helike.n_HighestResolved[nelem] = iso.n_HighestResolved[ipHE_LIKE][nelem];
		helike.nCollapsed[nelem] = nc - helike.n_HighestResolved[nelem];
		helike.numLevels[nelem] = helike.n_HighestResolved[nelem]*helike.n_HighestResolved[nelem] +
			helike.n_HighestResolved[nelem] + 1 + helike.nCollapsed[nelem];
		helike.lgFullSize[nelem] = TRUE;
		iso.numPrintLevels[ipHE_LIKE][nelem] = helike.numLevels[nelem] - helike.nCollapsed[nelem];
	}
	else
	{
		helike.numLevels[nelem] = iso.numLevels[ipHE_LIKE][nelem];
		helike.nCollapsed[nelem] = iso.nCollapsed[ipHE_LIKE][nelem];
		helike.n_HighestResolved[nelem] = iso.n_HighestResolved[ipHE_LIKE][nelem];
		helike.lgFullSize[nelem] = FALSE;
	}

	/* This 100 is to be kept parallel with the initialization in zero.c	*/
	/* TODO fix this to go to much larger n */
	iso.nLyman[ipHE_LIKE] = MIN2( nc, 100);

	ASSERT( helike.numLevels[nelem] <= iso.numLevels[ipHE_LIKE][nelem] );
	ASSERT( helike.nCollapsed[nelem] <= iso.nCollapsed[ipHE_LIKE][nelem] );
	ASSERT( helike.n_HighestResolved[nelem] <= iso.n_HighestResolved[ipHE_LIKE][nelem] );
	
	return;
}

static void printCustomAs( void )
{
	long i, ipLo, familyIndex, nHi, upperIndex, nelem;
	nelem = ipHELIUM;

	if ( (ofp = fopen("CustomAs.txt","w")) == NULL )
		return;
	fprintf(ofp,"ipLo\tnLo\tlLo\tsLo\tlHi\t" );
	for( i = 2; i<= iso.n_HighestResolved[ipHE_LIKE][ipHELIUM]; i++ )
	{
		fprintf( ofp,"%li\t", i);
	}
	fprintf( ofp, "\n" );
		
	for( ipLo = 0; ipLo< iso.numLevels[ipHE_LIKE][ipHELIUM] - iso.nCollapsed[ipHE_LIKE][ipHELIUM] - 1; ipLo++ )
	{
		if( N_(ipLo) < iso.n_HighestResolved[ipHE_LIKE][ipHELIUM] )
		{
			for( familyIndex = N_(ipLo)*(N_(ipLo)+1)+1; familyIndex < (N_(ipLo)+1)*(N_(ipLo)+2)+1; familyIndex++ )
			{
				if( N_(familyIndex) <= iso.n_HighestResolved[ipHE_LIKE][ipHELIUM] )
				{
					if( ( abs(L_(ipLo)-L_(familyIndex)) == 1 ) && (S_(ipLo) == S_(familyIndex)) )
					{
						fprintf( ofp, "%li\t%li\t%li\t%li\t%li\t",
							ipLo, N_(ipLo), L_(ipLo), S_(ipLo), L_(familyIndex) );

						for( nHi = 2; nHi < N_(familyIndex); nHi++ )
						{
							fprintf( ofp, "\t" );
						}
			
						for( nHi = N_(familyIndex); nHi <= iso.n_HighestResolved[ipHE_LIKE][ipHELIUM]; nHi++ )
						{
							upperIndex = QuantumNumbers2Index[ipHELIUM][nHi][L_(familyIndex)][S_(familyIndex)];

							fprintf( ofp, "%.4e\t", EmisLines[ipHE_LIKE][nelem][upperIndex][ipLo].Aul );
						}
						fprintf( ofp, "\n" );
					}
				}
			}
		}
	}
}

/*prt_He_like_DeparCoef routine to print departure coefficients for he-like species */
void prt_He_like_DeparCoef(long int nelem )
{
	long in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;
	
	for( is = 0; is<=1; ++is)
	{
		char chSpin[2][9]={"singlets","triplets"};

		ipFirstCollapsed= helike.numLevels[nelem]-helike.nCollapsed[nelem];
		nResolved = iso.quant_desig[ipHE_LIKE][nelem][ipFirstCollapsed-1].n;
		ASSERT( nResolved == helike.n_HighestResolved[nelem] );
		ASSERT(nResolved > 0 );
		
		/* give element number and spin */
		fprintf(ioQQQ," He-like %s  %s departure\n",
			elementnames.chElementSym[nelem],
			chSpin[is]);

		/* header with the l states */
		fprintf(ioQQQ," n\\l=>         ");
		for ( i =0; i < nResolved ; ++i)
		{
			fprintf(ioQQQ,"%2ld       ",i);
		}
		fprintf(ioQQQ,"\n");
	
		/* loop over prin quant numbers, one per line, with l across */
		for( in = is+1; in <= nResolved; ++in)
		{
			fprintf(ioQQQ," %2ld           ",in);
			for( il = 0; il < in; ++il)
			{
				if ( (in==2) && (il==1) && (is==1) )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][nelem][ipHe2p3P0]));
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][nelem][ipHe2p3P1]));
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][nelem][ipHe2p3P2]));
				}
				else
				{
					ipLo = QuantumNumbers2Index[nelem][in][il][is];
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][nelem][ipLo]));
				}
			}
			fprintf(ioQQQ,"\n");
		}
	}		
	/* above loop was over spin, now do collapsed levels, no spin or ang momen */
	for( il = ipFirstCollapsed; il < helike.numLevels[nelem]; ++il)
	{
		in = iso.quant_desig[ipHE_LIKE][nelem][il].n;

		/* prin quan number of collapsed levels */
		fprintf(ioQQQ," %2ld           ",in);
		fprintf( ioQQQ,PrintEfmt("%9.2e", iso.DepartCoef[ipHE_LIKE][nelem][il]));
		fprintf(ioQQQ,"\n");
	}		
	
	return;
}

/*prt_He_like_Pops routine to print level populations for he-like species */
void prt_He_like_Pops(long int nelem )
{
	long in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;

	for( is = 0; is<=1; ++is)
	{
		char chSpin[2][9]={"singlets","triplets"};

		ipFirstCollapsed= helike.numLevels[nelem]-helike.nCollapsed[nelem];
		nResolved = iso.quant_desig[ipHE_LIKE][nelem][ipFirstCollapsed-1].n;
		ASSERT( nResolved == helike.n_HighestResolved[nelem] );
		ASSERT(nResolved > 0 );
	
		/* give element number and spin */
		fprintf(ioQQQ," He-like %s  %s populations\n",
			elementnames.chElementSym[nelem],
			chSpin[is]);

		/* header with the l states */
		fprintf(ioQQQ," n\\l=>         ");
		for ( i =0; i < nResolved ; ++i)
		{
			fprintf(ioQQQ,"%2ld       ",i);
		}
		fprintf(ioQQQ,"\n");

		/* loop over prin quant numbers, one per line, with l across */
		for( in = is+1; in <= nResolved; ++in)
		{
			fprintf(ioQQQ," %2ld           ",in);

			for( il = 0; il < in; ++il)
			{
				if ( (in==2) && (il==1) && (is==1) )
				{
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2p3P0]));
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2p3P1]));
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2p3P2]));
				}
				else
				{
					ipLo = QuantumNumbers2Index[nelem][in][il][is];
					fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][nelem][ipLo]));
				}
			}
			fprintf(ioQQQ,"\n");
		}
	}
	/* above loop was over spin, now do collapsed levels, no spin or ang momen */
	for( il = ipFirstCollapsed; il < helike.numLevels[nelem]; ++il)
	{
		in = iso.quant_desig[ipHE_LIKE][nelem][il].n;
		/* prin quan number of collapsed levels */
		fprintf(ioQQQ," %2ld           ",in);
		fprintf( ioQQQ,PrintEfmt("%9.2e", iso.Pop2Ion[ipHE_LIKE][nelem][il]));
		fprintf(ioQQQ,"\n");
	}

	return;
}
/*lint +e662 creation of  out of bound pointer */

#undef PARALLEL
/*lint +e662 creation of  out of bound pointer */
/*lint +e661 creation of  out of bound pointer */









