/* 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 */
/*HeatSum evaluate heating and secondary ionization for current conditions */
/*HeatZero is called by ConvBase */
#include "cddefines.h"
#include "physconst.h"
#include "thermal.h"
#include "heavy.h"
#include "hextra.h"
#include "trace.h"
#include "secondaries.h"
#include "converge.h"
#include "called.h"
#include "coolheavy.h"
#include "iso.h"
#include "co.h"
#include "hmi.h"
#include "dense.h"
#include "ionbal.h"
#include "phycon.h"
#include "numderiv.h"
#include "atomfeii.h"
#include "atmdat.h"
#include "cooling.h"
#include "grainvar.h"
/* this is the faintest relative heating we will print */
#define	FAINT_HEAT	0.02
/*lint -e661 possible creation of out of bounds pointer */
/*lint -e662 possible creation of out of bounds pointer */
/*lint -e771 out of bounds pointer */

#define PRT_DERIV	FALSE
void HeatSum( void )
{
	/* use to dim some vectors */
#	ifdef NDIM
#		undef NDIM
#	endif
#	define NDIM 40

	/* secondary ionization and exciation due to compton scattering */
	float cosmic_ray_ionization_rate , 
		pair_production_ionization_rate , 
		fast_neutron_ionization_rate , scmpla;

	/* ratio of electron to nuclei density, could be >1 because of he, set in eden_sum 
	double ElecFrac;*/

	/* ionization and excitation rates from hydrogen itself */
	double SeconIoniz_iso[NISO] , 
		SeconExcit_iso[NISO] ;

	long int i, 
	  ion,
	  ipnt, 
	  ipsave[NDIM], 
	  j, 
	  jpsave[NDIM], 
	  limit ,
	  nelem ;
	double HeatingLo ,
		HeatingHi ,
		secmet ,
		smetla;
	long ipISO,
		ns;
	static long int nhit=0, 
	  nzSave=0;
	double photo_heat_2lev_atoms,
		photo_heat_ISO_atoms ,
		OtherHeat , 
		deriv, 
		oldfac, 
		save[NDIM];
	static double oldheat=0., 
	  oldtemp=0.;
	double secmetsave[LIMELM];

	float SaveOxygen1 , SaveCarbon1 ;

	/*total number of colliders (atoms, ions, molecules) in gas, no electrons */
	double AtomicCollidDensity;

	/* routine to sum up total heating, and print agents if needed
	 * it also computes the true deriv, dH/dT */
	float xe;
	double Cosmic_ray_heat_eff ,
		Cosmic_ray_sec2prim;
	float sec2prim_par_1;
	float sec2prim_par_2;
	
#	ifdef DEBUG_FUN
	fputs( "<+>HeatSum()\n", debug_fp );
#	endif

	/*******************************************************************
	 *
	 * reevaluate the secondary ionization and excitation rates 
	 *
	 *******************************************************************/
	/* >>chng 03 apr 29, move eval of AtomicCollidDensity to here from PresTotCurrent 
	 * since only used for secondary ionization */
	/* this is total neutral particle density for collisional ionization */
	AtomicCollidDensity = 0.;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		AtomicCollidDensity += dense.xIonDense[nelem][0];
	}
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf(ioQQQ," xIonDense AtomicCollidDensity tot\t%.3e",AtomicCollidDensity);
			for( nelem=0; nelem < LIMELM; nelem++ )
			{
				fprintf(ioQQQ,"\t%.2e",dense.xIonDense[nelem][0]);
			}
			fprintf(ioQQQ,"\n");
		}
	}

	/* now add all the heavy molecules */
	for( i=0; i < NUM_HEAVY_MOLEC; i++ )
	{
		AtomicCollidDensity += co.hevmol[i];
	}

	AtomicCollidDensity += (float)(hmi.Hmolec[ipMH2p] + hmi.Hmolec[ipMHm] + hmi.Hmolec[ipMH3p] + 
	/*2.*hmi.Hmolec[ipMH2g]);*/
	2.*hmi.H2_total);

	/* ElecFrac is electron fraction, used for secondary ionization efficiency
	ElecFrac = dense.eden/dense.xNucleiTotal; */
	/* >>chng 03 apr 30, ElecFrac had been above, relative to total number of nuclei,
	 * which includes ions.  Should have been total number of atoms, which is AtomicCollidDensity 
	ElecFrac = dense.eden/SDIV(AtomicCollidDensity);*/

	/* ElecFrac is electron fraction */
	/* XE = EDEN/(total neutrals, roughly H)
	 * analytic fits to Shull and Van Steenberg (1985; Ap.J. 298, 268).
	 * lgSecOFF turns off secondary ionizations, sets heating effic to 100% */

	/* at very low ionization - as per 
	 * >>>refer	sec	ioniz	Xu and McCray 1991, Ap.J. 375, 190.
	 * everything goes to asymptote that is not present in Shull and
	 * Van Steenberg - do this by not letting ElecFrac fall below 1e-4 */
	xe = (float)(MAX2(dense.eden/dense.gas_phase[ipHYDROGEN],1e-4));

	if( secondaries.lgSecOFF || xe > 0.95 )
	{
		secondaries.heatef = 1.;
		secondaries.efionz = 0.;
		secondaries.exctef = 0.;
		Cosmic_ray_heat_eff = 1.;
		Cosmic_ray_sec2prim = 0.;
	}
	/* >>chng 03 apr 29, only eval one time per zone since drove oscillations
	 * in He+ - He0 ionization front in very high Z models, like hizqso */
	else 
	{

		/* this is heating efficiency for high-energy photoejections.  It is the ratio
		 * of the final heat added to the energy of the photo electron.  Fully
		 * ionized, this is unity, and less than unity for neutral media.  
		 * It is a scale factor that multiplies the
		 * high energy heating rate */
		secondaries.heatef = (float)(0.9971*(1. - pow(1. - pow(xe,0.2663),1.3163)));

		/* secondary ionizations - this is the number of secondary ionizations
		 * produced for each ryd of photo electron energy.  It only multiplies the high-energy
		 * heating rate.  It is zero for a fully ionized gas */
		secondaries.sec2prim = (float)(0.3908*pow(1. - pow(xe,0.4092),1.7592));

		/* by dividing by the energy of one ryd, we can use this factor to multiply
		 * the heating rate in ergs */
		secondaries.efionz = (float)(secondaries.sec2prim/EN1RYD);

		/* This is new value, derived to approximate the curve given in Shull and
		 * van Steeberg for cosmic rays at 35 eV */		
		sec2prim_par_1 = (float)(-(1.251 + 195.932 * xe));
		sec2prim_par_2 = (float)((1 + 46.814 * xe - 44.969 * xe * xe));

		Cosmic_ray_sec2prim = (exp(sec2prim_par_1/SDIV( sec2prim_par_2)));

		/*  secondary excitation of L-alpha, actually all Lyman lines 
		 *
		 *  Note--This formula is derived for primary energies greater 
		 *  than 100 eV and is only an approximation.  This will
		 *  overpredict the secondary ionization of L-alpha.  We cannot make
		 *  fitting functions for this equation at low energies like we did
		 *  for the heating efficiency and for the number of secondaries
		 *  because the Shull and van Steeberg paper did not publish a similiar
		 *  curve for L-alpha */

		secondaries.exctef = (float)(0.4766*pow(1. - pow(xe,0.2735),1.5221));
		secondaries.exctef /= (float)EN1RYD;
		if( trace.lgTrace && trace.lgSecIon )
		{
			fprintf( ioQQQ, 
				"   HeatSum eval sec ion effic, xe = %.3e heatef %.3e sec2prim %.3e efioniz %.3e\n", 
				xe,
				secondaries.heatef ,
				secondaries.sec2prim ,
				secondaries.efionz );
		}

		/* cosmic ray heating, counts as non-ionizing heat since already result of cascade,
		* this was 35 eV per secondary ionization */
		/* We want to use the heating efficiency that is applicable to a 35 eV
		* primary electron, the current efficiency used is for 100eV cosmic ray
		* Here we use the correct value as given in Wolfire et al. 1995 */

		/* In general the equation for the cosmic ray heating rate is:*
		*															  *
		*															  *
		* CR_Heat_Rate = (density)*(cosmic ray ionization rate)*     *
		*				  (energy of primary electron)*               *
		*				  (Heating efficiency at that energy)		  *
		*															  *
		* The product of the last two terms gives the amount of heat *
		* added by the primary electron, and is dependent upon the   *
		* electron fraction.  We are using the same average primary  *
		* electron energy as Wolfire et al. (1995).  We do not,	  *
		* however, use their formula for the heating efficiency.     *
		* This is because the coefficients (f1, f2, and f3) were     *
		* originally intended for primary electron energies >100eV.  *
		* Instead we derived a heating efficiency based on the curves*
		* given in Shull and van Steenberg (1985).  We interpolated  *
		* for an energy of 35 eV the heating effieciency for electron*
		* fractions between 1e-4 and 1								  *
		**************************************************************/

		
		/*printf("Here is xe:\t%.3e\n", xe);*/
		Cosmic_ray_heat_eff =(float) (- 8.189 - 18.394 * xe - 6.608 * xe * xe * log(xe)
			+ 8.322 * exp( xe )  + 4.961 * sqrt(xe));
	}

	/*******************************************************************
	 *
	 * get total heating from all species
	 *
	 *******************************************************************/

	/* get total heating */
	photo_heat_2lev_atoms = 0.;
	photo_heat_ISO_atoms = 0.;

	/* add in effects of high-energy opacity of CO using C and O atomic opacities
	 * k-shell and valence photo of C and O in CO is not explicitly counted elsewhere
	 * this trick roughly accounts for it*/
	SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
	SaveCarbon1 = dense.xIonDense[ipCARBON][0];

	/* atomic C and O will include CO during the heating sum loop */
	dense.xIonDense[ipOXYGEN][0] += co.hevmol[ipCO];
	dense.xIonDense[ipCARBON][0] += co.hevmol[ipCO];

	/* this will hold cooling due to metal collisional ionization */
	CoolHeavy.colmet = 0.;
	/* metals secondary ionization, Lya excitation */
	secmet = 0.;
	smetla = 0.;
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		SeconIoniz_iso[ipISO] = 0.;
		SeconExcit_iso[ipISO] = 0.;
	}

	/* this loop includes hydrogenic, which is special case due to presence
	 * of substantial excited state populations */
	for( nelem=0; nelem<LIMELM; ++nelem)
	{
		secmetsave[nelem] = 0.;
		if( dense.lgElmtOn[nelem] )
		{
			/* sum heat over all stages of ionization that exist */
			/* first do the iso sequences,
			 * h-like and he-like are special because full atom always done, 
			 * can be substantial,  pops in excited states, with little in ground 
			 * (true near lte), these are done in following loop */

			limit = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );

			/* loop over all elements/ions done with as two-level systems */
			for( ion=dense.IonLow[nelem]; ion < limit; ion++ )
			{
				/* this will be heating for a single stage of ionization */
				HeatingLo = 0.;
				HeatingHi = 0.;
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					/* heating by various sub-shells */
					HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
					HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
				}

				/* total heating, all shells, for this stage of ionization */
				thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion]* 
					(HeatingLo + HeatingHi*secondaries.heatef);

				/* add to total heating */
				photo_heat_2lev_atoms += thermal.heating[nelem][ion];
				/*if( nzone>290 && thermal.heating[nelem][ion]/thermal.htot>0.01 )
					fprintf(ioQQQ,"buggg\t%li %li %.3f\n", nelem,ion,thermal.heating[nelem][ion]/thermal.htot);*/

				/* Cooling due to collisional ionization of heavy elements by thermal electrons
				 * CollidRate[nelem][ion][1] is cooling, erg/s/atom, eval in ion_collis */
				CoolHeavy.colmet += dense.xIonDense[nelem][ion]*ionbal.CollIonRate_Ground[nelem][ion][1];

				/* secondary ionization rate,  */
				secmetsave[nelem] += HeatingHi*secondaries.efionz* dense.xIonDense[nelem][ion];

				/* LA excit rate, =0 if ionized */
				smetla += HeatingHi*secondaries.exctef* dense.xIonDense[nelem][ion];
			}
			secmet += secmetsave[nelem];

			/* this branch loop over all ions done with full iso sequence */
			limit = MAX2( limit, dense.IonLow[nelem] );
			for( ion=MAX2(0,limit); ion < dense.IonHigh[nelem] ; ion++ )
			{
				/* this is the iso sequence */
				ipISO = nelem-ion;
				/* this will be heating for a single stage of ionization */
				HeatingLo = 0.;
				HeatingHi = 0.;
				/* the outer shell contains the compton recoil part */
				for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
				{
					/* heating by low energy, and then high energy, light */
					HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
					HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
				}

				/* net heating */
				thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion+1]*
				  iso.Pop2Ion[ipISO][nelem][0]*(HeatingLo + HeatingHi*secondaries.heatef);

				/* add to total heating */
				photo_heat_ISO_atoms += thermal.heating[nelem][ion];

				/* secondary ionization rate,  */
				SeconIoniz_iso[ipISO] += HeatingHi*secondaries.efionz* 
					iso.Pop2Ion[ipISO][nelem][0]*dense.xIonDense[nelem][ion+1]/AtomicCollidDensity;

				/* LA excit rate, =0 if ionized */
				SeconExcit_iso[ipISO] += HeatingHi*secondaries.exctef* 
					iso.Pop2Ion[ipISO][nelem][0]*dense.xIonDense[nelem][ion+1]/AtomicCollidDensity;

				ASSERT( SeconIoniz_iso[ipISO]>=0. && 
					SeconExcit_iso[ipISO]>=0.);
			}

			/* make sure stages of ionization with no abundances,
			 * also has no heating */
			for( ion=0; ion<dense.IonLow[nelem]; ion++ )
			{
				ASSERT( thermal.heating[nelem][ion] == 0. );
			}
			for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ion++ )
			{
				ASSERT( thermal.heating[nelem][ion] == 0. );
			}
		}
	}
	if( trace.lgTrace && trace.lgSecIon )
	{
		double savemax=0.;
		long int ipsavemax=-1;
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			if( secmetsave[nelem] > savemax )
			{
				savemax = secmetsave[nelem];
				ipsavemax = nelem;
			}
		}
		fprintf( ioQQQ, 
			"   HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n", 
			ipsavemax,
			savemax/SDIV(secmet),
			secmet);
	}

	/* now reset the abundances */
	dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
	dense.xIonDense[ipCARBON][0] = SaveCarbon1;

	/* convert secmet to proper final units */
	/*fprintf(ioQQQ,"bugggg\t%li\t%.3e\t%.3e\t%.3e\n", nzone , 
		secmet , AtomicCollidDensity , secmet / AtomicCollidDensity );*/
	secmet /= AtomicCollidDensity;
	smetla /= AtomicCollidDensity;

	/* >>chng 01 dec 20, do full some over all secondaries */
	/* bound compton recoil heating */
	/* >>chng 02 mar 28, save heating in this var rather than heating[0][18] 
	 * since now saved in photo heat 
	 * this is only used for a printout and in lines, not as heat since already counted*/
	ionbal.CompRecoilHeatLocal = 0.;
	for( nelem=0; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<nelem+1; ++ion )
		{
			ionbal.CompRecoilHeatLocal += 
				ionbal.CompRecoilHeatRate[nelem][ion]*secondaries.heatef*dense.xIonDense[nelem][ion];
		}
	}

	/* find heating due to charge transfer */
	thermal.heating[0][24] = ChargTranSumHeat();

	/* heating due to pair production */
	thermal.heating[0][21] = 
		ionbal.PairProducPhotoRate[2]*secondaries.heatef*(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM]);
	/* last term above is number of nuclei in helium */

	/* this is heating due to fast neutrons, assumed to secondary ionize */
	thermal.heating[0][20] = 
		ionbal.xNeutronHeatRate*secondaries.heatef*dense.gas_phase[ipHYDROGEN];

	/* turbulent heating, assumed to be a non-ionizing heat agent, added here */
	thermal.heating[0][20] += ionbal.ExtraHeatRate;

	/*printf("efficiency is:\t%.3e\nElectron fraction is:\t%.5e\n",Cosmic_ray_heat_eff, xe);*/

	thermal.heating[1][6] += 
		ionbal.CosRayHeatRate*
		((dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]) + 
		  dense.xIonDense[ipHELIUM][0])* 
		Cosmic_ray_heat_eff;
	/* >>chng 04 mar 24, do not multiply by further sec2prim since heat is 35 eV per primary */
	/*	secondaries.sec2prim;*/

	/* now sum up all heating agents not included in sum over normal bound levels above */
	OtherHeat = 0.;
	for( nelem=0; nelem<LIMELM; ++nelem)
	{
		/* we now have ionization solution for this element,sum heat over
		 * all stages of ionization that exist */
		/* >>>chng 99 may 08, following loop had started at nelem+3, and so missed [1][0],
		 * which is excited states of hydrogenic species.  increase this limit */
		for( i=nelem+1; i < LIMELM; i++ )
		{
			OtherHeat += thermal.heating[nelem][i];
		}
	}

	thermal.htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms;

	/* following checks whether total heating is strange, if we are not in search phase */
	if( called.lgTalk && !conv.lgSearch )
	{
		/* print this warning if not constant temperature and neg heat */
		if( thermal.htot < 0. && !thermal.lgTSetOn)
		{
			fprintf( ioQQQ, 
				" Total heating is <0; is this model collisionally ionized? zone is %li\n",
				nzone );
		}
		else if( thermal.htot == 0. )
		{
			fprintf( ioQQQ, 
				" Total heating is 0; is the density small? zone is %li\n",
				nzone);
		}
	}

	/* add on line heating to this array, heatl was evaluated in sumcool
	 * not zero, because of roundoff error */
	if( thermal.heatl/thermal.ctot < -1e-15 )
	{
		fprintf( ioQQQ, " HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n", 
		  thermal.heatl, thermal.ctot );

		fprintf( ioQQQ, " this is zone%4ld\n", nzone );
		ShowMe();
		puts( "[Stop in sumheat]" );
		cdEXIT(EXIT_FAILURE);
	}

	/*******************************************************************
	 *
	 * secondary ionization and excitation rates 
	 *
	 *******************************************************************/

	/* the terms cosmic_ray_ionization_rate & scmpla contain energies added in highen.  
	 * The only significant one is usually bound compton heating except when 
	 * cosmic rays are present */

	/* add on heating due to pair production, ionization then excitation */
	pair_production_ionization_rate = (float)(
		ionbal.PairProducPhotoRate[2]*secondaries.efionz*
		(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM])/AtomicCollidDensity);

	scmpla = (float)(
		ionbal.PairProducPhotoRate[2]*secondaries.exctef*1.3333F*
		(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM])/AtomicCollidDensity);

	/* this is heating due to fast neutrons, assumed to secondary ionize */
	fast_neutron_ionization_rate = (float)(
		ionbal.xNeutronHeatRate*secondaries.efionz*dense.gas_phase[ipHYDROGEN]/AtomicCollidDensity);

	scmpla += (float)(
		ionbal.xNeutronHeatRate*secondaries.exctef*1.3333F*dense.gas_phase[ipHYDROGEN]/
		AtomicCollidDensity);

	/* cosmic ray ionization */
	/* >>>chng 99 apr 29, term in PhotoRate was not present */
	cosmic_ray_ionization_rate = (float)(
		/* this term is cosmic ray ionization, set in highen, did not multiply by
		 * collider density in highen, so do not divide by it here */
		/* >>chng 99 jun 29, added efionz to cosmic ray rate, also multiply
		 * by number of secondaries per primary*/
		ionbal.CosRayIonRate*Cosmic_ray_sec2prim +
		/* this is the cosmic ray heating rate */
		ionbal.CosRayHeatRate*secondaries.efionz);

	/* cosmic ray Lya excitation rate */
	scmpla += (float)(
		ionbal.CosRayHeatRate*secondaries.exctef*1.3333f*
		/* multiply by atomic H and He densities */
		(dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHELIUM][0])/AtomicCollidDensity);


	/* find total suprathermal collisional ionization rate
	 * CSUPRA is H0 col ionize rate from non-thermal electrons (inv sec)
	 * X12tot is LA excit rate, inv sec
	 * SECCMP evaluated in HIGHEN, =ioniz rate for cosmic rays, sec bound */

	/* option to force secondary ionization rate, normally false */
	if( secondaries.lgCSetOn )
	{
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				secondaries.csupra[nelem][ion] = secondaries.SetCsupra*secondaries.csupra_effic[nelem][ion];
			}
		}
		secondaries.x12tot = secondaries.SetCsupra;
		/* this implements set leiden hac cosmic rays off command */
		if (hextra.lgLeidenCsupraHack)
		{
			secondaries.csupra[ipOXYGEN][0] = 0;
			secondaries.csupra[ipCARBON][0] = 0;
		}
	}
	else
	{
		double csupra;
		double facold , facnew;
		/* AtomicCollidDensity is total neutral particle density */
		/* >>chng 02 may 30, dense neutral model had eden - csupra oscil -
		 * try taking average of old and new */
		/*if( secondaries.csupra[ipHYDROGEN][0] / SDIV( hydro.DestRate[ipHYDROGEN][ipH1s] ) > 0.1 )*/
		if( secondaries.csupra[ipHYDROGEN][0] / SDIV( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] ) > 0.1 )
		{
			/* supra are dominant H destruction - make small changes */
			facold = 0.9;
		}
		else
		{
			/* secondaries are not important - only use new */
			facold = 0.;
		}
		facnew = 1. - facold;
		csupra = (secondaries.csupra[ipHYDROGEN][0]* facold + facnew*
			(cosmic_ray_ionization_rate + 
			pair_production_ionization_rate +
			fast_neutron_ionization_rate +
			SeconIoniz_iso[ipH_LIKE] + 
			SeconIoniz_iso[ipHE_LIKE] + 
			secmet ));

		/*TODO	2	find correct high-energy limit for these */
		/* now fill in ionization rates for all elements and ion stages */
		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				secondaries.csupra[nelem][ion] = (float)csupra*secondaries.csupra_effic[nelem][ion];
			}
		}

		/* now do same for Ly-alpha */
		secondaries.x12tot = (float)( secondaries.x12tot*facold + facnew*
			(scmpla + 
			SeconExcit_iso[ipH_LIKE] + 
			SeconExcit_iso[ipHE_LIKE] + 
			smetla));
	}

	if( trace.lgTrace && secondaries.csupra[ipHYDROGEN][0] > 0. )
	{
		fprintf( ioQQQ, 
			"   HeatSum retrn CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n", 
		  secondaries.csupra[ipHYDROGEN][0], 
		  (cosmic_ray_ionization_rate+pair_production_ionization_rate+fast_neutron_ionization_rate)/secondaries.csupra[ipHYDROGEN][0], 
		  SeconIoniz_iso[ipH_LIKE] / secondaries.csupra[ipHYDROGEN][0], 
		  SeconIoniz_iso[ipHE_LIKE]/secondaries.csupra[ipHYDROGEN][0], 
		  secmet/secondaries.csupra[ipHYDROGEN][0] ,
		  xe );
	}

	/*******************************************************************
	 *
	 * get derivative of total heating
	 *
	 *******************************************************************/

	/* now get derivative of heating due to photoionization, 
	 * >>chng 96 jan 14
	 *>>>>>NB heating decreasing with increasing temp is negative dH/dT */
	thermal.dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms)/phycon.te;
	/* >>chng 04 feb 28, add this correction factor - when gas totally neutral heating
	 * does not depend on temperature - when ionized depends on recom coef - this
	 * smoothly joins two limits */
	thermal.dHeatdT *= dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  0 %.3e  %.3e  %.3e \n", 
		thermal.dHeatdT,
		photo_heat_2lev_atoms,
		photo_heat_ISO_atoms);

	/* oldfac was factor used in old implimentation
	 * any term using it should be rethought */
	oldfac = -0.7/phycon.te;

	/* carbon monoxide */
	thermal.dHeatdT += thermal.heating[0][9]*oldfac;

	/* ly alpha collisoinal heating */
	thermal.dHeatdT += thermal.heating[0][10]*oldfac;

	/* line heating */
	thermal.dHeatdT += thermal.heating[0][22]*oldfac;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  1 %.3e \n", thermal.dHeatdT);

	/* free free heating, propto T^-0.5
	 * >>chng 96 oct 30, from heating(1,12) to CoolHeavy.brems_heat_total,
	 * do cooling separately assume CoolHeavy.brems_heat_total goes as T^-3/2
	 * dHTotDT = dHTotDT + heating(1,12) * (-0.5/te) */
	thermal.dHeatdT += CoolHeavy.brems_heat_total*(-1.5/phycon.te);

#if 0
	/* grain heating by photo
	 * >>chng 97 jul 02, change deriv in following, better fit to what came out */
	thermal.dHeatdT += thermal.heating[0][13]*0.1/phycon.te;

	/* grain heating by collsional */
	thermal.dHeatdT += thermal.heating[0][14]*oldfac;
#else
	/* >>chng 04 aug 07, use better estimate for heating derivative; needed in PDRs, PvH */
	/* this includes PE, thermionic, and collisional heating of the gas by the grains */
	thermal.dHeatdT += gv.dHeatdT;
#endif

	/* helium triplets heating */
	thermal.dHeatdT += thermal.heating[1][2]*oldfac;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  2 %.3e \n", thermal.dHeatdT);

	/* hydrogen molecule collisional deexcitation heating */
	/* >>chng 04 jan 25, add max to prevent cooling from entering here */
	/*thermal.dHeatdT += MAX2(0.,thermal.heating[0][8])*oldfac;*/
	if( hmi.HeatH2Dexc_used > 0. )
		thermal.dHeatdT += hmi.deriv_HeatH2Dexc_used;

	/* >>chng 96 nov 15, had been oldfac below, wrong sign
	 * H- H minus heating, goes up with temp since rad assoc does */
	thermal.dHeatdT += thermal.heating[0][15]*0.7/phycon.te;

	/* H2+ heating */
	thermal.dHeatdT += thermal.heating[0][16]*oldfac;

	/* Balmer continuum and all other excited states
	 * - T goes up so does pop and heating
	 * but this all screwed up by change in eden */
	thermal.dHeatdT += thermal.heating[0][1]*oldfac;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  3 %.3e \n", thermal.dHeatdT);

	/* all three body heating, opposite of coll ion cooling */
	thermal.dHeatdT += thermal.heating[0][3]*oldfac;

	/* bound electron recoil heating
	thermal.dHeatdT += ionbal.CompRecoilHeatLocal*oldfac; */

	/* compton heating */
	thermal.dHeatdT += thermal.heating[0][19]*oldfac;

	/* extra heating source, usually zero */
	thermal.dHeatdT += thermal.heating[0][20]*oldfac;

	/* pair production */
	thermal.dHeatdT += thermal.heating[0][21]*oldfac;
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  4 %.3e \n", thermal.dHeatdT);

	/* deriv of heating due to collisions of H lines, heating(1,24) */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem<LIMELM; ++nelem)
		{
			thermal.dHeatdT += iso.dLTot[ipISO][nelem];
		}
	}

	/* heating due to large FeII atom, zero unless atom feii is entered,
	 * FeII.feheat was entered into thermal.heating[25][27] */
	if( FeII.feheat > 0. )
	{
		thermal.dHeatdT += FeII.dfcool;
	}
	if( PRT_DERIV )
		fprintf(ioQQQ,"DEBUG dhdt  5 %.3e \n", thermal.dHeatdT);

	/* possibility of getting emperical heating derivative
	 * normally false, set true with 'set numerical derivatives' command */
	if( NumDeriv.lgNumDeriv )
	{
		/*lint -e777 test floats for equality */
		if( ((nzone > 2 && nzone == nzSave) && 
			oldtemp != phycon.te) && nhit > 4 )
		/*lint +e777 test floats for equality */
		{
			/* hnit is number of trys on this zone - use to stop numerical problems
			 * do not evaluate numerical deriv until well into solution */
			deriv = (oldheat - thermal.htot)/(oldtemp - phycon.te);
			thermal.dHeatdT = deriv;
		}
		else
		{
			deriv = thermal.dHeatdT;
		}

		/* this happens when new zone starts */
		if( nzone != nzSave )
		{
			nhit = 0;
		}
		nzSave = nzone;
		nhit += 1;
		oldheat = thermal.htot;
		oldtemp = phycon.te;
	}

	if( trace.lgTrace && trace.lgHeatBug )
	{
		ipnt = 0;
		/* this loops through the 2D array that contains all agents counted in htot */
		for( i=0; i < LIMELM; i++ )
		{
			for( j=0; j < LIMELM; j++ )
			{
				if( thermal.heating[i][j]/thermal.htot > FAINT_HEAT )
				{
					ipsave[ipnt] = i;
					jpsave[ipnt] = j;
					save[ipnt] = thermal.heating[i][j]/thermal.htot;
					/* increment ipnt, but do not let it get too big */
					ipnt = MIN2((long)NDIM-1,ipnt+1);
				}
			}
		}

		/* now check for possible line heating not counted in 1,23
		 * there should not be any significant heating source heree
		 * since they would not be counted in deriv correctly */
		for( i=0; i < thermal.ncltot; i++ )
		{
			if( thermal.heatnt[i]/thermal.htot > FAINT_HEAT )
			{
				ipsave[ipnt] = -1;
				jpsave[ipnt] = i;
				save[ipnt] = thermal.heatnt[i]/thermal.htot;
				ipnt = MIN2((long)NDIM-1,ipnt+1);
			}
		}

		fprintf( ioQQQ, 
		  "    HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n", 
		  thermal.htot, 
		  phycon.te, 
		  thermal.dHeatdT ,	
		  OtherHeat , 
		  photo_heat_2lev_atoms, 
		  photo_heat_ISO_atoms);

		fprintf( ioQQQ, "  " );
		for( i=0; i < ipnt; i++ )
		{
			/*lint -e644 these three are initialized above */
			fprintf( ioQQQ, "   [%ld][%ld]%6.3f",
				ipsave[i], 
				jpsave[i],
				save[i] );
			/*lint +e644 these three are initialized above */
			/* throw a cr every n numbers */
			if( !(i%8) && i>0 && i!=(ipnt-1) )
			{
				fprintf( ioQQQ, "\n  " );
			}
		}
		fprintf( ioQQQ, "\n" );
	}

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

/* =============================================================================*/
/* HeatZero is called by ConvBase */
void HeatZero( void )
{
	long int i , j;

	for( i=0; i < LIMELM; i++ )
	{
		for( j=0; j < LIMELM; j++ )
		{
			thermal.heating[i][j] = 0.;
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->HeatZero()\n", debug_fp );
#	endif
	return;
}
/*lint +e661 possible creation of out of bounds pointer */
/*lint +e662 possible creation of out of bounds pointer */
/*lint +e771 out of bounds pointer */
#undef PRT_DERIV
