/* 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 */
/*iso_cool compute net cooling due to hydrogenc atom species, ground state 
 * photoionization of hydrogenic species done in sumheat */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "phycon.h"
#include "dense.h"
#include "thermal.h"
#include "cooling.h"
#include "iso.h"
/*TODO	2	- if pc lint ever fixes this bug in their product, remove this -e */
/*lint -e661 Possible access of out-of-bounds pointer*/
/*lint -e662  (Warning -- Possible creation of out-of-bounds pointer  */
void iso_cool(
	   /* iso sequence, 0 for hydrogenic */
	   long int ipISO , 
		/* nelem is charge -1, so 0 for H itself */
		long int nelem)
{
	long int ipHi, 
	  ipbig,
	  ipLo, 
	  n;
	double RecCoolExtra, 
	  dhexdt, 
	  dCdT_all, 
	  edenHCorr_IonAbund, 
	  edenIonAbund, 
	  hex, 
	  HeatExcited,
	  *hexx,  /* will become save arrays dimd nhlvl+1*/
	  *hlone ,/* will become save arrays dimd nhlvl+1*/
	  ThinCoolingCaseB, 
	  ThinCoolingSum;

	double thin;

	double *SavePhotoHeat ,
		*SaveInducCool ,
		*SaveRadRecCool , 
		biggest =0.;
	double heat_max;
	long int nlo_heat_max  , nhi_heat_max;

	/* these are labels for iso line cooling, sum of ALL hydrogen lines,
	 * along the iso sequence */
	static char chHlin[LIMELM][5]=
	{"Hl 0","Hl 1","Hl 2","Hl 3","Hl 4","Hl 5","Hl 6","Hl 7","Hl 8","Hl 9",
	 "Hl10","Hl11","Hl12","Hl13","Hl14","Hl15","Hl16","Hl17","Hl18","Hl19",
	 "Hl20","Hl21","Hl22","Hl23","Hl24","Hl25","Hl26","Hl27","Hl28","Hl29"};

	/* these are labels for iso collisional ionization cooling,
	 * along the iso sequence */
	static char chHion[LIMELM][5]=
	{"Hi 0","Hi 1","Hi 2","Hi 3","Hi 4","Hi 5","Hi 6","Hi 7","Hi 8","Hi 9",
	 "Hi10","Hi11","Hi12","Hi13","Hi14","Hi15","Hi16","Hi17","Hi18","Hi19",
	 "Hi20","Hi21","Hi22","Hi23","Hi24","Hi25","Hi26","Hi27","Hi28","Hi29"};

	/* these are labels for iso photoionization heating,
	 * along the iso sequence */
	static char chHp[LIMELM][5]=
	{"Hp 0","Hp 1","Hp 2","Hp 3","Hp 4","Hp 5","Hp 6","Hp 7","Hp 8","Hp 9",
	 "Hp10","Hp11","Hp12","Hp13","Hp14","Hp15","Hp16","Hp17","Hp18","Hp19",
	 "Hp20","Hp21","Hp22","Hp23","Hp24","Hp25","Hp26","Hp27","Hp28","Hp29"};

	/* these are labels for iso induced recombination heating,
	 * along the iso sequence */
	static char chHn[LIMELM][5]=
	{"Hn 0","Hn 1","Hn 2","Hn 3","Hn 4","Hn 5","Hn 6","Hn 7","Hn 8","Hn 9",
	 "Hn10","Hn11","Hn12","Hn13","Hn14","Hn15","Hn16","Hn17","Hn18","Hn19",
	 "Hn20","Hn21","Hn22","Hn23","Hn24","Hn25","Hn26","Hn27","Hn28","Hn29"};

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

	/*TODO	2	this routine dumps all heating and cooling into only a very few lables,
	 * it would be best to break these out into individual labels that show element,
	 * iso sequence, and agent */

	/* validate the incoming data */
	ASSERT( nelem >= ipISO );
	ASSERT( ipISO <NISO );
	ASSERT( nelem < LIMELM );

	/* for zero abundance we need to set some produced variables to zero */
	if( dense.xIonDense[nelem][nelem+1-ipISO] <= 0. )
	{
		/* all global variables must be zeroed here or set below */
		iso.coll_ion[ipISO][nelem] = 0.;
		iso.cLya_cool[ipISO][nelem] = 0.;
		iso.cLyrest_cool[ipISO][nelem] = 0.;
		iso.cBal_cool[ipISO][nelem] = 0.;
		iso.cRest_cool[ipISO][nelem] = 0.;
		iso.xLineTotCool[ipISO][nelem] = 
		iso.RadRecCool[ipISO][nelem] = 0.;
		iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = 0.;
		iso.dLTot[ipISO][nelem] = 0.;
		iso.RecomInducCool_Rate[ipISO][nelem] = 0.;

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

	/* now make some space */
	hexx = (double *)MALLOC( (unsigned)(iso.numLevels[ipISO][nelem])*sizeof(double) );
	hlone = (double *)MALLOC( (unsigned)(iso.numLevels[ipISO][nelem])*sizeof(double) );
	SavePhotoHeat = (double *)MALLOC( (unsigned)(iso.numLevels[ipISO][nelem])*sizeof(double) );
	SaveInducCool = (double *)MALLOC( (unsigned)(iso.numLevels[ipISO][nelem])*sizeof(double) );
	SaveRadRecCool = (double *)MALLOC( (unsigned)(iso.numLevels[ipISO][nelem])*sizeof(double) );

	if( hexx==NULL || hlone==NULL || SavePhotoHeat==NULL || 
		SaveInducCool==NULL || SaveRadRecCool==NULL)
	{
		fprintf( ioQQQ," could not MALLOC space for scratch arrays\n");
		puts( "[Stop in iso_cool]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* compute the total cooling due to the hydrogen atom, and its derivative
	 * halfte is 1/(2T) */

	/***********************************************************************
	 *                                                                     *
	 * collisional ionization cooling, less three-body recombination  heat *
	 *                                                                     *
	 ***********************************************************************/

	/* hex will be net collisional ionization cooling, units erg/cm^3/s */
	hex = 0.;
	dhexdt = 0.;

	/* collisional ionization cooling, three body heating 
	 * highest level is fastest collision by far, its coupling to
	 * continuum can overwhelm heating, so only count nearly up to the top */
	/* >>chng 02 apr 25, from numLevels to numLevels-1 - departure coef for highest can be large */
	for( n=0; n < iso.numLevels[ipISO][nelem]-1; ++n )
	{
		/* total collisional ionization cooling less three body heating */
		hexx[n] = 
		  iso.xIsoLevNIonRyd[ipISO][nelem][n]*iso.ColIoniz[ipISO][nelem][n]*dense.EdenHCorr*
		  (iso.Pop2Ion[ipISO][nelem][n] -iso.PopLTE[ipISO][nelem][n]*dense.eden);
		hex += hexx[n];
		/* the derivative of the cooling */
		/* >>chng 01 sep 23, get rid of HCionTe */
		/* need extra factor of temp of 1 ryd since div by square of temp in ryd */
		dhexdt += hexx[n]*(iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd) - 
		  thermal.halfte);
	}

	/* convert to physical units, need to convert ryd to ergs, 
	 * and bring to density per vol not per ion */
	hex *= EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO];
	dhexdt *= EN1RYD*dense.xIonDense[nelem][nelem+1-ipISO];

	/* save net collisional ionization cooling less H-3 body heating
	 * for inclusion in printout */
	iso.coll_ion[ipISO][nelem] = hex;
	/* dump the coolant onto the cooling stack */
	CoolAdd(chHion[nelem],0,MAX2(0.,iso.coll_ion[ipISO][nelem]));
	thermal.dCooldT += dhexdt;

	/* heating[0][3] is all three body heating, opposite of coll ion cooling,
	 * would be unusual for this to be non-zero since would require excited
	 * state departure coefficients to be greater than unity */
	thermal.heating[0][3] += MAX2(0.,-iso.coll_ion[ipISO][nelem]);

	/***********************************************************************
	 *                                                                     *
	 * hydrogen recombination free-bound free bound cooling                *
	 *                                                                     *
	 ***********************************************************************/
	
	/* this is the product of the ion abundance times the electron density,
	 * will multiply level pops which are stored relative to ion
	 * EdenHCorr is used in level pops, so should be used here too */
	edenIonAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];

	/* this is the product of the ion abundance times the electron density with
	 * a small correction for the presence of neutrals, which should be used in
	 * reactions that involve collisions between states, but NOT radiative recombination
	 * but should be used in collisional ionization / three body recombination */
	edenHCorr_IonAbund = dense.EdenHCorr*dense.xIonDense[nelem][nelem+1-ipISO];

	/* >>chng 02 apr 23 */
	/* new logic, for any iso seq */
	/* now do case b sum to compare with exact value below */
	iso.RadRecCool[ipISO][nelem] = 0.;
	ThinCoolingSum = 0.;
	if( ipISO == ipH_LIKE )
	{
		/* do ground with special approximate fits to Ferland et al. '92 */
		n = 0;
		thin = HydroRecCool(
			/* n is the prin quantum number on the physical scale */
			1 , 
			/* nelem is the charge on the C scale, 0 is hydrogen */
			nelem);
	}
	else
	{
		/* this is the cooling before correction for optical depths */
		 thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
			/* arg is the scaled temperature, T * n^2 / Z^2, 
			 * n is prin quant number, Z is charge, 1 for H */
			HCoolRatio( 
			phycon.te * POW2( (double)iso.quant_desig[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))*
			/* convert results to energy per unit vol */
			phycon.te * BOLTZMANN;
	}
	/* the cooling, corrected for optical dephs */
	SaveRadRecCool[n] = iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] * thin;
	/* this is now total free-bound cooling */
	iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[n] * edenIonAbund ;

	/* radiative recombination cooling for all excited states */
	for( n=1; n < iso.numLevels[ipISO][nelem]; n++ )
	{
		/* this is the cooling before correction for optical depths */
		 thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
			/* arg is the scaled temperature, T * n^2 / Z^2, 
			 * n is prin quant number, Z is charge, 1 for H */
			HCoolRatio( 
			phycon.te * POW2( (double)iso.quant_desig[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))*
			/* convert results to energy per unit vol */
			phycon.te * BOLTZMANN;

		/* the cooling, corrected for optical dephs */
		SaveRadRecCool[n] = iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] * thin * edenIonAbund;
		/* this is now total free-bound cooling */
		iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[n];

		/* keep track of case b sum for topoff below */
		ThinCoolingSum += thin;
	}
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC  )
		{
			if( nelem==1 && ipISO==ipHE_LIKE)
			{
				for( n=0; n < (iso.numLevels[ipISO][nelem] - 1); n++ )
				{
					fprintf(ioQQQ,"\t%.2f",SaveRadRecCool[n]/ThinCoolingSum) ;
				}
				fprintf(ioQQQ,"\n") ;
			}
		}
	}

	/* following is expression for case b sum of 
	 * optically thin radiative recombination cooling - we will need to add any remainder to
	 * the sum from above - high precision is needed to get STE result to converge
	 * close to equilibrium - as in lte.in and induchhe.in */
	/* >>chng 02 apr 26, do not print for he-like - they are always negative at low t */
	if( ipISO == ipH_LIKE )
	{
		/* these expressions are only valid for hydrogenic sequence */
		if( nelem == 0 )
		{
			/*expansion for hydrogen itself */
			ThinCoolingCaseB = (-25.859117 + 
			0.16229407*phycon.telogn[0] + 
			0.34912863*phycon.telogn[1] - 
			0.10615964*phycon.telogn[2])/(1. + 
			0.050866793*phycon.telogn[0] - 
			0.014118924*phycon.telogn[1] + 
			0.0044980897*phycon.telogn[2] + 
			6.0969594e-5*phycon.telogn[3]);
		}
		else
		{
			/* same expansion but for hydrogen ions */
			ThinCoolingCaseB = (-25.859117 + 
			0.16229407*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) + 
			0.34912863*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) - 
			0.10615964*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) )/(1. + 
			0.050866793*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) - 
			0.014118924*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) + 
			0.0044980897*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) + 
			6.0969594e-5*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),4) );
		}
		
		/* now convert to linear cooling coefficient */
		ThinCoolingCaseB = POW3(1.+nelem-ipISO)*pow(10.,ThinCoolingCaseB)/(phycon.te/POW2(1.+nelem-ipISO) ) ;

		/* this is the error, expect positive since do not include infinite number of levels */
		RecCoolExtra = ThinCoolingCaseB - ThinCoolingSum;
	}
	else
	{
		ThinCoolingCaseB = 0. ;
		RecCoolExtra = 0.;
	}

	/* don't let the extra be negative - should be positive if H-like, negative for
	 * he-like only due to real difference in recombination coefficients */
	RecCoolExtra = MAX2(0., RecCoolExtra );

	/* add error onto total - this is significant for approach to ste */
	iso.RadRecCool[ipISO][nelem] += RecCoolExtra* edenIonAbund *iso.RadRecomb[ipISO][nelem][iso.numLevels[ipISO][nelem]-1][ipRecNetEsc];

	/***********************************************************************
	 *                                                                     
	 * heating  by photoionization of                                      
	 * excited states of all species                            
	 *                                                                     
	 ***********************************************************************/

	/* photoionization of excited levels
	 * HPhotHeat(n,nelem) is photoionization heating due to level n,
	 * evaluated in iso_photo */
	HeatExcited = 0.;
	ipbig = -1000;
	for( n=1; n < (iso.numLevels[ipISO][nelem] - 1); ++n )
	{
		SavePhotoHeat[n] = dense.xIonDense[nelem][nelem+1-ipISO]*iso.Pop2Ion[ipISO][nelem][n]*
			iso.PhotoHeat[ipISO][nelem][n];
		HeatExcited += SavePhotoHeat[n];
		if( SavePhotoHeat[n] > biggest )
		{
			biggest = SavePhotoHeat[n];
			ipbig = (int)n;
		}
	}
	{
		/* this block to debug heating due to photo of var n */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC  && ipISO==0 && nelem==0  && nzone > 700)
		{
			/* this was not done above */
			SavePhotoHeat[ipH1s] = dense.xIonDense[nelem][nelem+1-ipISO]*iso.Pop2Ion[ipISO][nelem][ipH1s]*
			  iso.PhotoHeat[ipISO][nelem][ipH1s];
			fprintf(ioQQQ,"ipISO = %li nelem=%li ipbig=%li biggest=%g HeatExcited=%.2e ctot=%.2e\n",
				ipISO , nelem,
				ipbig , 
				biggest,
				HeatExcited ,
				thermal.ctot);
			for(n=ipH1s; n< (iso.numLevels[ipISO][nelem] - 1); ++n )
			{
				fprintf(ioQQQ,"DEBUG phot heat%2li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 
					n,
					SavePhotoHeat[n]/HeatExcited,
					dense.xIonDense[nelem][nelem+1-ipISO],
					iso.Pop2Ion[ipISO][nelem][n],
					iso.PhotoHeat[ipISO][nelem][n],
					iso.gamnc[ipISO][nelem][n]);
			}
		}
	}

	/* HFreBndNet is net cooling due to 
	 * HRadRecCool is total radiative recombination cooling sum to all levels,
	 * net is this less n>=2 photoionization heating */
	iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = iso.RadRecCool[ipISO][nelem] - HeatExcited;
	/*fprintf(ioQQQ,"DEBUG Hn2\t%.3e\t%.3e\n",
		-iso.RadRecCool[ipISO][nelem]/SDIV(thermal.htot),
		HeatExcited/SDIV(thermal.htot));*/

	/* heating[1][0] is all excited state photoionization heating from ALL species,
	 * this is set to zero in coolr before loop where this is called.
	 * >>chng 96 jan 25, was -HeatExcited  */
	thermal.heating[0][1] += MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipISO][nelem]);

	/* net free-bound cooling minus free-free heating,
	 * these have labels like "Hp 0","Hp 1" */
	CoolAdd(chHp[nelem], 0, MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipISO][nelem]));

	/* if rec coef goes as T^-.8, but times T, so propto t^.2 */
	thermal.dCooldT += 0.2*iso.FreeBnd_net_Cool_Rate[ipISO][nelem]*phycon.teinv;

	/***********************************************************************
	 *                                                                     *
	 * hydrogen induced recombination cooling                              *
	 *                                                                     *
	 ***********************************************************************/

	iso.RecomInducCool_Rate[ipISO][nelem] = 0.;
	for( n=0; n < (iso.numLevels[ipISO][nelem] - 1); ++n )
	{
		/* >>chng 02 jan 22, removed cinduc, replace with RecomInducCool */
		SaveInducCool[n] = iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]*edenIonAbund;
		iso.RecomInducCool_Rate[ipISO][nelem] += SaveInducCool[n];
		thermal.dCooldT += SaveInducCool[n]*
			(iso.xIsoLevNIonRyd[ipISO][nelem][n]/phycon.te_ryd - 1.5)*phycon.teinv;
	}
	{
		/* print rec cool, induc rec cool, photo heat */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && ipISO==0 && nelem==5 /**/ )
		{
			fprintf(ioQQQ," ipISO=%li nelem=%li ctot = %.2e\n",
				ipISO,
				nelem,
				thermal.ctot);
			fprintf(ioQQQ,"sum\t%.2e\t%.2e\t%.2e\n", 
				HeatExcited,
				iso.RadRecCool[ipISO][nelem],
				iso.RecomInducCool_Rate[ipISO][nelem]);
			fprintf(ioQQQ,"sum\tp ht\tr cl\ti cl\n");

			for(n=0; n< (iso.numLevels[ipISO][nelem] - 1); ++n )
			{
				fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e  \n", 
					n,
					SavePhotoHeat[n],
					SaveRadRecCool[n],
					SaveInducCool[n] ,
					iso.RecomInducCool_Coef[ipISO][nelem][n],
					iso.PopLTE[ipISO][nelem][n],
					iso.RecomInducRate[ipISO][nelem][n]);
			}
			fprintf(ioQQQ," \n");
		}
	}
	/* this cooling has labels like "Hn 0","Hn 1" - induced rec cooling */
	CoolAdd(chHn[nelem],0,iso.RecomInducCool_Rate[ipISO][nelem]);

	/* lines are remainder of this routine */
	/* this will be the only loop in the final code
	 * find total collisional energy exchange */
	iso.xLineTotCool[ipISO][nelem] = 0.;
	dCdT_all = 0.;
	heat_max = 0.;
	nlo_heat_max = -1;
	nhi_heat_max = -1;

	for( ipLo=0; ipLo < iso.numLevels[ipISO][nelem]-2; ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipISO][nelem]-1; ipHi++ )
		{
			hlone[ipHi] = 
			  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
			  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
			  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
			  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
			  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

			iso.xLineTotCool[ipISO][nelem] += hlone[ipHi];
			{
				/* print rec cool, induc rec cool, photo heat */
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC && ipISO==ipHE_LIKE && nelem==ipHELIUM /**/ )
				{
					fprintf(ioQQQ,"iso_coolll %li %li %.2e %.2e %.2e %.2e\n",
						ipHi , ipLo, 
						iso.Pop2Ion[ipISO][nelem][ipHi] ,
						EmisLines[ipISO][nelem][ipHi][ipLo].ColUL , 
						EmisLines[ipISO][nelem][ipHi][ipLo].cs , 
						hlone[ipHi]/SDIV( thermal.htot ) );
				}
			}

			/* next get derivative */
			if( hlone[ipHi] > 0. )
			{
				/* usual case, this line was a net coolant */
				/* for derivative taking the exponential from ground gave better
				 * representation of effects */
				dCdT_all += hlone[ipHi]*
				  (EmisLines[ipISO][nelem][ipHi][ipH1s].EnergyK*thermal.tsq1 - thermal.halfte);
			}
			else
			{
				/* this line acted to heat the gas, remember which one it was,
				 * the following three vars never are used, but could be for
				 * debugging */
				if( hlone[ipHi] < heat_max )
				{
					heat_max = hlone[ipHi];
					nlo_heat_max = ipLo;
					nhi_heat_max = ipHi;
				} 
				dCdT_all -= hlone[ipHi]*thermal.halfte;
			}
		}
	}
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/\
		if( DEBUG_LOC )
		{
			if( ipISO==ipHE_LIKE && nelem==ipHELIUM )
				fprintf(ioQQQ,"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max );
		}
	}

	/***********************************************************************
	 *                                                                     *
	 * Lyman Lya all by itself                                             *
	 *                                                                     *
	 ***********************************************************************/

	/* lyman alpha cooling - total cooling is net rad loss */
	ipHi = ipH2p;
	ipLo = ipH1s;

	/* EmisLines[ipISO].ColUL is downward collision rate coefficient*/

	ipHi = ipH2p;
	ipLo = ipH1s;
	hlone[ipHi] = 
	  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
	  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
	  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
	  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
	  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

	ipHi = ipH2s;
	ipLo = ipH1s;
	hlone[ipHi] = 
	  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
	  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
	  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
	  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
	  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

	iso.cLya_cool[ipISO][nelem] = hlone[ipH2s]+hlone[ipH2p] ;

	/***********************************************************************
	 *                                                                     *
	 * Lyman line collisional heating/cooling for lines higher than Lya    *
	 *                                                                     *
	 ***********************************************************************/

	iso.cLyrest_cool[ipISO][nelem] = 0.;

	for( n=3; n < iso.numLevels[ipISO][nelem]; n++ )
	{
		hlone[n] = EmisLines[ipISO][nelem][n][ipH1s].ColUL*
		  (iso.Pop2Ion[ipISO][nelem][ipH1s]*iso.Boltzmann[ipISO][nelem][n][ipH1s]*POW2((double)n) - 
		  iso.Pop2Ion[ipISO][nelem][n])* edenHCorr_IonAbund*
		  EmisLines[ipISO][nelem][n][ipH1s].EnergyErg;
		/* sum energy in higher lyman lines */
		iso.cLyrest_cool[ipISO][nelem] += hlone[n];
	}

	/***********************************************************************
	 *                                                                     *
	 * Balmer lines                                                        *
	 *                                                                     *
	 ***********************************************************************/

	/* Balmer line collisional heating/cooling and derivative */
	iso.cBal_cool[ipISO][nelem] = 0.;
	/*for( ipHi=3; ipHi <= (iso.numLevels[ipISO][nelem] - 1); ipHi++ )*/
	for( ipHi=3; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ )
	{
		/* single line cooling */
		ipLo = ipH2s;
		hlone[ipHi] = 
		  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
		  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
		  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
		  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
		  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

		ipLo = ipH2p;
		hlone[ipHi] += 
		  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
		  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
		  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
		  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
		  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

		/* this is only used to add to line array */
		iso.cBal_cool[ipISO][nelem] += hlone[ipHi];

	}

	/***********************************************************************
	 *                                                                     *
	 * all hydrogen lines from Paschen on up, but not Balmer               *
	 *                                                                     *
	 ***********************************************************************/

	iso.cRest_cool[ipISO][nelem] = 0.;
	for( ipLo=3; ipLo < iso.numLevels[ipISO][nelem]-1; ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ )
		{
			hlone[ipHi] = 
			  EmisLines[ipISO][nelem][ipHi][ipLo].ColUL*(iso.Pop2Ion[ipISO][nelem][ipLo]*
			  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
			  iso.stat[ipISO][nelem][ipHi]/iso.stat[ipISO][nelem][ipLo] - 
			  iso.Pop2Ion[ipISO][nelem][ipHi])*edenHCorr_IonAbund*
			  EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg;

			iso.cRest_cool[ipISO][nelem] += hlone[ipHi];

		}
	}

	/***********************************************************************
	 *                                                                     *
	 * add total line heating or cooling into stacks, derivatives          *
	 *                                                                     *
	 ***********************************************************************/

	/* above can be either heating or coolant
	 * must add this to total heating or cooling, as appropriate */
	if( iso.xLineTotCool[ipISO][nelem] > 0. )
	{
		/* hydrogen is a net coolant */
		CoolAdd(chHlin[nelem],0,iso.xLineTotCool[ipISO][nelem]);
		thermal.dCooldT += dCdT_all;
		iso.dLTot[ipISO][nelem] = 0.;
	}
	else
	{
		/* hydrogen is a net heat source, thermal.heating[0][23]was set to 0 in coolr*/
		thermal.heating[0][23] -= iso.xLineTotCool[ipISO][nelem];
		CoolAdd(chHlin[nelem],0,0.);
		iso.dLTot[ipISO][nelem] = -dCdT_all;
	}
#	if 0
	if( 0 && nzone>=134 && ipISO==1 && nelem==1 )
	{
		fprintf(ioQQQ,"isocoolll\t%.4e\t%.2e\t%.2e\n", 
			phycon.te, iso.xLineTotCool[ipISO][nelem],
			dense.xIonDense[ipHELIUM][1] );
	}
#	endif

	/* debug print for understanding contributors to HLineTotCool */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/\
		if( DEBUG_LOC )
		{
			if( nelem == 0 )
			{
				fprintf(ioQQQ,"%.2e la %.2f restly %.2f barest %.2f hrest %.2f\n",
					iso.xLineTotCool[ipISO][nelem] ,
					iso.cLya_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
					iso.cLyrest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
					iso.cBal_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
					iso.cRest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] );

#				if 0
				if( iso.cLya_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] > 0.9 )
				fprintf(ioQQQ," 1s 2s 2p %.2e %.2e %.2e\n",
					iso.Pop2Ion[ipISO][nelem][ipH1s] , 
					iso.Pop2Ion[ipISO][nelem][ipH2s] ,
					iso.Pop2Ion[ipISO][nelem][ipH2p]);

#				endif
			}
		}
	}
	{
		/* this is an option to print out each rate affecting each level in strict ste,
		 * this is a check that the rates are indeed in detailed balance */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		enum {LTEPOP=TRUE} ;
		/*@+redef@*/
		/* special debug print for gas near ste */
		if( DEBUG_LOC  && (nelem==1 || nelem==0) )
		{
			/* ltepop is option to assume ste for rates */
			if( LTEPOP )
			{
				for(n=ipH1s; n<iso.numLevels[ipISO][nelem]-1; ++n )
				{
					fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
						iso.gamnc[ipISO][nelem][n] *iso.PopLTE[ipISO][nelem][n], 
						/* induced recom, intergral times hlte */
						/*iso.RadRecomb[ipISO][nelem][n][ipRecRad]+iso.rinduc[ipISO][nelem][n]  ,*/
						/* >>chng 02 jan 22, remove rinduc, replace with RecomInducRate */
						iso.RadRecomb[ipISO][nelem][n][ipRecRad]+ 
							iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]  ,
						/* induced rec */
						iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]  ,
						/* spontaneous recombination */
						iso.RadRecomb[ipISO][nelem][n][ipRecRad] ,
						/* heating, followed by two processes that must balance it */
						iso.PhotoHeat[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n], 
						iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]+SaveRadRecCool[n] ,
						/* induced rec cooling, integral times hlte */
						iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] ,
						/* free-bound cooling per unit vol */
						SaveRadRecCool[n] );
				}
			}
			else
			{
				for(n=ipH1s; n<iso.numLevels[ipISO][nelem]-1; ++n )
				{
					fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
						iso.gamnc[ipISO][nelem][n] *dense.xIonDense[nelem][nelem+1-ipISO]*iso.Pop2Ion[ipISO][nelem][n], 
						/* induced recom, intergral times hlte */
						iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund+
							iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund ,
						iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund ,
						iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund ,
						/* heating, followed by two processes that must balance it */
						SavePhotoHeat[n], 
						SaveInducCool[n]+SaveRadRecCool[n]*edenIonAbund ,
						/* induced rec cooling, integral times hlte */
						SaveInducCool[n] ,
						/* free-bound cooling per unit vol */
						SaveRadRecCool[n]*edenIonAbund );
				}
			}
		}
	}

	/* now free up the space */
	free( hexx );
	free( hlone );
	free( SavePhotoHeat );
	free( SaveInducCool );
	free( SaveRadRecCool );

#	ifdef DEBUG_FUN
	fputs( " <->iso_cool()\n", debug_fp );
#	endif
	return;
}
/*lint +e661 Possible access of out-of-bounds pointer*/
/*lint +e662 (Warning -- Possible creation of out-of-bounds pointer  */

