/*HydroCool compute net cooling due to hydrogenc atom species, ground state 
 * photoionization of hydrogenic species done in sumheat */
#include <limits.h>
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "tepowers.h"
#include "hstat.h"
#include "phycon.h"
#include "hdeltat.h"
#include "heating.h"
#include "hphoto.h"
#include "teinv.h"
#include "powi.h"
#include "hpheat.h"
#include "hfbnet.h"
#include "ionfracs.h"
#include "tcool.h"
#include "hcolrt.h"
#include "elecden.h"
#include "coladd.h"

void HydroCool(
		/* ipZ is charge -1, so 0 for H itself */
		long int ipZ)
{
	long int ipHi, 
	  ipLo, 
	  n;
	double RecCoolError, 
	  dbarest, 
	  tcoola, 
	  dhexdt, 
	  *dhone,  /* will become save arrays dimd nhlvl+1*/
	  dhrest, 
	  drestly, 
	  edenIonAbund, 
	  ChargeSquared,
	  hex, 
	  *hexx,  /* will become save arrays dimd nhlvl+1*/
	  *hlone, /* will become save arrays dimd nhlvl+1*/
	  ThinCoolingCaseB, 
	  ThinCoolingSum;

	double *SavePhotoHeat ,
		*SaveInducCool ,
		*SaveRadRecCool , 
		biggest =0.;

	/* these are labels for hydrogen line cooling */
	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 hydrogen collisional ionization cooling */
	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 hydrogen photoionization heating */
	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 hydrogen induced recombination heating */
	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( "<+>HydroCool()\n", debug_fp );
#	endif

	/* validate the incoming data */
	assert( ipZ >= 0 );
	assert( ipZ < LIMELM );

	/* for zero abundance we need to set some produced variables to zero */
	if( xIonFracs[ipZ+2][ipZ] <= 0. )
	{
		/* all global variables must be zeroed here or set below */
		hydro.chion[ipZ] = 0.;
		hydro.coola[ipZ] = 0.;
		hydro.restly[ipZ] = 0.;
		hydro.barest[ipZ] = 0.;
		hydro.hrest[ipZ] = 0.;
		hydro.HLineTotCool[ipZ] = 
		hydro.HRadRecCool[ipZ] = 0.;
		HFBNet.HeatExcited[ipZ] = 0.;
		HFBNet.HFreBndNet[ipZ] = 0.;
		HFBNet.crcind[ipZ] = 0.;
		/* hydrogenic species are special cases since their
		 * ground state heating is evaluated here rather than in bidiag */
		if( ipZ <= 1 )
		{
			/* heating by ground of hydrogen and Helium ion
			 * heavies were done in bidiag */
			HeatingCom.heating[ipZ][ipZ] = 0.;
		}

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

	/* now make some space */
	dhone = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );
	hexx = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );
	hlone = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );
	SavePhotoHeat = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );
	SaveInducCool = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );
	SaveRadRecCool = (double *)malloc( (unsigned)(nhlevel+1)*sizeof(double) );

	if( dhone==NULL || hexx==NULL || hlone==NULL || SavePhotoHeat==NULL || 
		SaveInducCool==NULL || SaveRadRecCool==NULL)
	{
		fprintf( ioQQQ," could not malloc space for scratch arrays\n");
		puts( "[Stop in HydroCool]" );
		exit(1);
	}

	/* this is energy factor that will appear in some places below */
	ChargeSquared = POW2(ipZ+1.);

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

	/* highest level is fastest collision by far, its coupling to
	 * continuum can overwhelm heating, so only count nearly up to the top */
	for( n=IP1S; n <= (nhlevel - 2); n++ )
	{
		/* total collisional cooling less three body heating for the atom */
		hexx[n] = 
			hydro.HLevNIonRyd[ipZ][n]*hydro.hcolnc[ipZ][n]*ElecDen.EdenHCorr*
		  (hydro.hn[ipZ][n] -hcolrt.hlte[ipZ][n]*phycon.eden);
		hex += hexx[n];
		/* the derivative of the cooling */
		dhexdt += hexx[n]*(HdeltaT.HCionTe[ipZ][n]*tcool.tsq1 - 
		  tcool.halfte);
	}

	/* convert to physical units, need to convert ryd to ergs, 
	 * and bring to density per vol not per ion */
	hex *= EN1RYD*xIonFracs[ipZ+2][ipZ];
	dhexdt *= EN1RYD*xIonFracs[ipZ+2][ipZ];

	/* save net hydrogen collisional ionization cooling less H-3 body heating
	 * for inclusion in printout */
	hydro.chion[ipZ] = hex;
	/* dump the coolant on the cooling stack */
	coladd(chHion[ipZ],0,MAX2(0.,hydro.chion[ipZ]));
	tcool.dCooldT += dhexdt;

	/* heating(1,4) is all three body heating, opposite of coll ion cooling */
	HeatingCom.heating[3][0] += MAX2(0.,-hydro.chion[ipZ]);

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

	/* this is the product of the ion abundance times the electron density,
	 * will multiply level pops which are stored relative to ion */
	edenIonAbund = phycon.eden*xIonFracs[ipZ+2][ipZ];

	/* lyman alpha cooling - total cooling is net rad loss */
	ipHi = IP2P;
	ipLo = IP1S;

	/* hydro.hcol is downward collision rate, not defined in hydrogenic arrays*/
	hydro.coola[ipZ] = 
		edenIonAbund* HydroLines[ipZ][ipHi][ipLo].EnergyErg *
		(hydro.hn[ipZ][IP1S]*hydro.hlbolt[ipZ][IP1S][IP2P]*
		(hydro.hcol[ipZ][IP1S][IP2S] + 3.*hydro.hcol[ipZ][IP1S][IP2P]) - 
	   hydro.hn[ipZ][IP2S]*hydro.hcol[ipZ][IP1S][IP2S] - 
		hydro.hn[ipZ][IP2P]*hydro.hcol[ipZ][IP1S][IP2P]);

	if( hydro.coola[ipZ] > 0. )
	{
		/* net coolant, exponential dominates deriv */
		tcoola = hydro.coola[ipZ]*(1.184e5*ChargeSquared*tcool.tsq1 - tcool.halfte);
	}
	else
	{
		/* net heating, te^-1/2 dominates */
		tcoola = -hydro.coola[ipZ]*tcool.halfte;
	}

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

	drestly = 0.;
	hydro.restly[ipZ] = 0.;
	for( n=3; n <= (nhlevel - 1); n++ )
	{
		hlone[n] = hydro.hcol[ipZ][IP1S][n]*
		  (hydro.hn[ipZ][IP1S]*hydro.hlbolt[ipZ][IP1S][n]*POW2((float)n) - 
		  hydro.hn[ipZ][n])* edenIonAbund*
		  HydroLines[ipZ][n][IP1S].EnergyErg;
		/* sum energy in higher lyman lines */
		hydro.restly[ipZ] += hlone[n];
		if( hlone[n] > 0. )
		{
			/* HdeltaTe[ipLo][hi], ip1s is 1s, ip2p is 2, etc */
			drestly += hlone[n]*
			  (HydroLines[ipZ][n][IP1S].EnergyK*tcool.tsq1 - tcool.halfte);
			  /*>>chng 99 jan 28, get rid of hdeltate, use stored temps instead */
			  /*(HdeltaT.HdeltaTe[IP1S][n]*ChargeSquared*tcool.tsq1 - tcool.halfte);*/
		}
		else
		{
			drestly += -hlone[n]*tcool.halfte;
		}
	}

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

	/* Balmer line collisional heating/cooling and derivative */
	hydro.barest[ipZ] = 0.;
	dbarest = 0.;
	for( ipHi=3; ipHi <= (nhlevel - 1); ipHi++ )
	{
		/* single line cooling */
		hlone[ipHi] = 
			hydro.hcol[ipZ][IP2P][ipHi]*((hydro.hn[ipZ][IP2S] + 
		  hydro.hn[ipZ][IP2P])*hydro.hlbolt[ipZ][IP2P][ipHi]*
		  POW2((float)ipHi) /4. - hydro.hn[ipZ][ipHi])*
		  edenIonAbund* HydroLines[ipZ][ipHi][IP2S].EnergyErg;

		hydro.barest[ipZ] += hlone[ipHi];
		/* count boltzmann factor from ground - not a typo */
		if( hlone[ipHi] > 0. )
		{
			dhone[ipHi] = 
			  /*>>chng 99 jan 28, get rid of hdeltate, use stored temps instead */
				/*hlone[ipHi]*(HdeltaT.HdeltaTe[IP1S][ipHi]*ChargeSquared*tcool.tsq1 - tcool.halfte);*/
				hlone[ipHi]*(HydroLines[ipZ][ipHi][IP1S].EnergyK*tcool.tsq1 - tcool.halfte);
		}
		else
		{
			/* heating agent */
			dhone[ipHi] = -hlone[ipHi]*ChargeSquared*tcool.halfte;
		}
		dbarest += dhone[ipHi];
	}

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

	hydro.hrest[ipZ] = 0.;
	dhrest = 0.;
	for( ipLo=3; ipLo <= (nhlevel - 2); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi <= (nhlevel - 1); ipHi++ )
		{
			hlone[ipHi] = 
			  hydro.hcol[ipZ][ipLo][ipHi]*(hydro.hn[ipZ][ipLo]*
			  hydro.hlbolt[ipZ][ipLo][ipHi]*hstat.HStatWght[ipHi]/
			  hstat.HStatWght[ipLo] - hydro.hn[ipZ][ipHi])*edenIonAbund*
			  HydroLines[ipZ][ipHi][ipLo].EnergyErg;

			hydro.hrest[ipZ] += hlone[ipHi];
			if( hlone[ipHi] > 0. )
			{
				/* for derivative taking the exponential from ground gave better
				 * representation of effects */
				dhrest += hlone[ipHi]*
		   	  /*>>chng 99 jan 28, get rid of hdeltate, use stored temps instead */
				  (HydroLines[ipZ][ipHi][IP1S].EnergyK*tcool.tsq1 - tcool.halfte);
				  /*(HdeltaT.HdeltaTe[IP1S][ipHi]*ChargeSquared*tcool.tsq1 - tcool.halfte);*/
			}
			else
			{
				dhrest += -hlone[ipHi]*tcool.halfte;
			}
		}
	}

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

	/* this is total hydrogen line cooling or heating */
	hydro.HLineTotCool[ipZ] = 
	  hydro.coola[ipZ] + hydro.restly[ipZ] + hydro.barest[ipZ] + 
	  hydro.hrest[ipZ];

	/* above can be either heating or coolant
	 * must add this to total heating or cooling, as appropriate */
	if( hydro.HLineTotCool[ipZ] > 0. )
	{
		/* hydrogen is a net coolant, HeatingCom.heating[23][0] set to 0 in coolr*/
		coladd(chHlin[ipZ],0,hydro.HLineTotCool[ipZ]);
		tcool.dCooldT += tcoola + drestly + dhrest + dbarest;
		hydro.dHLTot[ipZ] = 0.;
	}
	else
	{
		/* hydrogen is a net heat source, must add since all hydrogenic added here */
		HeatingCom.heating[23][0] -= hydro.HLineTotCool[ipZ];
		coladd(chHlin[ipZ],0,0.);
		hydro.dHLTot[ipZ] = -tcoola - drestly - dbarest - dhrest;
	}

	/***********************************************************************
	 *                                                                     *
	 * hydrogen recombination free-bound free bound cooling                *
	 *                                                                     *
	 ***********************************************************************/
	
	/* this is ground, not ncluded in case b sum.  0 in ground in vector, 1 in function */
	SaveRadRecCool[0] = HydroRecCool(1,ipZ)*hydro.hrec[ipZ][0][ipRecNetEsc];
	hydro.HRadRecCool[ipZ] = SaveRadRecCool[0];

	/* now do case b sum to compare with exact value below */
	ThinCoolingSum = 0.;

	/* in this loop HydroRecCool(2,ipZ) is both 2s and 2p */
	for( n=2; n <= (nhlevel - 1); n++ )
	{
		/* evaluate and hold optically thin spontaneous radiative cooling*/
		SaveRadRecCool[n] = HydroRecCool(n,ipZ);
		/* this is optically thin sum of cooling, will subtract from case b sum
		 * so that we can "top off" the model atom's cooling */
		ThinCoolingSum += SaveRadRecCool[n];
		/* this includes escaping fraction */
		hydro.HRadRecCool[ipZ] += SaveRadRecCool[n]*
			hydro.hrec[ipZ][n][ipRecNetEsc];
	}
	/* this is really for debugging - [2] from above is both 2s and 2p, [1] undefined
	 * will make them independent and correct, with [1] defined */
	SaveRadRecCool[IP2S] = SaveRadRecCool[IP2P] / 3.;
	SaveRadRecCool[IP2P] -= SaveRadRecCool[IP2S];

	/* following is expression for case b sum of 
	 * optically thin radiative recombination cooling */
	if( ipZ == 0 )
	{
		/*expansion for hydrogen itself */
		ThinCoolingCaseB = (-25.859117 + 
		  0.16229407*TePowers.telogn[0] + 
		  0.34912863*TePowers.telogn[1] - 
		  0.10615964*TePowers.telogn[2])/(1. + 
		  0.050866793*TePowers.telogn[0] - 
		  0.014118924*TePowers.telogn[1] + 
		  0.0044980897*TePowers.telogn[2] + 
		  6.0969594e-5*TePowers.telogn[3]);
	}
	else
	{
		/* same expansion but for hydrogen ions */
		ThinCoolingCaseB = (-25.859117 + 
		  0.16229407*(TePowers.telogn[0]-TePowers.sqlogz[ipZ]) + 
		  0.34912863*POW2(TePowers.telogn[0]-TePowers.sqlogz[ipZ]) - 
		  0.10615964*powi( (TePowers.telogn[0]-TePowers.sqlogz[ipZ]),3) )/(1. + 
		  0.050866793*(TePowers.telogn[0]-TePowers.sqlogz[ipZ]) - 
		  0.014118924*POW2(TePowers.telogn[0]-TePowers.sqlogz[ipZ]) + 
		  0.0044980897*powi( (TePowers.telogn[0]-TePowers.sqlogz[ipZ]),3) + 
		  6.0969594e-5*powi( (TePowers.telogn[0]-TePowers.sqlogz[ipZ]),4) );
	}
	
	/* now convert to linear cooling coefficient */
	ThinCoolingCaseB = POW3(1.+ipZ)*pow(10.,ThinCoolingCaseB)/(phycon.te/POW2(1.+ipZ) ) ;

	/* this is the error, expect positive since do not include infinite number of levels */
	RecCoolError = ThinCoolingCaseB - ThinCoolingSum;

	/* leave this message in for now.  very large atoms tend to overestimate cooling
	 * due to approximate form of cooling above level 15 */
	if( RecCoolError/ThinCoolingSum < -.04 )
	{
		fprintf( ioQQQ, " RecCoolError lt 0 in HydroCool, frac error=%e ipZ is %li\n", 
		  RecCoolError/ThinCoolingSum , ipZ );
	}
	hydro.HRadRecCool[ipZ] += RecCoolError*hydro.hrec[ipZ][nhlevel][ipRecNetEsc];

	/* this is now total free-bound cooling */
	hydro.HRadRecCool[ipZ] *= edenIonAbund;

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

	/* photoionization of excited levels
	 * HPhotHeat(n,ipZ) is photoionization heating due to level n,
	 * evaluated in routine hydrgn */
	HFBNet.HeatExcited[ipZ] = 0.;
	{
		enum {DEBUG = FALSE};
		int ipbig=-1000;
		for( n=IP2S; n <= (nhlevel - 1); n++ )
		{
			SavePhotoHeat[n] = xIonFracs[ipZ+2][ipZ]*hydro.hn[ipZ][n]*
			  HPHeat.HPhotHeat[ipZ][n];
			HFBNet.HeatExcited[ipZ] += SavePhotoHeat[n];
			if( SavePhotoHeat[n] > biggest )
			{
				biggest = SavePhotoHeat[n];
				ipbig = n;
			}
		}
		if( DEBUG  )
		{
			/* this was not done above */
			SavePhotoHeat[IP1S] = xIonFracs[ipZ+2][ipZ]*hydro.hn[ipZ][IP1S]*
			  HPHeat.HPhotHeat[ipZ][IP1S];
			fprintf(ioQQQ," ipZ=%ld ipbig=%d biggest=%g\n",ipZ,
				ipbig , biggest );
			for(n=IP1S; n<= (nhlevel - 1); ++n )
			{
				fprintf(ioQQQ," %g", SavePhotoHeat[n]/HFBNet.HeatExcited[ipZ]);
			}
			fprintf(ioQQQ," \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 */
	HFBNet.HFreBndNet[ipZ] = hydro.HRadRecCool[ipZ] - HFBNet.HeatExcited[ipZ];

	/* 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  */
	HeatingCom.heating[1][0] += MAX2(0.,-HFBNet.HFreBndNet[ipZ]);

	/* net free-bound cooling minus free-free heating,
	 * these have labels like "Hp 0","Hp 1" */
	coladd(chHp[ipZ], 0, MAX2(0.,HFBNet.HFreBndNet[ipZ]));

	/* if rec coef goes as T^-.8, but times T, so propto t^.2 */
	tcool.dCooldT += 0.2*HFBNet.HFreBndNet[ipZ]*teinvCom.teinv;


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

	HFBNet.crcind[ipZ] = 0.;
	for( n=IP1S; n <= (nhlevel - 1); n++ )
	{
		SaveInducCool[n] = HPhoto.cinduc[ipZ][n]*edenIonAbund;
		HFBNet.crcind[ipZ] += SaveInducCool[n];
		tcool.dCooldT += SaveInducCool[n]*
			(HdeltaT.HCionTe[ipZ][n]/phycon.te - 1.5)*teinvCom.teinv;
	}
	{
		enum {DEBUG = FALSE};
		if( DEBUG  )
		{
			fprintf(ioQQQ," ipZ=%ld \n",ipZ);
			for(n=IP1S; n<= (nhlevel - 1); ++n )
			{
				fprintf(ioQQQ," %g", SaveInducCool[n]/HFBNet.crcind[ipZ]);
			}
			fprintf(ioQQQ," \n");
		}
	}
	/* this cooling has labels like "Hn 0","Hn 1" - induced rec cooling */
	coladd(chHn[ipZ],0,HFBNet.crcind[ipZ]);

	{
		/* 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 */
		enum {DEBUG = FALSE};
		enum {LTEPOP=TRUE} ;
		/* special debug print for gas near ste */
		if( DEBUG  && (ipZ==1 || ipZ==0) )
		{
			/* ltepop is option to assume ste for rates */
			if( LTEPOP )
			{
				for(n=IP1S; n<nhlevel; ++n )
				{
					fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", ipZ,n,
						HPhoto.hgamnc[ipZ][n] *hcolrt.hlte[ipZ][n], 
						/* induced recom, intergral times hlte */
						hydro.hrec[ipZ][n][ipRecRad]+HPhoto.rinduc[ipZ][n]  ,
						/* induced rec */
						HPhoto.rinduc[ipZ][n]  ,
						/* spontaneous recombination */
						hydro.hrec[ipZ][n][ipRecRad] ,
						/* heating, followed by two processes that must balance it */
						HPHeat.HPhotHeat[ipZ][n]*hcolrt.hlte[ipZ][n], 
						HPhoto.cinduc[ipZ][n]+SaveRadRecCool[n] ,
						/* induced rec cooling, integral times hlte */
						HPhoto.cinduc[ipZ][n] ,
						/* free-bound cooling per unit vol */
						SaveRadRecCool[n] );
				}
			}
			else
			{
				for(n=IP1S; n<nhlevel; ++n )
				{
					fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", ipZ,n,
						HPhoto.hgamnc[ipZ][n] *xIonFracs[ipZ+2][ipZ]*hydro.hn[ipZ][n], 
						/* induced recom, intergral times hlte */
						hydro.hrec[ipZ][n][ipRecRad]*edenIonAbund+HPhoto.rinduc[ipZ][n] *edenIonAbund ,
						hydro.hrec[ipZ][n][ipRecRad]*edenIonAbund ,
						HPhoto.rinduc[ipZ][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( dhone );
	free( hexx );
	free( hlone );
	free( SavePhotoHeat );
	free( SaveInducCool );
	free( SaveRadRecCool );

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