/* 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 */
/*ConvTempEdenIoniz determine  temperature, called by ConPresTempEdenIoniz,
 * calls ConvEdenIoniz to get electron density and ionization */
/*MakeDeriv derive numerical derivative of heating minus cooling */
/*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
/*lgConvTemp returns true if heating-cooling is converged */
/*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
#include "cddefines.h"
#include "hmi.h"
#include "thermal.h"
#include "iso.h"
#include "hydrogenic.h"
#include "colden.h"
#include "h2.h"
#include "stopcalc.h"
#include "pressure.h"
#include "dense.h"
#include "trace.h"
#include "phycon.h"
#include "conv.h"

/* >>chng 01 apr 06, from 10 to 20, also put damper at 0.5 below,
 * as part of strategy of letting code bracket solution */
/* >>chng 04 mar 17, from 20 to 40, now doing much finer changes in
 * many conditions, including smaller dT, so allow more loops */
/* >>chng 05 mar 25, from 40 to 60, thermal fronts */
#define	LIMDEF	60

/*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
static double CoolHeatError( double temp );

/* use brent's method to find heating-cooling match */
static double TeBrent(
	  double x1, 
	  double x2);

/*MakeDeriv derive numerical derivative of heating minus cooling */
static void MakeDeriv(
  /* which job to do, either "zero" or "incr" */
  const char *job, 
  /* result, the numerical derivative */
  double *DerivNumer);

/*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
static void PutHetCol(double te, 
  double htot, 
  double ctot);

/*ConvTempEdenIoniz determine  temperature, called by ConPresTempEdenIoniz,
 * calls ConvEdenIoniz to get electron density and ionization 
 * returns 0 if ok, 1 if disaster */
int ConvTempEdenIoniz( void )
{
	char chDEType;
	long int limtry,
		nTemperature_loop ;
	/* followin is zero because there are paths were not set following abort */
	double CoolMHeat=0., 
	  Damper, 
	  dtl, 
	  fneut;
	int kase;
	double DerivNumer;
	static double OldCmHdT = 0.;
	long int nTemp_oscil;

	/* derivative of cooling less heating wrt temperature */
	double dCoolmHeatdT;

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

	/* determine ionization and temperature, called by PIonte
	 * nCall tells routine which call this is, 0 for first one, not changed here
	 * nTemperature_loop counts how many loops performed */

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
	}
	if( trace.lgTrConvg>=2 )
	{
		fprintf( ioQQQ, "  ConvTempEdenIoniz1 called. Te:%.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n", 
			phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
			100./MAX2(1e-36,thermal.htot), dense.eden );
	}

	if ( strcmp( conv.chSolverTemp , "new" ) == 0 )
	{
		/* this has never worked */
		double t1 , t2 ,
			error1 , error2;
		double TeNew;
		double factor = 1.02;

		t1 = phycon.te;
		error1 = CoolHeatError( t1 );
		t2 = t1 * factor;
		error2 = CoolHeatError( t2 );
		TeNew = TeBrent( t1 , t2 );
	}

	else if( strcmp( conv.chSolverTemp , "simple") == 0 )
	{

		/* main routine to find ionization and temperature
		 * following is flag to check for temp oscillations
		 * it could be set true in main temp loop */
		conv.lgTOscl = FALSE;
		/* ots rates not oscillating */
		conv.lgOscilOTS = FALSE;

		/* scale factor for changes in temp - make smaller if oscillations */
		Damper = 1.;

		/* flag for d(cool - heat)/dT changing sign, an internal error
		 * >>chng 97 sep 03, only turn false one time, remains true if ever true */
		if( nzone < 1 )
		{
			conv.lgCmHOsc = FALSE;
		}

		/* option to make numerical derivatives of heating cooling */
		MakeDeriv("zero",&DerivNumer);

		/* this is will initial limit to number of tries at temp */
		limtry = LIMDEF;

		/* used to keep track of sign of dt, to check for oscillations */
		dtl = 0.;

		/* this counts how may thermal loops we go through
		 * used in calling routine to make sure that at least
		 * two solutions are performed */
		nTemperature_loop = 0;
		nTemp_oscil = 0;

		/* this is change in temp, which is used early in following loop, must be preset */
		thermal.dTemper = 1e-3f*phycon.te;

		/* this is the start of the heating-cooling match loop,
		 * this exits by returning in the middle */
		while ( TRUE )
		{

			/* >>chng 02 dec 03 rjrw, insert these so dynamics cooling rate is calculated
			 *	 before step is chosen -- but how much duplication results? 
			 * must have current estimate of heating and cooling before temperature is changed.
			 * at bottom of following loop these two are called again to see whether heatin cooling
			 * have converged */
			if( ConvEdenIoniz() ) 
			{
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvTempEdenIoniz()\n", debug_fp );
#				endif
				return 1;
			}
			PresTotCurrent();

			/* we will try to make the following zero */
			CoolMHeat = thermal.ctot - thermal.htot;
			/*fprintf(ioQQQ,"DEBUG grn 13 14\t%e\t%e\t%e\n",
				phycon.te, 
				thermal.heating[0][13],
				thermal.heating[0][14] );*/

			if( thermal.lgTSetOn )
			{
				/* constant temp - declare this a success */
				CoolMHeat = 0.;
			}

			if( thermal.lgTLaw )
			{
				/* temperature law has been specified */
				if( thermal.lgTeBD96 )
				{
					/* Bertoldi & Drain 96 temp law specified by TLAW BD96 command */
					phycon.te = thermal.T0BD96 / (1.f + thermal.SigmaBD96 * colden.colden[ipCOL_HTOT] );
				}
				else if( thermal.lgTeSN99 )
				{
					/* Sternberg & Neufeld 99 temp law specified by TLAW SN99 command,
					 * this is equation 16 of 
					 * >>refer	H2	temperature law	Sternberg, A., & Neufeld, D.A. 1999, ApJ, 516, 371-380 */
					phycon.te = thermal.T0SN99 / 
						/*(float)(1.f + 9.*pow(2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 4.) );*/
						(float)(1.f + 9.*pow(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 4.) );
				}
				else
					TotalInsanity();
				tfidle(FALSE);
			}

			/* >>chng 02 dec 09, this block moved to just after loop entry above
			 * to check on heating cooling and
			 * exit if match is ok - prevents unncessary reevaluate of ConvEdenIoniz and PresTotCurrent */
			/* this is test for valid thermal solution,
			 * lgConvTemp returns TRUE if diff in heat cool is small */
			/* >>chng 02 dec 04, add check on size of dTemper rel to te,
			 * can become too small to increment a float */
			/* >>chng 04 dec 13, for very cold gas the dTemper < 1e-6 Te limit was hit
			 * during normal search.  dont let this limit be tested here, rather use
			 * fltepsilon to set limit to dTemper */
			if( (lgConvTemp() /*|| fabs(thermal.dTemper) < 1e-6*phycon.te*/) 
					&& conv.lgConvIoniz )
			{
				/* we have a good solution */
				if( trace.lgTrace )
					fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );

				/* test for thermal stability, set flag if unstable, 
				 * later increment number of zones only once per zone */
				/* this may become a check on thermal stability since neg slope is
				* unstable, but we may not have a valid soln, so only set flag */
				if( thermal.dCooldT < 0. )
				{
					thermal.lgUnstable = TRUE;
				}
				else
				{
					thermal.lgUnstable = FALSE;
				}

				/* remember the highest and lowest temperatures */
				thermal.thist = (float)MAX2(phycon.te,thermal.thist);
				thermal.tlowst = (float)MIN2(phycon.te,thermal.tlowst);

				/* say that we have a valid solution */
				conv.lgConvTemp = TRUE;

				if( trace.lgTrConvg>=2 )
				{
					fprintf( ioQQQ, "  ConvTempEdenIoniz2 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n", 
					  phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
					  100./MAX2(1e-36,thermal.htot), dense.eden );
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvTempEdenIoniz()\n", debug_fp );
#				endif
				/*******************************************************
				 * 
				 * this is return for valid solution 
				 *
				 *******************************************************/
				return 0 ;
			}
			else if (nTemperature_loop >= limtry || lgAbort )
			{
				/* >>chng 02 dec 17, add test on ots fields oscillating
				 * while ionization is not converged */
				/* >>chng 04 may 25, do not test on lgOscilOTS since we do test
				 * on lgConvIoniz */
				if( nTemperature_loop == LIMDEF && !conv.lgTOscl && 
					conv.lgConvIoniz /*&& !conv.lgOscilOTS*/ )
				{
					/* limdef is set to 40 in a macro above
					 * >>chng 97 aug 10, from 2 to 4 to let keep going when not oscillation
					 * this helped in getting over huge jumpts in temperature at metal fronts */
					limtry = LIMDEF*4;
					if( trace.lgTrConvg>2 )
					{
						fprintf( ioQQQ, "  ConvTempEdenIoniz9 increases limtry to%3ld\n", 
										 limtry );
					}
				}
				else
				{
					/* looped too long, so bail out anyway */
#					ifdef DEBUG_FUN
					fputs( " <->ConvTempEdenIoniz()\n", debug_fp );
#					endif
					break;
				}
			}

			/* get numerical heating-cooling derivative, but we may not use it */
			MakeDeriv("incr",&DerivNumer);

			/* remember this temp, heating, and cooling */
			PutHetCol(phycon.te,thermal.htot,thermal.ctot);

			/* the flag conv.lgConvEden was set during above call to ConvEdenIoniz, and
			 * if is not true then eden failed, and we abort, 
			 * but don't abort on very first try, when flag has not been set 
			 * >> chng 02 dec 10 rjrw flag is now set on first try */
			if( !conv.lgConvEden )
			{
				/* set this to zero so that will not flag temperature failure, since we are
				 * aborting due to eden failure */
				CoolMHeat = 0.;
				if( trace.lgTrConvg>=2 )
				{
					fprintf( ioQQQ, 
						"  ConvTempEdenIoniz3 breaks since lgConvEden returns failed electron density.\n"); 
				}
				break;
			}

			/* save old deriv to check for oscillations */
			if( nzone < 1 )
			{
				/* >>chng 96 jun 07, was set to dCoolmHeatdT, is not defined initially */
				OldCmHdT = thermal.ctot/phycon.te;
			}

#			define USENUMER FALSE
			if( USENUMER )
			{
				dCoolmHeatdT = DerivNumer;
				chDEType = 'n';
			}
			else if( phycon.te < 1e4 && thermal.dCooldT < 0. )
			{
				/* use numerical guess, HTOT nearly equal to CTOT, drive Te lower
				 * this is to overcome thermal inflection points */
				dCoolmHeatdT = (double)(thermal.htot)/phycon.te*2.;
				chDEType = 'd';
			}
			else
			{
				/* >>chng 97 mar 18, all below growth code just did the following 
				 * use analytical estimate of change in heat - cooling wrt temp */
				dCoolmHeatdT = thermal.dCooldT - thermal.dHeatdT;
				chDEType = 'a';
			}

			/* check whether deriv changing sign - an internal error or
			 * numerical instability */
			if( OldCmHdT*dCoolmHeatdT < 0. )
			{
				conv.lgCmHOsc = TRUE;

				/* derivative can change sign when heating and cooling balance,
				* and processes have exactly same temp dependence
				* happened in mods of ism.in test
				* >>chng 97 sep 02, added following test for oscillating derivs */
				dCoolmHeatdT = thermal.ctot/phycon.te;
				chDEType = 'o';
			}
			OldCmHdT = dCoolmHeatdT;

			/* if not first pass through then numerical deriv should not be zero */
			/* required change in temperature, this may be too large */
			thermal.dTemper = (float)(CoolMHeat/dCoolmHeatdT);
			kase = 0;

			/* try to stop numerical oscillation if dTemper is changing sign */
			/* when first time through and far from soln (when zones a bit too thick) first
			 * step can be in wrong direction, and all subsequent are in right direciton.
			 * do not call this an oscillation - only declare oscillation when second occurs,
			 * nTemperature_loop is zero first time through, (dtl is also zero),
			 * is 1 first time through with memory of this zone, (when false oscillation occurs)
			 * and is >1 thereafter */
			/* >>chng 05 mar 25, had counted as oscil when loop count greater than 1, change to 3
			 * to allow things to settle down.  Also, for first three evaluations, do not let te
			 * chamge by too much, until lower solvers settle down */
			if( dtl*thermal.dTemper < 0. && nTemperature_loop > 3)
			{
				conv.lgTOscl = TRUE;
				/* make DT smaller and smaller */
				/* >>chng 01 mar 16, from 0.5 to 0.75 
				Damper *= 0.75;*/
				Damper *= 0.5;
				nTemp_oscil = nTemperature_loop;
			}
			else if( Damper < 1. && nTemperature_loop-nTemp_oscil > 10 )
			{
				/* >>chng 05 mar 26, some cases there is noise in heating-cooling
				 * in first few evaluations and this causes te oscillations.  but it may
				 * settle down and become well behaved.  this allows that to happen */
				conv.lgTOscl = FALSE;
				if( !((nTemperature_loop-nTemp_oscil)%5) )
					Damper = MIN(1., Damper*1.5 );
			}

			/* SIGN: sign of sec arg and abs val of first */
			/* >>chng 96 nov 08, from 0.1 to 0.08 in fneut */
			/*fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.Hmolec[ipMH2g])/dense.gas_phase[ipHYDROGEN];*/
			fneut = (dense.xIonDense[ipHYDROGEN][0] + 2.*hmi.H2_total)/dense.gas_phase[ipHYDROGEN];
			/*>>chng 04 jan 11, when big H2 is on and heating is important,
			 * only allow small changes in temp since must converge pops */
			/* >>chng 04 jan 15, chng ratio of heating in h2 to total as the criteria,
			 * to the error.  temp changes will become much larger when this is
			 * not tripped, and h2 can start to oscillate wildly. 
			 * NB this limit should be kept parallel with similar test for convergence
			 * of heating in h2 */
			if( 0 && h2.lgH2ON && fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > conv.HeatCoolRelErrorAllowed )
			{
				/* heating due to collisional deexcitation of the large H2 molecule
				 * is very important.  this rate depends strongly on the temperature
				 * due to Boltzmann factors, and large changes in temperature
				 * inhibit convergence of the level pops in the big H2.  That
				 * atom has tests on converging the heating, so want only small
				 * temperature changes when H2 heating is important. */
				double absdt = fabs(thermal.dTemper);
				thermal.dTemper = (sign(MIN2(absdt,0.0025*phycon.te),
				  thermal.dTemper));
				kase = 1;
			}
			else if( fneut > 0.08 && hydro.lgHColionImp )
			{
				/* lgHColionImp is true if model collisionally ionized
				 * dont let TE change by too much if lgHColionImp is set by HLEVEL */
				double absdt = fabs(thermal.dTemper);
				thermal.dTemper = (sign(MIN2(absdt,0.005*phycon.te),
				  thermal.dTemper));
				kase = 2;
			}
#			if 0
			/* >>chng 05 mar 25, for first few steps through solvers in a new zone, don't let temp
			 * change by too much initially, to avoid a te oscil which sets the oscil flag */
			else if( 0 && nTemperature_loop <= 3)
			{
				/* >>chng 05 mar 25, add this branch, an early pass through these solvers, 
				 * do not let te change by too much to avoid te oscillations */
				double absdt = fabs(thermal.dTemper);
				thermal.dTemper = (sign(MIN2(absdt,0.002*phycon.te),
				  thermal.dTemper));
				kase = 3;
			}
#			endif
			else
			{
				/* this branch, we will use the change in temperature that came from the deriv,
				 * but do not let temp change by more than 2 percent */
				double absdt = fabs(thermal.dTemper);
				thermal.dTemper = (sign(MIN2(absdt,0.02*phycon.te),
				  thermal.dTemper));
				kase = 4;
			}

			/* >>chng 04 mar 14, add this test, cooling flow clouds had temp failure
			 * when dt as large as 4% were allowed with the big H2 */
			/* in no case do we want the temperature to change by more than 1% when the big
			 * H2 molecule is turned on - it needs small temp changes to converge */
			if( 0 && h2.lgH2ON && fabs(thermal.dTemper)/ phycon.te > 0.01 )
			{
				double absdt = fabs(thermal.dTemper);
				thermal.dTemper = (sign(MIN2(absdt,0.01*phycon.te),
				  thermal.dTemper));
				kase = 5;
			}

			/* >>chng 03 mar 16, to following logic, previous had been impossible
			 * to hit - only do fine changes in Te when at half-way point in i-front
			 * or if model is in collisional recombination equilibirum
			 * if so then small change in temp causees big change in equilibrium
			 * lgHColRec set in hlevel, only true is totl rec >> rad rec */
			/* are most recombinations due to 3-body? */
			if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.9 &&
				/* is hydrogen more than 30% ionized */
				dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.01 )
			{
				if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 )
					/* >>chng 04 feb 07, add test for H mostly collisionally ionized */
					/* is hydrogen less than 70% ionized? */
					/* >>chng 04 feb 07, from 70% to 90% ionized - oscill in hard
					* x-ray continua (laser at 80 ryd) when i front hit */
				{
					double absdt = fabs(thermal.dTemper);
					/* >>chng 03 nov 03 , fac from 0.001 to 0.004 */
					/* >>chng 04 feb 24 , fac back from 0.004 to 0.001,
					 * this is needed for col den 25, high density (14) blr models
					 * with flux of 22 - small changes in temperature drive large
					 * changes in electron density due to three body recomb
					 * this small change is ok due to test above, which makes sure
					 * that we are in collsiional recombination limit 
					 * sign returns sign of sec par times value of first */
					thermal.dTemper = (sign(MIN2(absdt,0.001*
						phycon.te),thermal.dTemper));
					kase = 6;
				}
			}

			/* make steps smaller if oscillations, caused by overshots, occur */
			thermal.dTemper *= Damper;

			/* do not let dTemper/Te grow smaller than 10 * flt epsilon, the smaller number that
			 * can be added to a float and still work */
			if( fabs(thermal.dTemper)/phycon.te <10.*DBL_EPSILON )
			{
				kase = 7;
				thermal.dTemper = (sign( 10.*DBL_EPSILON*phycon.te , thermal.dTemper));
			}

			/* save last change in temperature */
			dtl = thermal.dTemper;

			/* THIS IS IT !!! this is it !!! this is where te changes. */
			if( !thermal.lgTLaw )
			{
				/* usual case - valid thermal solution */
				phycon.te = phycon.te - thermal.dTemper;
				tfidle(FALSE);
			}

			if( trace.lgTrConvg>=2 )
			{
				fprintf( ioQQQ, 
					"  ConvTempEdenIoniz4 a %4li T:%.4e dt:%10.2e%7.2f%% C:%10.2e H:%10.2e"
					" CoolMHeat/h:%7.2f%% dC-H/dT:%.2e kase:%i%c damp%.2f\n", 
				  nTemperature_loop, 
				  phycon.te, 
				  thermal.dTemper, 
				  thermal.dTemper/(phycon.te+thermal.dTemper)*100,
				  thermal.ctot, 
				  thermal.htot, 
				  CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. ,
				  dCoolmHeatdT,
				  kase,
				  /* which type of deriv we have */
				  chDEType,
				  Damper);
				/* dCoolmHeatdT is the actual number used for getting dT
				 * thermal.dCooldT is just cooling */
			}

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, 
					" ConvTempEdenIoniz TE:%.5e dT%.3e HTOT:%.5e CTOT:%.5e dCoolmHeatdT:%c%.4e C-H%.4e HDEN%.2e kase %i\n", 
				  phycon.te, 
				  thermal.dTemper, 
				  thermal.htot, 
				  thermal.ctot, 
				  chDEType, 
				  dCoolmHeatdT, 
				  CoolMHeat, 
				  dense.gas_phase[ipHYDROGEN],
				  kase );
			}

#			if 0

			/* >>chng 02 jun 11, move this call down to here, from above,
			 * so that temperature is changed, and heating cooling evaluated then */
			/* temperature has changed, so we need to reconverge the ionization and electron density; 
			 * this calls ionize until lgIonDone is true,
			 * heating and cooling are reevaluated in ConvBase 
			 * at bottom of this call chain */
			if( ConvEdenIoniz() )
			{
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvTempEdenIoniz()\n", debug_fp );
#				endif
				return 1;
			}

			/* must reevaluate pressure since temperature changed */
			/* this sets values of pressure.PresTotlCurr, also calls tfidle */
			PresTotCurrent();

			CoolMHeat = thermal.ctot - thermal.htot;

			/* >>chng 02 dec 09, this block moved to just after loop entry above
			 * to check on heating cooling and
			 * exit if match is ok - prevents unncessary reevaluate of ConvEdenIoniz and PresTotCurrent */
			/* if te falls below 3K quit now */
			if( phycon.te <= 3. )
			{
				/* these will force all routines out to quit
				 * when lowest temp read in values below 3 were not allowed */
				phycon.te = 2.9998f;
				tfidle(FALSE);
				/* this sets values of pressure.PresTotlCurr, also calls tfidle,
				 * must do because we just changed the temperature */
				PresTotCurrent();
				CoolMHeat = 0.;
				conv.lgConvIoniz = TRUE;
			}

			if( trace.lgTrConvg>=2 )
			{
				fprintf( ioQQQ, 
					"  ConvTempEdenIoniz5 b %4li T:%.4e dt:%9.2e%7.1f%% C:%10.2e H:%10.2e CoolMHeat/h:%7.1f%% dC-H/dT:%.2e kase:%i damp:%.2f\n", 
				  nTemperature_loop, 
				  phycon.te, 
				  thermal.dTemper, 
				  thermal.dTemper/(phycon.te+thermal.dTemper)*100,
				  thermal.ctot, 
				  thermal.htot, 
				  CoolMHeat/MIN2(thermal.ctot , thermal.htot )*100. ,
				  dCoolmHeatdT,
				  kase ,
				  Damper);
				/* dCoolmHeatdT is the actual number used for getting dT
				 * thermal.dCooldT is just cooling */
			}

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, 
					" ConvTempEdenIoniz TE:%11.5e dT%10.3e HTOT:%11.5e CTOT:%11.5e dCoolmHeatdT:%c%.4e C-H%.4e kase %li\n", 
				  phycon.te, thermal.dTemper, thermal.htot, thermal.ctot, chDEType, 
				  dCoolmHeatdT, CoolMHeat, 
				  kase );
			}

			/* this is test for valid thermal solution,
			 * lgConvTemp returns TRUE if diff in heat cool is small */
			/* >>chng 02 dec 04, add check on size of dTemper rel to te,
			 * can become too small to increment a float */
			/* if te falls below 3K quit now */
			if( phycon.te <= 3. ||
				((lgConvTemp() || fabs(thermal.dTemper) < 1e-6*phycon.te) 
				&& conv.lgConvIoniz) )
			{
				/* we have a good solution */
				if( trace.lgTrace )
					fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );

				/* remember the highest and lowest temperatures */
				thermal.thist = (float)MAX2(phycon.te,thermal.thist);
				thermal.tlowst = (float)MIN2(phycon.te,thermal.tlowst);

				/* test for thermal stability, set flag if unstable, 
				 * later increment number of zones only once per zone */
				if( thermal.dCooldT < 0. )
				{
					thermal.lgUnstable = TRUE;
				}
				else
				{
					thermal.lgUnstable = FALSE;
				}

				/* say that we have a valid solution */
				conv.lgConvTemp = TRUE;

				if( trace.lgTrConvg>=2 )
				{
					fprintf( ioQQQ, "  ConvTempEdenIoniz6 converged. Te:%11.4e Htot:%11.4e Ctot:%11.4e error:%9.1e%%, eden=%11.4e\n", 
					  phycon.te, thermal.htot, thermal.ctot, (thermal.htot - thermal.ctot)*
					  100./MAX2(1e-36,thermal.htot), dense.eden );
				}
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvTempEdenIoniz()\n", debug_fp );
#				endif
				/*******************************************************
				 * 
				 * this is return for valid solution 
				 *
				 *******************************************************/
				return 0;
			}
#			endif

			/* increment loop counter */
			++nTemperature_loop;
		}

		/* fall through, exceeded limit to number of solutions,
		 * >>chng 01 mar 14, or no electron density convergence
		 * no temperature convergence */
		if( fabs(CoolMHeat/thermal.htot ) > conv.HeatCoolRelErrorAllowed )
		{
			conv.lgConvTemp = FALSE;
			if( trace.lgTrConvg>=2 )
			{
				fprintf( ioQQQ, "  ConvTempEdenIoniz7 fell through loop, heating cooling not converged, setting conv.lgConvTemp = FALSE\n");
			}
		}
		else
		{
			/* possible that ionization has failed, and we fell through to here
			 * with a valid thermal soln */
			conv.lgConvTemp = TRUE;
			if( trace.lgTrConvg>=2 )
			{
				fprintf( ioQQQ, "  ConvTempEdenIoniz8 fell through loop, heating cooling is converged (ion not?), setting conv.lgConvTemp = TRUE\n");
			}
		}
	}
	else
	{
		fprintf( ioQQQ, "ConvTempEdenIoniz finds insane solver%s \n" , conv.chSolverTemp );
		ShowMe();
		puts( "[Stop in ConvTempEdenIoniz]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*MakeDeriv derive numerical derivative of heating minus cooling */
static void MakeDeriv(
  const char *job, 
  double *DerivNumer)
{
	static long int nstore;
	static double OldTe , OldCool , OldHeat;

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

	if( strcmp(job,"zero") == 0 )
	{
		nstore = 0;
		OldTe = 0.;
		OldCool = 0.;
		OldHeat = 0.;
	}

	else if( strcmp(job,"incr") == 0 )
	{

		if( nstore > 0 && !thermal.lgTSetOn && fabs(phycon.te - OldTe)>SMALLFLOAT )
		{
			/* get numerical derivatives */
			double dCdT = (thermal.ctot - OldCool)/(phycon.te - OldTe);
			double dHdT = (thermal.htot - OldHeat)/(phycon.te - OldTe);
			*DerivNumer = dCdT - dHdT;
#			if 0
			fprintf(ioQQQ,"derivderiv\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
			  nzone, dCdT , thermal.dCooldT , 
			  thermal.ctot , OldCool,phycon.te , OldTe
			  /*(thermal.ctot - OldCool),(phycon.te - OldTe),dHdT, thermal.dHeatdT*/);
#			endif
		}

		else
		{
			/* if early in soln, or constant temperature, return zero */
			*DerivNumer = 0.;
		}

		OldTe = phycon.te;
		OldCool = thermal.ctot;
		OldHeat = thermal.htot;
		++ nstore;

	}
	else
	{
		fprintf( ioQQQ, " MakeDeriv called with insane job=%4.4s\n", 
		  job );
		puts( "[Stop in makederiv]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*PutHetCol save heating, cooling, and temperature in stack for numerical derivatives */
static void PutHetCol(double te, 
  double htot, 
  double ctot)
{

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

	if( nzone == 0 )
	{
		/* code is searching for first temp now - do not remember these numbers */
		thermal.ipGrid = 0;
	}
	else
	{
		if( thermal.ipGrid >= NGRID )
			thermal.ipGrid = 0;

		ASSERT( thermal.ipGrid >= 0 );
		ASSERT( thermal.ipGrid < NGRID );
		ASSERT( nzone>=0 );

		thermal.TeGrid[thermal.ipGrid] = (float)te;
		thermal.HtGrid[thermal.ipGrid] = (float)htot;
		thermal.ClGrid[thermal.ipGrid] = (float)ctot;
		thermal.nZonGrid[thermal.ipGrid] = nzone;

		++thermal.ipGrid;
	}

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

/*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
static double CoolHeatError( double temp )
{

	phycon.te = temp;
	tfidle(FALSE);

	/* converge the ionization and electron density; 
	 * this calls ionize until lgIonDone is true */
	/* NB should NOT set insanity - but rather return error condition */
	if( ConvEdenIoniz() )
		lgAbort = TRUE;

	/* >>chng 01 mar 16, evaluate pressure here since changing and other values needed */
	/* reevaluate pressure */
	/* this sets values of pressure.PresTotlCurr, also calls tfidle */
	PresTotCurrent();

	/* remember this temp, heating, and cooling */
	PutHetCol(phycon.te,thermal.htot,thermal.ctot);

	return thermal.ctot - thermal.htot;
}

/*#define	EPS	3.e-8*/
#define	ITERMAX	100
/* use brent's method to find heating-cooling match */
static double TeBrent(
	  double x1, 
	  double x2)
{
	long int iter;
	double c, 
	  d, 
	  e, 
	  fx1, 
	  fx2, 
	  fc, 
	  p, 
	  q, 
	  r, 
	  s, 
	  xm;

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

	/* first evaluate function at two boundaries, and check
	 * that we have bracketed the target */
	/*NB big logical error - this routine ignores return value of
	 * routine it calls, so will continue despite disaster */
	fx1 = CoolHeatError(x1);
	fx2 = CoolHeatError(x2);

	if( fx1*fx2 > 0. )
	{
		/* both values either positive or negative, this is insane */
		fprintf( ioQQQ, " TeBrent called without proper bracket - this is impossible\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in TeBrent]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* we have bracketed the target */
	c = x2;
	fc = fx2;
	iter = 0;
	/* these are sanity checks, to get code past lint */
	e = DBL_MAX;
	d = DBL_MAX;
	while(  iter < ITERMAX )
	{
		if( (fx2 > 0. && fc > 0.) || (fx2 < 0. && fc < 0.) )
		{
			c = x1;
			fc = fx1;
			d = x2 - x1;
			e = d;
		}
		if( fabs(fc) < fabs(fx2) )
		{
			x1 = x2;
			x2 = c;
			c = x1;
			fx1 = fx2;
			fx2 = fc;
			fc = fx1;
		}

		xm = .5*(c - x2);

		/* soln converged if residual less than tol, or hit zero */
		if( fabs(xm) <= thermal.htot*conv.HeatCoolRelErrorAllowed || fx2 == 0. )
			break;

		if( fabs(e) >= thermal.htot*conv.HeatCoolRelErrorAllowed && fabs(fx1) > fabs(fx2) )
		{
			double aa , bb ;
			s = fx2/fx1;
			/*lint -e777 test floats for equality */
			if( x1 == c )
			/*lint +e777 test floats for equality */
			{
				p = 2.*xm*s;
				q = 1. - s;
			}
			else
			{
				q = fx1/fc;
				r = fx2/fc;
				p = s*(2.*xm*q*(q - r) - (x2 - x1)*(r - 1.));
				q = (q - 1.)*(r - 1.)*(s - 1.);
			}
			if( p > 0. )
				q = -q;

			p = fabs(p);
			aa = fabs(thermal.htot*conv.HeatCoolRelErrorAllowed*q);
			bb = fabs(e*q);
			if( 2.*p < MIN2(3.*xm*q-aa,bb) )
			{
				e = d;
				d = p/q;
			}
			else
			{
				d = xm;
				e = d;
			}
		}
		else
		{
			d = xm;
			e = d;
		}
		x1 = x2;
		fx1 = fx2;
		if( fabs(d) > thermal.htot*conv.HeatCoolRelErrorAllowed )
		{
			x2 += d;
		}
		else
		{
			x2 += sign(thermal.htot*conv.HeatCoolRelErrorAllowed,xm);
		}
		fx2 = CoolHeatError(x2);

		++iter;
	}

	/* check whether we fell through (hit limit to number of iterations)
	 * of broke out when complete */
	if( iter == ITERMAX )
	{
		/* both values either positive or negative, this is insane */
		fprintf( ioQQQ, " TeBrent did not converge in %i iterations\n",ITERMAX );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in TeBrent]" );
		cdEXIT(EXIT_FAILURE);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->TeBrent()\n", debug_fp );
#	endif
	return( x2 );

}

/* returns true if heating-cooling is converged */
int lgConvTemp(void)
{
	int lgRet;
	/* >>chng 03 sep 19, add check on StopCalc.TeLowest */
	if( fabs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed ||
		thermal.lgTSetOn || phycon.te <= StopCalc.TeLowest )
	{
		/* announce that temp is converged if relative heating - cooling mismatch
		 * is less than the relative heating cooling error allowed, or
		 * if this is a constant temperature model */
		lgRet = TRUE;
		conv.lgConvTemp = TRUE;
	}
	else
	{
		/* big mismatch, this has not converged */
		lgRet = FALSE;
		conv.lgConvTemp = FALSE;
	}
	return lgRet;
}
#undef	LIMDEF

