/* 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 */
/*ConvInitSolution drive search for initial temperature, for illuminated face */
#include "cddefines.h"
#include "trace.h"
#include "stopcalc.h"
#include "struc.h"
#include "rfield.h"
#include "dense.h"
#include "heavy.h"
#include "wind.h"
#include "geometry.h"
#include "thermal.h"
#include "radius.h"
#include "phycon.h"
#include "pressure.h"
#include "ionbal.h"
#include "conv.h"
/* XMULTP is scale factor to multiply temperature */
#define	XMULTP	3.0

int ConvInitSolution(void)
{
	long int i, 
	  ionlup, 
	  isear, 
	  j, 
	  looplimit,
	  nelem ,
	  nflux_old,
	  nelem_reflux ,
	  ion_reflux;

	double cnew, 
	  cold, 
	  dt, 
	  fac2, 
	  factor, 
	  hnew, 
	  hold, 
	  Cool_min_Heat, 
	  t1, 
	  t2, 
	  tinc, 
	  tmax, 
	  tmin, 
	  tnew, 
	  told;

	static double TempPrevIteration = -DBL_MAX;
	double derrdt;

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

	tnew = DBL_MAX;
	cnew = DBL_MAX;
	hnew = DBL_MAX;
		
	/* this counts number of times ConvBase is called by PressureChange, in current zone */
	conv.nPres2Ioniz = 0;
	/* this counts how many times ConvBase is called in this iteration after startr,
	 * and is flag used by ConvBase to understand it is being called the first time*/
	conv.nTotalIoniz = 0;
	/* ots rates not oscillating */
	conv.lgOscilOTS = FALSE;

	lgAbort = FALSE;
	dense.lgEdenBad = FALSE;
	dense.nzEdenBad = 0;
	/* these are variables to remember the biggest error in the
	 * electron density, and the zone in qhich it occurred */
	conv.BigEdenError = 0.;
	conv.AverEdenError = 0.;
	conv.BigHeatCoolError = 0.;
	conv.AverHeatCoolError = 0.;
	conv.BigPressError = 0.;
	conv.AverPressError = 0.;

	/* these are equal if set dr was used, and we know the first dr */
	/*lint -e777 test float equality */
	if( radius.sdrmin == radius.sdrmax )
	{
		radius.drad = MIN2( radius.sdrmax  ,  radius.drad );
		radius.drad_x_fillfac = radius.drad * geometry.FillFac;
	}
	/*lint +e777 test float equality */

	ASSERT( dense.xIonDense[ipHYDROGEN][0] >0 && dense.xIonDense[ipHYDROGEN][1]>= 0.);

	if( trace.lgTrConvg )
	{
		fprintf( ioQQQ, "  \n" );
		fprintf( ioQQQ, " ConvInitSolution entered \n" );
	}

	/********************************************************************
	 *
	 * force an initial temperature, also constant temperature 
	 *
	 *********************************************************************/
	if( thermal.ConstTemp != 0. )
	{
		/* force te to certain value */
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ConvInitSolution called, forcing TE to%10.2e\n", 
			  thermal.ConstTemp );
		}
		phycon.te = thermal.ConstTemp;
		tfidle(FALSE);

		/* evaluate current pressure for this constant temperature model */
		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* now remember these inital values */
		pressure.PresRamInit = pressure.PresRamCurr;
		pressure.PresTotlInit = pressure.PresTotlCurr;

		/* say this is search for first solution, let Ne change a lot */
		conv.lgSearch = TRUE;

		ionlup = 0;
		conv.lgConvIoniz = FALSE;

		/* >>chng 96 oct 11, loop at least 3x to make sure ionization converged */
		/* >>chng 02 jul 11, slight rewrite of logic, keep iterating as long as
		 * 1) have not aborted,
		 * 2) have not gone over loop three times, and
		 * 3) ionization not converged with ionlup <= 7 
		 */
		while( (ionlup<3) || (!lgAbort && (ionlup <= 3) && !conv.lgConvIoniz && (ionlup <= 7) ) )
		{
			/* >>chng 02 jun 10, this calls tfidle, among other things */
			/* this sets values of pressure.PresTotlCurr */
			PresTotCurrent();
			if( ConvEdenIoniz() )
			{
#				ifdef DEBUG_FUN
				fputs( " <->ConvInitSolution()\n", debug_fp );
#				endif
				/* this is an error return, calculation will immediately stop */
				return(1);
			}
			if( ionlup<3 )
			{
				/* >>chng 03 may 18, set first zone thickness here */
				/* set thickness of first zone, this affects the pumping rates due
				 * to correction for attenuation across zone, so must be known */
				radius_first();
			}
			/* only used to keep forcing new optical depths below */
			++ionlup;
		}

		if( trace.lgTrace || trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ========================================================================================================================\n");
			fprintf( ioQQQ, " ConvInitSolution: search set false 1 Te=%.3e================================================================================\n" , phycon.te);
			fprintf( ioQQQ, " ========================================================================================================================\n");
		}
		conv.lgSearch = FALSE;
	}

	/********************************************************************
	 *
	 * this is second or higher iteration, reestablish original temperature
	 *
	 *********************************************************************/
	else if( iteration != 1 )
	{
		/* this is second or higher iteration on multi-iteration model */
		if( trace.lgTrace || trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ConvInitSolution called, ITER=%2ld resetting Te to %10.2e\n", 
			  iteration, TempPrevIteration );
		}

		if( trace.lgTrace || trace.lgTrConvg )
		{
			fprintf( ioQQQ, " search set true\n" );
		}

		/* search phase must be turned on so that variables such as the ors rates,
		 * secondary ionization, and auger yields, can converge more quickly to 
		 * proper values */
		conv.lgSearch = TRUE;

		/* this is the temperature and pressure from the previous iteration */
		phycon.te = TempPrevIteration;
		tfidle(FALSE);

		/* the inital pressure should now be valid */
		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* now remember inital values */
		pressure.PresRamInit = pressure.PresRamCurr;
		pressure.PresTotlInit = pressure.PresTotlCurr;

		/* >>chng 03 aug 21, rm zero of nTotalIoniz */
		/* this counts how many times ConvBase is called in this iteration after startr,
		 * and is flag used by ConvBase to understand it is being called the first time
		conv.nTotalIoniz = 0;*/

		/* now get final temperature */
		if( ConvTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

		if( trace.lgTrace || trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ========================================================================================================================\n");
			fprintf( ioQQQ, " ConvInitSolution: search set false 2 Te=%.3e========================================================================================\n" , phycon.te);
			fprintf( ioQQQ, " ========================================================================================================================\n");
		}
		conv.lgSearch = FALSE;

		/* now get final temperature */
		if( ConvTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

		/* the inital pressure should now be valid */
		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* now remember these inital values */
		pressure.PresRamInit = pressure.PresRamCurr;
		pressure.PresTotlInit = pressure.PresTotlCurr;

		/* remember the current temperature for next time */
		TempPrevIteration = phycon.te;

	}

	else
	{
		/********************************************************************
		 *
		 * do first te from scratch 
		 *
		 *********************************************************************/

		/* say that we are in search phase */
		conv.lgSearch = TRUE;

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

		/* approach equilibrium from either high or low TE */
		if( thermal.lgTeHigh )
		{
			/* approach from high TE */
			phycon.te = 1e6;
			tfidle(FALSE);
			ionbal.lgNoDec = FALSE;
			factor = 1./XMULTP;
		}
		else
		{
			/* approach from low TE */
			ionbal.lgNoDec = TRUE;
			phycon.te = 4000.;
			tfidle(FALSE);
			factor = XMULTP;
		}

		/* find first temperature */
		thermal.lgColNeg = FALSE;

		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* now remember these inital pressures */
		pressure.PresTotlInit = pressure.PresTotlCurr;
		pressure.PresRamInit = pressure.PresRamCurr;

		thermal.htot = 1.;
		thermal.ctot = 1.;

		/* call cooling, heating, opacity, loop to convergence
		 * this is very first call to it, by default is at 4000K */
		if( trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ConvInitSolution calling ConvEdenIoniz1 with te=%10.3e\n", 
			  phycon.te );
		}

		for( ionlup=0; ionlup<2; ++ionlup )
		{
			/*tfidle(FALSE);*/
			if( ConvEdenIoniz() )
			{
#				ifdef DEBUG_FUN
				fputs( " <->ConvInitSolution()\n", debug_fp );
#				endif
				/* this is an error return, calculation will immediately stop */
				return(1);
			}
			/* >>chng 03 may 18, set first zone thickness here */
			/* set thickness of first zone, this affects the pumping rates due
			 * to correction for attenuation across zone, so must be known */
			radius_first();
		}

		if( lgAbort )
		{
			/* we hit an abort */
			fprintf( ioQQQ, " Search for inital conditions aborted - lgAbort set true.\n" );
			ShowMe();

#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

		/* >>chng 97 feb 18, define following three */
		told = phycon.te;
		cold = thermal.ctot;
		hold = thermal.htot;

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " RETURN TO ConvInitSolution 1, HTOT:%10.3e  CTOT:%10.3e  Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", 
			  thermal.htot, thermal.ctot, phycon.te, dense.eden );
		}

		/* Cool_min_Heat is cooling minus heating, in most cases negative means T too low */
		Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);

		/* if low te approach then check if already past soln (cooling > heating) */
		if( !thermal.lgTeHigh && Cool_min_Heat >= 0. )
		{
			/* this branck, low Te approach to thermal soln and we started out too hot,
			 * so will gradually lower the temperature */
			if( trace.lgTrConvg )
			{
				fprintf( ioQQQ, " ConvInitSolution entering T<ConvInitSolution loop, te htot ctot=%10.3e%10.3e%10.3e Cool_min_Heat=%9.1e\n", 
				  phycon.te, thermal.htot, thermal.ctot, Cool_min_Heat );
			}
			/* check range 10-4000K, TeLowest usually 2.8K, set with SET LOWEST command
			 * >>chng 96 may 30, logic had gone straight to 3K, which blew up dense clouds
			 * now go down more gently
			 * >>chng 97 mar 03, variable factor, smaller when close to 3K */
			fac2 = 0.8;
			while( phycon.te*fac2 > StopCalc.TeLowest && Cool_min_Heat >= 0. && !lgAbort )
			{
				/* >>chng 96 jun 14, to below, so that we can get soln near 3K
				 * do while( te.gt.TeLowest .and. Cool_min_Heat.ge.0. )
				 * >>chng 97 feb 18, define following three */
				told = phycon.te;
				cold = thermal.ctot;
				hold = thermal.htot;

				/* >>chng 97 mar 3, make fac2 small when close to lower bound */
				if( phycon.te < 5. )
				{
					fac2 = 0.97;
				}
				else if( phycon.te < 20. )
				{
					fac2 = 0.9;
				}
				else
				{
					fac2 = 0.8;
				}

				/* now change to new temp */
				phycon.te = phycon.te * fac2;
				tfidle(TRUE);

				/* this sets values of pressure.PresTotlCurr */
				PresTotCurrent();

				/* now remember inital pressures */
				pressure.PresTotlInit = pressure.PresTotlCurr;
				pressure.PresRamInit = pressure.PresRamCurr;

				/* call cooling, heating, opacity, loop to convergence */
				if( trace.lgTrConvg )
				{
					fprintf( ioQQQ, "\n calling ConvEdenIoniz2, te=%9.2e\n", 
					  phycon.te );
				}

				/*tfidle(FALSE);*/
				if( ConvEdenIoniz() )
				{
#					ifdef DEBUG_FUN
					fputs( " <->ConvInitSolution()\n", debug_fp );
#					endif
					/* this is an error return, calculation will immediately stop */
					return(1);
				}
				tnew = phycon.te;
				cnew = thermal.ctot;
				hnew = thermal.htot;

				if( trace.lgTrConvg )
				{
					fprintf( ioQQQ, " ConvEdenIoniz2 returns, te=%9.2e htot%9.2e ctot%9.2e\n", 
					  phycon.te, thermal.htot, thermal.ctot );
				}
				if( trace.lgTrace || trace.lgTrConvg )
				{
					fprintf( ioQQQ, " RETURN TO ConvInitSolution 2  HTOT:%10.3e  CTOT:%10.3e  Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", 
					  thermal.htot, thermal.ctot, phycon.te, dense.eden );
				}

				Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);
			}

			if( lgAbort )
			{
				/* we hit an abort */
				fprintf( ioQQQ, " Search for inital conditions aborted - lgAbort set true.\n" );
				ShowMe();
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvInitSolution()\n", debug_fp );
#				endif
				/* this is an error return, calculation will immediately stop */
				return(1);
			}
			else if( thermal.ctot > thermal.htot )
			{
				/* even at lowest temp, too much cooling; no solution possible */
				fprintf( ioQQQ, " Equilibrium temperture below T=%8.2e.  If this is OK, then reset with LOWEST command.  This is ConvInitSolution.  TeLowest was%9.2e\n", 
				  phycon.te, StopCalc.TeLowest );
				
#				ifdef DEBUG_FUN
				fputs( " <->ConvInitSolution()\n", debug_fp );
#				endif
				/* this is an error return, calculation will immediately stop */
				return(1);
			}

			/* >>chng 97 mar 3, from *0.9 to / fac2 */
			phycon.te = phycon.te / fac2;
			tfidle(TRUE);
			factor = 1./(0.5*(1. + fac2));
			Cool_min_Heat = fabs(cold-hold)/(cold - hold);
			if( trace.lgTrConvg )
			{
				fprintf( ioQQQ, " <4000K loop returns, te=%9.2e Cool_min_Heat%9.2e factor%9.2e\n", 
				  phycon.te, Cool_min_Heat, factor );
			}
		}

		tmin = 1e9;
		tmax = 0.;

		/* coming up to here te is either 4,000 (usually) or 10^6 */
		t1 = StopCalc.TeLowest/1.001;

		/* >>chng 96 dec 28, added div by factor to prevent checking on temp > 1e10 */
		t2 = StopCalc.TeHighest*1.001/factor;
		if( trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ConvInitSolution entering second loop te=%9.2e Cool_min_Heat%9.2e factor%9.2e t1, t2=%10.2e%10.2e\n", 
			  phycon.te, Cool_min_Heat, factor, t1, t2 );
		}
		/* loop until sign of heating - cooling mismatch changes, and we are in Te bounds */
		while( (thermal.ctot - thermal.htot)*Cool_min_Heat > 0. && phycon.te > t1 && phycon.te < t2 )
		{
			/* >>chng 97 feb 18, define following three */
			told = phycon.te;
			cold = thermal.ctot;
			hold = thermal.htot;
			phycon.te = phycon.te * factor;
			tfidle(TRUE);
			/*fprintf(ioQQQ,"DEBUG loop 2 te %.2e %.2e\n", phycon.te , factor);*/
			/* >>chng 05 jan 15, we are increasing temperature above default first guess so must also
			 * reset highest ionization stages - this fixed major logical error when sim is
			 * totally mechanically heated to coronal temperatures 
			 * problem id'd by Ronnie Hoogerwerf */
			if( factor > 1. )
			{
				for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
				{
					if( dense.lgElmtOn[nelem] )
					{
						dense.IonHigh[nelem] = nelem + 1;
					}
				}
			}

			/* this sets values of pressure.PresTotlCurr */
			PresTotCurrent();

			/* now remember inital pressures */
			pressure.PresTotlInit = pressure.PresTotlCurr;
			pressure.PresRamInit = pressure.PresRamCurr;

			tmin = MIN2(tmin,phycon.te);
			tmax = MAX2(tmax,phycon.te);
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "ConvInitSolution srh" );
			isear = 1;
			while( isear <= 6 && !conv.lgConvIoniz && !lgAbort )
			{
				/* call cooling, heating, opacity, loop to convergence */
				if( trace.lgTrConvg )
				{
					fprintf( ioQQQ, "\n calling ConvEdenIoniz3, te=%9.2e\n", 
					  phycon.te );
				}

				/* >>chng 02 jun 10, this calls tfidle, among other things */
				/* this sets values of pressure.PresTotlCurr */
				PresTotCurrent();
				if( ConvEdenIoniz() )
				{
#					ifdef DEBUG_FUN
					fputs( " <->ConvInitSolution()\n", debug_fp );
#					endif
					/* this is an error return, calculation will immediately stop */
					return(1);
				}
				if( trace.lgTrConvg )
				{
					fprintf( ioQQQ, " ConvEdenIoniz3 returns, te=%9.2e htot%9.2e ctot%9.2e\n", 
					  phycon.te, thermal.htot, thermal.ctot );
				}
				if( trace.lgTrace )
				{
					fprintf( ioQQQ, " RETURN TO ConvInitSolution 3  HTOT:%10.3e  ctot:%10.3e  Te:%11.4e eden%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", 
					  thermal.htot, thermal.ctot, phycon.te, dense.eden );
				}
				isear += 1;
			}
			tnew = phycon.te;
			cnew = thermal.ctot;
			hnew = thermal.htot;
		}

		/* fall through loop, te too high or low */
		if( phycon.te <= t1 || phycon.te >= t2 )
		{
			if( phycon.te <= t1 )
			{
				fprintf( ioQQQ, " ConvInitSolution finds TE below TE = %.3eK, HTOT=%12.4e CTOT=%12.4e\n", 
				  t1, thermal.htot, thermal.ctot );
			}
			else if( phycon.te >= t2 )
			{
				fprintf( ioQQQ, " ConvInitSolution finds TE above TE = %.3eK, HTOT=%12.4e last CTOT=%12.4e\n", 
				  t2, thermal.htot, thermal.ctot );
				fprintf( ioQQQ, 
					" The electron temperature is too high; the electrons are relativistic and little optical/UV radiation will result.\n");
				fprintf( ioQQQ, 
					" This is often due to the continuum shape being wrong - unphysically strong radio or gamma ray emission will do this.\n");
			}
			else
			{
				fprintf( ioQQQ, " ConvInitSolution finds insanity \n" );
				ShowMe();
			}
			
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error condition, calculation will immediately stop */
			return(1);
		}

		/* now search inside (also either side of) TE range where sign changed
		 * first turn on ionization stage trimming */
		ionbal.lgNoDec = FALSE;
		thermal.lgColNeg = TRUE;

		/* change sense of temp search so that we can "weave" into soln */
		tinc = 1./pow(factor,0.2);

		for( j=0; j < 2; j++ )
		{
			if( (thermal.ctot - thermal.htot > 1e-30) && thermal.ctot > 0. && !lgAbort )
			{
				Cool_min_Heat = fabs(thermal.ctot-thermal.htot)/(thermal.ctot - thermal.htot);
				if( trace.lgTrConvg )
				{
					fprintf( ioQQQ, "\n\n >>>entering second loop, te reset to%9.2e\n", 
					  phycon.te );
				}

				/* this sets values of pressure.PresTotlCurr */
				PresTotCurrent();

				/* now remember initial pressures */
				pressure.PresTotlInit = pressure.PresTotlCurr;
				pressure.PresRamInit = pressure.PresRamCurr;

				/* call cooling, heating, opacity, loop to convergence */
				isear = 1;
				conv.lgConvIoniz = FALSE;

				if( trace.lgTrace )
				{
					fprintf( ioQQQ, " RETURN TO ConvInitSolution 4  HTOT:%10.3e  CTOT:%10.3e  Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", 
					  thermal.htot, thermal.ctot, phycon.te, dense.eden );
				}
				i = 0;

				/* >>chng 97 jul 29, to recover fumbled initial ionization */
				looplimit = 12;
				told = phycon.te;
				cold = thermal.ctot;
				hold = thermal.htot;

				/* this loop is not necessariy executed (ctot and htot very close) */
				while( i < looplimit && (thermal.ctot - thermal.htot)*Cool_min_Heat > 0. && 
					!lgAbort )
				{
					i += 1;
					told = phycon.te;
					cold = thermal.ctot;
					hold = thermal.htot;
					phycon.te = phycon.te *tinc;
					tfidle(TRUE);

					/* this sets values of pressure.PresTotlCurr */
					PresTotCurrent();

					/* now remember initial pressures */
					pressure.PresTotlInit = pressure.PresTotlCurr;
					pressure.PresRamInit = pressure.PresRamCurr;

					/*tfidle(FALSE);*/
					if( ConvEdenIoniz() )
					{
#						ifdef DEBUG_FUN
						fputs( " <->ConvInitSolution()\n", debug_fp );
#						endif
						/* this is an error return, calculation will immediately stop */
						return(1);
					}

					if( trace.lgTrConvg )
					{
						fprintf( ioQQQ, " ConvEdenIoniz returned, te=%9.2e htot%9.2e ctot%9.2e\n", 
						  phycon.te, thermal.htot, thermal.ctot );
					}

					tnew = phycon.te;
					cnew = thermal.ctot;
					hnew = thermal.htot;

					if( trace.lgTrace )
					{
						fprintf( ioQQQ, " RETURN TO ConvInitSolution 4  HTOT:%10.3e  CTOT:%10.3e  Te:%11.4e EDEN%11.4e<<<<<<<<<<<<<<<<<<<<<<<<<\n", 
						  thermal.htot, thermal.ctot, phycon.te, dense.eden );
					}
				}

				/* >>chng 04 jan 25, rm this block, no longer necessary to prevent fully molecular limit 
				 * as Nick Abel pointed out */
#				if 0
				if( (i >= looplimit && (thermal.ctot - thermal.htot)*Cool_min_Heat > 0.) && 
					fabs(thermal.ctot/thermal.htot-1.) > 0.05 )
				{
					/* >>chng 96 dec 17, way out for dense, cold, molecular gas */
					if( phycon.te < 1e3 && co.hevmol[ipCO]/dense.gas_phase[ipCARBON] > 0.5 )
					{
						fprintf( ioQQQ, " cold molecular trap, te=%10.2e F(CO)=%10.2e\n", 
						  phycon.te, co.hevmol[ipCO]/dense.gas_phase[ipCARBON] );
						fprintf( ioQQQ, " Conditions in this cloud are not appropriate for this code.\n" );
						
#						ifdef DEBUG_FUN
						fputs( " <->ConvInitSolution()\n", debug_fp );
#						endif
						/* this is an error return, calculation will immediately stop */
						return(1);
					}
					else
					{
						fprintf( ioQQQ, " Insanity occurred in ConvInitSolution, i=0\n" );
						ShowMe();
						puts( "[Stop in ConvInitSolution]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
#				endif
				/* weave back in again */
				tinc = 1./pow(tinc,0.2);
			}
		}

		if( trace.lgTrace || trace.lgTrConvg )
		{
			fprintf( ioQQQ, " ========================================================================================================================\n");
			fprintf( ioQQQ, " ConvInitSolution: search set false 3 Te=%.3e========================================================================================\n" , phycon.te);
			fprintf( ioQQQ, " ========================================================================================================================\n");
		}
		conv.lgSearch = FALSE;

		/* now must be close enough to call regular ionization-cooling
		 * first do linear interpolation to guess temperature */
		/*lint -e777 float test equality */
		if( tnew != told && tnew > 0. && told > 0. )
		/*lint +e777 float test equality */
		{
			derrdt = ((cnew-hnew)-(cold-hold))/(tnew-told);
			dt = -(cold-hold)/derrdt;
			phycon.te = (told + dt);
			phycon.te = MAX2(phycon.te,MIN2(told,tnew));
			phycon.te = MIN2(phycon.te,MAX2(told,tnew));
			/* following must be defined for ram pressure to be computed in PresTotCurrent */
			struc.DenMass[0] = dense.xMassDensity;

			/* this sets values of pressure.PresTotlCurr */
			PresTotCurrent();

			/* now remember current pressures */
			pressure.PresTotlInit = pressure.PresTotlCurr;
			pressure.PresRamInit = pressure.PresRamCurr;
		}

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ConvInitSolution calls IONTE after interpolating a temperature of%10.3e Told=%10.2e Tnew=%10.2e dC-Hold%10.2e dC-Hnew%10.2e\n", 
			  phycon.te, told, tnew, cold - hold, cnew - hnew );
		}

		if( trace.lgTrConvg )
		{
			fprintf( ioQQQ, "\n\n ConvInitSolution calling ConvTempEdenIoniz with interpolated te=%10.3e\n", 
			  phycon.te );
		}

		if( ConvTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}
		TempPrevIteration = phycon.te;

		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " ConvInitSolution return, TE:%10.2e==================\n", 
			  phycon.te );
		}

	}

	/* this counts number of times ConvBase is called by PressureChange, in current zone
	 * these are reset here, so that we count from first zone not search */
	conv.nPres2Ioniz = 0;

	dense.lgEdenBad = FALSE;
	dense.nzEdenBad = 0;

	/* these are variables to remember the biggest error in the
	 * electron density, and the zone in qhich it occurred */
	conv.BigEdenError = 0.;
	conv.AverEdenError = 0.;
	conv.BigHeatCoolError = 0.;
	conv.AverHeatCoolError = 0.;
	conv.BigPressError = 0.;
	conv.AverPressError = 0.;

	/* now remember some things we may need even in first zone, these
	 * are normally set towards end of zone calc in RT_tau_inc */
	struc.testr[0] = (float)phycon.te;
	/* >>chng 02 May 2001 rjrw: add hden for dilution */
	struc.hden[0] = dense.gas_phase[ipHYDROGEN];
	/* pden is the total number of particles per unit vol */
	struc.DenParticles[0] = dense.pden;
	struc.heatstr[0] = thermal.htot;
	struc.coolstr[0] = thermal.ctot;
	struc.volstr[0] = (float)radius.dVeff;
	struc.drad_x_fillfac[0] = (float)radius.drad_x_fillfac;
	struc.histr[0] = dense.xIonDense[ipHYDROGEN][0];
	struc.hiistr[0] = dense.xIonDense[ipHYDROGEN][1];
	struc.ednstr[0] = (float)dense.eden;
	struc.o3str[0] = dense.xIonDense[ipOXYGEN][2];
	struc.DenMass[0] = dense.xMassDensity;
	struc.drad[0] = (float)radius.drad;

	/* check that nflux extends above IP of highest ionization species present.
	 * for collisional case it is possible for species to exist that are higher
	 * IP than the limit to the continuum.  Need continuum to encompass all
	 * possible emission - to account for diffuse emission
	 * NB 
	 * on second iteration of multi-iteration model this may result in rfield.nflux increasing
	 * which can change the solution */
	nflux_old = rfield.nflux;
	nelem_reflux = -1;
	ion_reflux = -1;
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		/* do not include hydrogenic species in following */
		for( i=0; i < nelem; i++ )
		{
			if( dense.xIonDense[nelem][i+1] > 0. )
			{
				if( Heavy.ipHeavy[nelem][i] > rfield.nflux )
				{
					rfield.nflux = Heavy.ipHeavy[nelem][i] ;
					nelem_reflux = nelem;
					ion_reflux = i;
				}
			}
		}
	}

	/* was the upper limit to the opacities updated? */
	if( nflux_old != rfield.nflux )
	{
		/* zero out parts of rfield arrays that were previously undefined */
		rfield_opac_zero( nflux_old-1 , rfield.nflux );

		/* if continuum reset up, we need to define gaunt factors through high end */
		/*tfidle(FALSE);*/
		/* this calls tfidle, among other things */
		/* this sets values of pressure.PresTotlCurr */
		PresTotCurrent();

		/* redo ionization and update opacities */
		if( ConvBase(1) )
		{
			/* this is catastropic failure */
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			return( 1 );
		}

		/* need to update continuum opacities */
		if( trace.lgTrace )
		{
			fprintf(ioQQQ," nflux updated from %li to %li, anu from %g to %g \n",
				nflux_old , rfield.nflux ,
				rfield.anu[nflux_old] , rfield.anu[rfield.nflux] );
			fprintf(ioQQQ," caused by element %li ion %li \n",
			  nelem_reflux ,ion_reflux );
		}
	}

	/* wind with negative velocity is flow from pdr - change density slightly
	 * and call pressure solver to see if it returns to original density */
	if( (strcmp(dense.chDenseLaw,"WIND") == 0) && wind.windv < 0. )
	{
		/* >>chng 04 apr 23, change pressure and let solver come back to correct
		 * pressure.  This trick makes sure we are correctly either super or sub sonic 
		 * since solver will jump from one branch to the other if necessary */
#		define PCHNG	0.98
		/* this ignores return condition, assume must be ok if got this far */
		if( ConvPresTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

		pressure.PresTotlInit *= PCHNG;
		if( ConvPresTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

		pressure.PresTotlInit /= PCHNG;
		if( ConvPresTempEdenIoniz() )
		{
#			ifdef DEBUG_FUN
			fputs( " <->ConvInitSolution()\n", debug_fp );
#			endif
			/* this is an error return, calculation will immediately stop */
			return(1);
		}

#		undef PCHNG
	}

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

	/* this is the only valid return and lgAbort should be false*/
	return( lgAbort );
}

