/home66/gary/public_html/cloudy/c08_branch/source/conv_temp_eden_ioniz.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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 "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 */
static const int 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;
        /* following 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;

        DEBUG_ENTRY( "ConvTempEdenIoniz()" );

        /* 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.nTrConvg>=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 );
                fprintf(ioQQQ , "DEBUG new temp solver error1 %.2e error2 %.2e TeNew %.2e\n",
                        error1 , error2 , TeNew );
        }

        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() ) 
                        {
                                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.lgTemperatureConstant )
                        {
                                /* constant temp - declare this a success */
                                CoolMHeat = 0.;
                        }

                        if( thermal.lgTLaw )
                        {
                                double TeNew = phycon.te;
                                /* temperature law has been specified */
                                if( thermal.lgTeBD96 )
                                {
                                        /* Bertoldi & Drain 96 temp law specified by TLAW BD96 command */
                                        TeNew = 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 */
                                        TeNew = thermal.T0SN99 / 
                                                /*(realnum)(1.f + 9.*pow(2.*hmi.Hmolec[ipMH2g]/dense.gas_phase[ipHYDROGEN], 4.) );*/
                                                (realnum)(1.f + 9.*pow(2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN], 4.) );
                                }
                                else
                                        TotalInsanity();

                                TempChange(TeNew ,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 unnecessary reevaluate of ConvEdenIoniz and PresTotCurrent */
                        /* this is test for valid thermal solution,
                         * lgConvTemp returns true if difference 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.  do not 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 solution, so only set flag */
                                if( thermal.dCooldT < 0. )
                                {
                                        thermal.lgUnstable = true;
                                }
                                else
                                {
                                        thermal.lgUnstable = false;
                                }

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

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

                                if( trace.nTrConvg>=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 );
                                }
                                /*******************************************************
                                 * 
                                 * 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 60 above
                                         * >>chng 97 aug 10, from 2 to 4 to let keep going when not oscillation
                                         * this helped in getting over huge jumps in temperature at metal fronts */
                                        limtry = LIMDEF*4;
                                        if( trace.nTrConvg>2 )
                                        {
                                                fprintf( ioQQQ, "  ConvTempEdenIoniz9 increases limtry to%3ld\n", 
                                                                                 limtry );
                                        }
                                }
                                else
                                {
                                        /* looped too long, so bail out anyway */
                                        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.nTrConvg>=2 )
                                {
                                        fprintf( ioQQQ, 
                                                "  ConvTempEdenIoniz3 breaks since lgConvEden returns failed electron density.\n"); 
                                }
                                break;
                        }

                        /* save old derivative 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 derivative 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 derivative should 
                         * not be zero required change in temperature, this may be too large */
                        thermal.dTemper = (realnum)(CoolMHeat/dCoolmHeatdT);
                        kase = 0;

                        /* try to stop numerical oscillation if dTemper is changing sign */
                        /* when first time through and far from solution (when zones a bit too thick) first
                         * step can be in wrong direction, and all subsequent are in right direction.
                         * 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 oscillations when loop count greater than 1, change to 3
                         * to allow things to settle down.  Also, for first three evaluations, do not let te
                         * change 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 = MIN2( 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
                                 * do not 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 equilibrium
                         * if so then small change in temp causes 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 collisional 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 */
                                TempChange(phycon.te - thermal.dTemper , false);
                        }

                        if( trace.nTrConvg>=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 );
                        }

                        /* 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.nTrConvg>=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 solution */
                        conv.lgConvTemp = true;
                        if( trace.nTrConvg>=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();
                cdEXIT(EXIT_FAILURE);
        }
        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;

        DEBUG_ENTRY( "MakeDeriv()" );

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

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

                if( nstore > 0 && !thermal.lgTemperatureConstant && 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 );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

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

        DEBUG_ENTRY( "PutHetCol()" );

        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] = (realnum)te;
                thermal.HtGrid[thermal.ipGrid] = (realnum)htot;
                thermal.ClGrid[thermal.ipGrid] = (realnum)ctot;
                thermal.nZonGrid[thermal.ipGrid] = nzone;

                ++thermal.ipGrid;
        }
        return;
}

/*CoolHeatError evaluate ionization, and difference in heating and cooling, for temperature temp */
STATIC double CoolHeatError( double temp )
{
        DEBUG_ENTRY( "CoolHeatError()" );

        TempChange(temp , 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 */
        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;

        DEBUG_ENTRY( "TeBrent()" );

        /* 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" );
                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);

                /* solution converged if residual less than tolerance, 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;
                        if( fp_equal( x1, c ) )
                        {
                                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" );
                cdEXIT(EXIT_FAILURE);
        }
        return( x2 );

}

/* returns true if heating-cooling is converged */
bool lgConvTemp(void)
{
        bool lgRet;
        DEBUG_ENTRY( "lgConvTemp()" );

        /* >>chng 03 sep 19, add check on phycon.TEMP_LIMIT_LOW */
        if( fabs(thermal.htot - thermal.ctot)/thermal.htot <= conv.HeatCoolRelErrorAllowed ||
                thermal.lgTemperatureConstant || phycon.te <= phycon.TEMP_LIMIT_LOW )
        {
                /* 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;
}

---