/* 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 */
/*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.*/
/*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 "elementnames.h"
#include "taulines.h"
#include "phycon.h"
#include "thermal.h"
#include "dense.h"
#include "iso.h"
#include "helike.h"
#include "helike_recom.h"
#include "helike_cs.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]; */

/*
static void printstuff(long flag);
*/
/*static long ipLev,globalZ;*/

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

/* This routine handles errors when that option is turned on (via the command
 * "atom he-like error generation" */
void putError(long int nelem,
			  long int ipHi,
			  long int ipLo,
			  long int whichData,
			  float errorToPut)
{
	if( helike.lgRandErrGen )
	{
		/* whichData is either IPRAD, IPCOLLIS, or IPENERGY */
		ASSERT( whichData <= 2 );
		ASSERT( nelem < LIMELM );
		ASSERT( ipHi <= iso.numLevels_max[ipHE_LIKE][nelem] );
		ASSERT( ipLo <= iso.numLevels_max[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,
		ipISO=ipHE_LIKE;
	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)  
				/* >>chng 04 dec 07, add this test, some models are very ionized may not include
				 * helike */
				/* >>chng 05 mar 21, chng from < nelem to <= nelem, for case of HeI rec, He0
				 * abundance does not matter, only He+ abundance */
				&& dense.IonLow[nelem] <= nelem )
			{
				/* truncate atom if physical conditions limit the maximum principal quantum number of a
				 * bound electron to a number less than the malloc'd size */
				iso_continuum_lower( ipHE_LIKE, nelem );
			
				/* generate gaussian errors */
				if( helike.lgRandErrGen && nzone==0 )
				{
					HeLikeError(nelem);
				}

				/* evaluate recombination rates */
				/* if the flag helike.lgHugeCaseB is true, cloudy exits in HeRecom(). */
				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_max[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.;
					}
				}
			}
			/* option to force ionization */
			if( dense.lgSetIoniz[nelem] )
			{
				dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem];
				dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem];
				if( dense.SetIoniz[nelem][nelem+1-ipISO] > SMALLFLOAT )
				{
					iso.Pop2Ion[ipISO][nelem][ipH1s] = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
					EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = dense.SetIoniz[nelem][nelem-ipISO] / dense.SetIoniz[nelem][nelem+1-ipISO];
				}
				else
				{
					iso.Pop2Ion[ipISO][nelem][ipH1s] = 0.;
					EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo = 0.;
				}
			}
		}
	}
	dense.IonLow[ipHELIUM] = lowsav;
	dense.IonHigh[ipHELIUM] = ihisav;

	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_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[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_max[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_max[ipHE_LIKE][ipHELIUM],
					iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[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_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[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 int 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 
/* 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.};
	double 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 = TeList[i];
		tfidle(FALSE);

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

	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
}


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

		ipFirstCollapsed= iso.numLevels_local[ipHE_LIKE][nelem]-iso.nCollapsed_local[ipHE_LIKE][nelem];
		nResolved = iso.quant_desig[ipHE_LIKE][nelem][ipFirstCollapsed-1].n;
		ASSERT( nResolved == iso.n_HighestResolved_local[ipHE_LIKE][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 < iso.numLevels_local[ipHE_LIKE][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 int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;

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

		ipFirstCollapsed= iso.numLevels_local[ipHE_LIKE][nelem]-iso.nCollapsed_local[ipHE_LIKE][nelem];
		nResolved = iso.quant_desig[ipHE_LIKE][nelem][ipFirstCollapsed-1].n;
		ASSERT( nResolved == iso.n_HighestResolved_local[ipHE_LIKE][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 < iso.numLevels_local[ipHE_LIKE][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 */



