/*ConvTempEdenIoniz determine ionization and temperature, called by ConPresTempEdenIoniz*/
#include <float.h>
#include "cddefines.h"
#include "sign.h"
#define	LIMDEF	10
#include "showme.h"
#include "hmi.h"
#include "tehigh.h"
#include "h.h"
#include "zonecnt.h"
#include "tfidle.h"
#include "boltgn.h"
#include "colion.h"
#include "tcool.h"
#include "tseton.h"
#include "called.h"
#include "trace.h"
#include "phycon.h"
#include "heat.h"
#include "toler.h"
#include "nustbl.h"
#include "dcneg.h"
#include "dtemp.h"
#include "hcolrec.h"
#include "bit32.h"
#include "holod.h"
#include "makederiv.h"
#include "puthetcol.h"
#include "converge.h"

void ConvTempEdenIoniz( void )
{
	char chDEType;
	long int limtry,
		nThLoop ;
	double CoolMHeat=-DBL_MAX, 
	  Damper, 
	  dtl, 
	  absdt ,
	  fneut;
	float DerivNumer;
	static double OldCmHdT = 0.;

#	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
	 * nThLoop counts how many loops performed */

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "\n ConvTempEdenIoniz called\n" );
	}

	/* there are some cases where this routine is called with a
	 * whole new temperature/density (start of new iterations for
	 * example).  These are called for protection and safely */

	/* update density - temperature variables which may have changed */
	tfidle();

	/* reevaluate gaunt factors and Boltzmann factors */
	boltgn();

	/* 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 3, only turn false one time, remains true if ever true */
	if( ZoneCnt.nzone < 1 )
	{
		conv.lgCmHOsc = FALSE;
	}

	/* option to make numerical derivatives of heating cooling */
	MakeDeriv("zero",0.,0.,&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 */
	nThLoop = 0;

	/* this is the start of the heating-cooling match loop */
	while( nThLoop < limtry )
	{

		/* converge the ionization and electron density; 
		 * this calls ionize until lgIonDone is true */
		ConvEdenIoniz();

		if( called.lgTalk )
		{
			if( heat.htot < 0. )
			{
				fprintf( ioQQQ, " Total heating is <0; is this model collisionally ionized?\n" );
			}
			else if( heat.htot == 0. )
			{
				fprintf( ioQQQ, " Total heating is 0; is the density small?\n" );
			}
		}

		fneut = (h.hi + 2.*hmi.htwo)/phycon.hden;

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

		if( tcool.dCooldT < 0. )
		{
			dcneg.lgdCNeg = TRUE;
		}
		else
		{
			dcneg.lgdCNeg = FALSE;
		}

		++nThLoop;

		if( tseton.lgTSetOn )
		{
			/* constant temp - declare this a success */
			CoolMHeat = 0.;
		}
		else
		{
			/* we will try to make the following zero */
			CoolMHeat = tcool.ctot - heat.htot;
		}

		/* get numerical derivative, but we may not use it */
		MakeDeriv("incr",tcool.ctot-heat.htot,phycon.te,&DerivNumer);

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

		if( phycon.te < 1e4 && tcool.dCooldT < 0. )
		{
			/* use numerical guess, HTOT nearly equal to CTOT, drive Te lower
			 * this is to overcome thermal inflection points */
			tcool.dCoolmHeatdT = (double)(heat.htot)/phycon.te*2.;
			chDEType = 'd';
		}
		else
		{
			/* >>chng 97 mar 18, all below growth code just did the following */
			tcool.dCoolmHeatdT = tcool.dCooldT - heat.dHTotDT;
			chDEType = 'a';
		}

		/* >>chng 96 Aug 15, alphas with fast compile underflow at low densities */
		if( tcool.dCoolmHeatdT == 0. )
		{
			chDEType = 'z';
			tcool.dCoolmHeatdT = tcool.ctot/phycon.te;
			if( tcool.dCoolmHeatdT == 0. && tcool.ctot == 0. )
			{
				fprintf( ioQQQ, " ConvTempEdenIoniz the cooling and its derivitive are zero.\n" );
				fprintf( ioQQQ, " ConvTempEdenIoniz I cannot perform the energy balance.\n" );
				if( phycon.hden <= 1e-3 && bit32.lgBit32 )
				{
					fprintf( ioQQQ, " Due to very low density on a 32-big cpu?\n" );
				}
				ShowMe();
				puts( "[Stop in ConvTempEdenIoniz]" );
				exit(1);
			}
			else if( tcool.dCoolmHeatdT == 0. && tcool.ctot != 0. )
			{
				fprintf( ioQQQ, " ConvTempEdenIoniz the cooling derivitive is zero.\n" );
				fprintf( ioQQQ, " ConvTempEdenIoniz I cannot perform the energy balance.\n" );
				if( phycon.hden <= 1e-3 && bit32.lgBit32 )
				{
					fprintf( ioQQQ, " Due to very low density on a 32-big cpu?\n" );
				}
				ShowMe();
				puts( "[Stop ConvTempEdenIoniz]" );
				exit(1);
			}
		}

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

		/* 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 */
		if( conv.lgCmHOsc )
		{
			tcool.dCoolmHeatdT = tcool.ctot/phycon.te;
			chDEType = 'o';
		}

		/* required change in temperature, this may be too large */
		dtemp.dTemper = (float)(CoolMHeat/tcool.dCoolmHeatdT);

		/* try to stop numerical oscillation if dTemper is changing sign */
		if( dtl*dtemp.dTemper < 0. )
		{
			conv.lgTOscl = TRUE;
			/* make DT smaller and smaller */
			Damper *= 0.5;
		}

		/* SIGN: sign of sec arg and abs val of first
		 * >>chng 96 nov 8, from 0.1 to 0.08 in fneut */
		if( fneut > 0.08 && colion.lgHColionImp )
		{
			/* lgHColionImp is true if model collisionally ionized
			 * dont let TE change by too much if lgHColionImp is set by HLEVEL */
			absdt = fabs(dtemp.dTemper),
			dtemp.dTemper = (float)(sign(MIN2(absdt,0.005*phycon.te),
			  dtemp.dTemper));
		}
		else if( nThLoop > 5 && !conv.lgTOscl )
		{
			/* need several changes but te not oscillating, increase dT */
			absdt = fabs(dtemp.dTemper),
			dtemp.dTemper = (float)(sign(MIN2(absdt,0.06*phycon.te),
			  dtemp.dTemper));
		}
		else if( nThLoop > 3 && !conv.lgTOscl )
		{
			/* need several changes but te not oscillating, increase dT */
			absdt = fabs(dtemp.dTemper),
			dtemp.dTemper = (float)(sign(MIN2(absdt,0.04*phycon.te),
			  dtemp.dTemper));
		}
		else
		{
			/* do not let temp change by more than 2 percent */
			absdt = fabs(dtemp.dTemper),
			dtemp.dTemper = (float)(sign(MIN2(absdt,0.02*phycon.te),
			  dtemp.dTemper));
		}

		/* check if model 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 */
		if( HColRec.lgHColRec )
		{
			if( HColRec.lgHCoolng )
			{
				if( conv.lgTOscl )
				{
					absdt = fabs(dtemp.dTemper),
					dtemp.dTemper = (float)(sign(MIN2(absdt,0.002*phycon.te),dtemp.dTemper));
				}
				else
				{
					/*   not oscillating */
					absdt = fabs(dtemp.dTemper),
					dtemp.dTemper = (float)(sign(MIN2(absdt,
					  0.005*phycon.te),dtemp.dTemper));
				}
			}
			else
			{
				absdt = fabs(dtemp.dTemper),
				dtemp.dTemper = (float)(sign(MIN2(absdt,0.01*
				  phycon.te),dtemp.dTemper));
			}
		}

		dtemp.dTemper *= (float)Damper;
		dtl = dtemp.dTemper;

		/* THIS IS IT !!! this is it !!! this is the first of two places where te changes. */
		phycon.te = phycon.te - dtemp.dTemper;

		/* update density - temperature variables which may have changed */
		tfidle();

		/* reevaluate gaunt factors and Boltzmann factors */
		boltgn();


		/* 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;
			CoolMHeat = 0.;
			conv.lgConvIoniz = TRUE;
		}

		if( trace.lgTrConvg>1 )
		{
			fprintf( ioQQQ, 
				"   ConvTempEdenIoniz%4li TE%9.2e dt:%8.2e%7.1f%% C:%10.2e H:%10.2e CoolMHeat:%7.1f%%\n", 
			  nThLoop, 
			  phycon.te, 
			  dtemp.dTemper, 
			  dtemp.dTemper/(phycon.te+dtemp.dTemper)*100,
			  tcool.ctot, 
			  heat.htot, 
			  CoolMHeat/heat.htot*100. );
			/* dCoolmHeatdT is the actual number used for getting dT
			 * tcool.dCooldT is just cooling */
		}

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, 
				" ConvTempEdenIoniz TE:%11.5e dT%10.3e HTOT:%11.5e CTOT:%11.5e dCoolmHeatdT:%c%11.4e C-H%11.4e lgHColionImp:%1c HDEN%10.2e\n", 
			  phycon.te, dtemp.dTemper, heat.htot, tcool.ctot, chDEType, 
			  tcool.dCoolmHeatdT, CoolMHeat, TorF(colion.lgHColionImp), 
			  phycon.hden );
		}

		/* this is test for valid thermal solution,
		 * lgIonDone = true if ionization routines got OK solution
		 * lgIonDone = false if problems
		 * TOLER can be set in READR, normally=0.02 */
		if( fabs(CoolMHeat/(heat.htot+holod.feheat)) < tolerCom.toler && 
		  conv.lgConvIoniz )
		{
			/* we have a good solution */
			/* back off a bit from the last change since it was probably too big */
			/* this is it - the second change, undoing a bit of the first */
			/* >>chng 99 nov 18, back off a bit, since temp in ism.in showed big jitter due to
			 * 1 sigma changes in T in nearly isothermal model was large */
			phycon.te = phycon.te + dtemp.dTemper/2.f;

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " ConvTempEdenIoniz returns ok.\n" );
			}

			tehigh.thist = (float)MAX2(phycon.te,tehigh.thist);
			tehigh.tlowst = (float)MIN2(phycon.te,tehigh.tlowst);

			/* test for thermal stability
			 * >>chng 97 jul 31, set flag if unstable, increment only once per zone
			 * if( lgdCNeg ) nustbl = nustbl + 1 */
			if( dcneg.lgdCNeg )
			{
				nustbl.lgUnstable = TRUE;
			}
			else
			{
				nustbl.lgUnstable = FALSE;
			}

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

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

		/* no solution, we must keep trying,
		 * this is option to increase limit to number of iterations
		 * if the solution is not oscillating
		 * nThLoop is number of times through heating-cooling-ionization loop */
		if( nThLoop == LIMDEF && !conv.lgTOscl )
		{
			/* limdef is set to 10 in a data statement 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>1 )
			{
				fprintf( ioQQQ, "   ConvTempEdenIoniz increases limtry to%3ld\n", 
				  limtry );
			}
		}
	}

	/* fall through, exceeded limit to number of solutions,
	 * no temperature convergence */
	if( fabs(CoolMHeat/(heat.htot+holod.feheat)) > tolerCom.toler)
	{
		conv.lgConvTemp = FALSE;
	}
	else
	{
		/* possible that ionization has failed, and we fell through to here
		 * with a valid thermal soln */
		conv.lgConvTemp = TRUE;
	}

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