/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */
/*lgConvEden returns true if electron density is converged */
#include "cddefines.h"
#include "dense.h"
#include "trace.h"
#include "thermal.h"
#include "bevington.h"
#include "phycon.h"
#include "converge.h"
/* limit to how many loops we will do */
/* >>chng 04 sep 27, from 25 to 35 */
#define	LOOPMAX	35

/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden 
 * returns 0 if ok, 1 if abort */
int ConvEdenIoniz(void)
{
	long int LoopEden ,
		/* this will be LOOPMAX at first, then doubled if no oscillations occur*/
		LoopLimit ,
		LoopBig;
	int lgFailLastTry;
	int kase = -1;
	static int lgSlope_oscil=FALSE;

	static int lgEden_ever_oscil = FALSE;
	static double eden_assumed_old ,
		eden_assumed_minus_true_old,
		eden_assumed_minus_true,
		slope_ls=0., slope_ls_lastzone=0.;
	static int lgEvalSlop_ls=0;
	double EdenEntry;
	static double slope=-1.;
	double 
		/* upper and lower limits to the range of the change we want to consider */
		EdenUpperLimit,
		EdenLowerLimit ,
		change_eden_rel2_tolerance ,
		PreviousChange ;

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

	/* this routine is called by ConvTempIonz, it calls ConvIonize
	 * and changes the electron density until it converges */

	/* save entry value of eden */
	EdenEntry = dense.eden;

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "  \n" );
		fprintf( ioQQQ, " ConvEdenIoniz entered \n" );
	}

	if( trace.lgTrConvg>=3 )
	{
		fprintf( ioQQQ, 
			"   ConvEdenIoniz3 called, entering eden loop using solver %sw.\n",
			conv.chSolverEden);
	}

	/* keep track of temperature solver in this zone */
	if( !conv.nTotalIoniz )
	{
		/* one time initialization of variables */
		conv.hist_temp_ntemp = -1;
		conv.hist_temp_nzone = -1;
		conv.hist_temp_limit = 0;
	}
	/* this will save the  history of density - pressure relationship
	 * for the current zone */
	if( nzone!=conv.hist_temp_nzone )
	{
		/* first time in this zone - reset counters */
		conv.hist_temp_nzone = nzone;
		/* this counter will be updated after vars are saved so will be
		 * total number of saved values */
		conv.hist_temp_ntemp = 0;
	}
	/* do we need to allocate, or reallocate, memory for the vectors */
	if( conv.hist_temp_limit==0 )
	{
		conv.hist_temp_limit = 200;
		conv.hist_temp_temp = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
		conv.hist_temp_heat = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
		conv.hist_temp_cool = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
		conv.chNotConverged = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
		conv.chConvEden = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
		conv.chConvIoniz = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
	}
	/* ran out of space - allocate more */
	if( (conv.hist_temp_ntemp+1) >= conv.hist_temp_limit )
	{
		conv.hist_temp_limit *= 3;
		if( (conv.hist_temp_temp = (double *)realloc( conv.hist_temp_temp, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double))))==NULL )
			BadMalloc();
		if( (conv.hist_temp_heat = (double *)realloc( conv.hist_temp_heat, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double))))==NULL )
			BadMalloc();
		if( (conv.hist_temp_cool = (double *)realloc( conv.hist_temp_cool, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double))))==NULL )
			BadMalloc();
	}

	/* >>chng 04 feb 11, add option to remember current density and pressure */
	conv.hist_temp_temp[conv.hist_temp_ntemp] = phycon.te;
	conv.hist_temp_heat[conv.hist_temp_ntemp] = thermal.htot;
	conv.hist_temp_cool[conv.hist_temp_ntemp] = thermal.ctot;
	++conv.hist_temp_ntemp;

	if( conv.nPres2Ioniz < 5 )
		lgEden_ever_oscil = FALSE;

	/* this flag says wether we converged the density but failed after final tweek,
	 * need to use smaller steps */
#	define PRT_FAIL_LAST_TRY	FALSE
	lgFailLastTry = FALSE;
	if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set FALSE\n");

	/* >>chng 03 mar 17, add this big loop, only go through one time if fully
	 * converged at end, but if eden jumps off in final touchup,
	 * do whole thing again */
	LoopBig = 0;
	while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) )
	{
		/* the old default solver */
		if( strcmp( conv.chSolverEden , "simple" )== 0 )
		{
			static double damp;
			long int nLoopOscil;
			LoopEden = 0;
			conv.nConvIonizFails = 0;

			/* this will be increased by 2x if, at the end, no oscillations have occurred */
			/* >>chng 04 aug 04, from 2x to 4x, to following through on eden front */
			LoopLimit = LOOPMAX;

			/* will be set true if sign of change, ever changes */
			conv.lgEdenOscl = FALSE;
			nLoopOscil = 0;
#			define PRT_EDEN_OSCIL	FALSE
			if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE2\n");
			lgFailLastTry = FALSE;
			if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set FALSE\n");

			/* will be used to look for oscillations in the electron density */
			PreviousChange = 0.;

			strcpy( conv.chConvIoniz, " NONE!!!!!" );

			/* we have not yet  */
			eden_assumed_old = 0.;
			eden_assumed_minus_true_old = 0.;

			damp = 1.;
			/* if working with zone 0 (nzone=0) reset slope to zero */
			if( !nzone )
				slope = 0.;

			/* this is electron density convergence loop */
			do 
			{
				double slopenew , slopeold ;
				/* this flag will be set false below if electron density not within tolerance 
				 * after ionization is recomputed */
				conv.lgConvEden = TRUE;

				/* converge the current level of ionization, this calls eden_sum which updates EdenTrue */
				if( ConvIoniz() || conv.lgAbort )
				{
#					ifdef DEBUG_FUN
					fputs( " <->ConvEdenIoniz()\n", debug_fp );
#					endif
					return 1;
				}
				/* this is the new error in the electron density, the difference between
				 * the assumed and true electron density - we want this to be zero */
				eden_assumed_minus_true = dense.eden - dense.EdenTrue;

				{
					static long int limEdenHistory=1000;
					static int lgMustMalloc = TRUE;
					int lgHistUpdate=FALSE;
					static double *eden_error_history=NULL, *eden_assumed_history=NULL ,
						*stdarray;
					static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1;
					if( lgMustMalloc )
					{
						lgMustMalloc = FALSE;
						lgSlope_oscil = FALSE;
						if( (eden_error_history = 
							(double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ))==NULL )
							BadMalloc();
						if( (eden_assumed_history = 
							(double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ))==NULL )
							BadMalloc();
						if( (stdarray = 
							(double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory ))==NULL )
							BadMalloc();
					}
					if( nzoneUsed!=nzone || iterUsed!=iteration )
					{
						/* first evaluation this zone */
						if( nzone==0 )
						{
							lgEvalSlop_ls = TRUE;
							slope_ls = 1.;
						}
						/* reset some vars at start of new zone */
						nEval = 0;
						nzoneUsed = nzone;
						iterUsed = iteration;
						lgSlope_oscil = FALSE;
						slope_ls_lastzone = slope_ls;
					}
					/* do not evaluate this during search phase or if slope is oscillating */
					else if( !conv.lgSearch && !lgSlope_oscil )
					{
						/* may cycle through here with exactly same edensity before tripping if  */
						int ip = MAX2(0,nEval-1);
						if( 
							/* no need to do this if already converged */
							fabs(dense.eden - dense.EdenTrue)/dense.eden > conv.EdenErrorAllowed  &&
							/* do it if never evaluated */
							(!nEval ||
							/* this test - don't want tiny corrections, within the tolerance, which
							 * can have noise since at the level code is converging, affecting the slope */
							fabs((dense.eden - dense.EdenTrue)-eden_error_history[ip]*dense.gas_phase[ipHYDROGEN] )/dense.EdenTrue > 1e-5 ||
							fabs(dense.eden-eden_assumed_history[ip]*dense.gas_phase[ipHYDROGEN])/ dense.EdenTrue > 1e-5 ))
						{
							/* use relative fraction for constant pressure case */
							eden_error_history[nEval] = (dense.eden - dense.EdenTrue)/dense.gas_phase[ipHYDROGEN];
							eden_assumed_history[nEval] = dense.eden/dense.gas_phase[ipHYDROGEN];
							stdarray[nEval] = 0.;
							++nEval;
							lgHistUpdate = TRUE;
							/*fprintf(ioQQQ,"DEBUG history\t%li\t%e\t%e\t%e\n",
								nEval,
								dense.gas_phase[ipHYDROGEN],
								eden_assumed_history[nEval],
								eden_error_history[nEval]);*/
						}
						if( nEval>=limEdenHistory )
						{
							/* used a lot of points - must create more space */
							limEdenHistory *=10;
							if( (eden_error_history = (double *)realloc( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)))==NULL )
								BadMalloc();
							if( (eden_assumed_history = (double *)realloc( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)))==NULL )
								BadMalloc();
							if( (stdarray = (double *)realloc( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double)))==NULL )
								BadMalloc();
						}
						if( nEval>3 && lgHistUpdate )
						{
							double y_intercept , stdx , stdy ,reg_coef,
								slope_save;
							slope_save = slope_ls;
							/* enough data to do least squares */
							if( linfit(eden_assumed_history, 
								eden_error_history, 
								stdarray, 
								nEval, 
								/* this says same weight for all */
								0, 
								&y_intercept, 
								&stdy, 
								&slope_ls, 
								&stdx, 
								&reg_coef) )
							{
								fprintf(ioQQQ," ConvEdenIoniz: linfit returns impossible condition.\n");
								ShowMe();
								puts( "[Stop in ConvEdenIoniz]" );
								cdEXIT(EXIT_FAILURE);
							}
							lgEvalSlop_ls = TRUE;
							/* check for slope changing sign - if does, use last zone's slope */
							if( slope_ls*slope_save<0. )
							{
								slope_ls = slope_ls_lastzone;
								lgSlope_oscil = TRUE;
							}
							/*fprintf(ioQQQ,"DEBUG LS\t%.2f\t%li\t%.3e\t%.3e\t%.3e\n",
								fnzone ,
								nEval ,
								slope_ls,
								y_intercept , 
								slope);*/
						}
					}
				}
				slopeold = slope;
				/* may not be set below, but could print it. set to zero as flag for this case */
				slopenew = 0.;
				if( fabs(eden_assumed_minus_true_old) > 0. )
				{
					if( fabs( eden_assumed_old-dense.eden ) > SMALLFLOAT )
					{
#						define	OLDFAC	0.75
						slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-dense.eden );
						/* >>chng 04 sep 15 do not update slope if changing sign */
#						if 0
						if( slope !=0. && slope*slopenew>=0.)
							slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
						else if( slope==0.)
							slope = slopenew;
#						endif
						if( slope !=0. )
							slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
						else 
							slope = slopenew;
#						undef	OLDFAC
					}
					/*else
						slope = 0.;*/
				}
				/* following block is to not let electron density change by
				* too much
				* conv.EdenErrorAllowed is allowed error
				* the default value of conv.EdenErrorAllowed is 0.01 in zerologic
				* is changed with the SET Eden Error command
				* eden_assumed_old was incoming value of eden
				* EdenTrue is correct value from eden_sum for current ionization
				* new eden will be set using these, trying to prevent oscillations */

				/* is error larger than the tolerance we are trying to beat? */
				if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
					conv.EdenErrorAllowed )
				{
					double eden_proposed;
					/* diference is assumed and true electron density is larger
					* than tolerance, we are not yet converged, default is 0.01 */
					conv.lgConvEden = FALSE;
					strcpy( conv.chConvIoniz, "Ne big chg" );

					/* SEARCH is TRUE if this is only search for first solution */
					if( conv.lgSearch )
					{
						dense.eden = pow(10.,
							(log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
					}
					else
					{
						/* >>chng 03 jul 07, add case for grains contain significant
						 * fraction of electrons */
						/* >>chng 03 nov 28, add this test on fraction of electrons from ions,
						 * this branch is to identify limit where molecules and grains have
						 * most of the electrons */
						/* >>chng 04 sep 14, change from all ions to just metals */
						/* this was patch to fix steep slope */
						if( 0 && dense.eden_from_metals > 0.5 )
						{
							static long int nzone_eval=-1;
							static int lgOscilkase10 = FALSE;
							/* this flag says be careful, if kase is this then
							 * small changes in assumed eden result in large changes
							 * in returned eden */
#							define	KASE_EDEN_NOT_ION	10
							if( nzone != nzone_eval )
							{
								nzone_eval = nzone;
								lgOscilkase10 = FALSE;
							}

							/* >>chng 03 dec 25, check whether oscillations have occurred */
							if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. && nzone_eval > 6 )
								lgOscilkase10 = TRUE;

							/* few of the electrons are due to ions, most are grains and
							* or molecules, be very careful */
							change_eden_rel2_tolerance = 0.2;
							if( lgOscilkase10 )
							{
								/* eden is oscillating, make changes small */
								change_eden_rel2_tolerance = 0.05;
							}
							kase = KASE_EDEN_NOT_ION;
						}
						/* was the sign of the change in the electron density changing,
						 * or has it ever changed? Also use if we are close to converged? */
						/* >>chng 04 aug 04, rm test on eden converged, since fooled solver when
						 * initially very close to soln, but near eden front, stopped aggressive
						 * persuit of true eden, pdr_leiden_v2 */
						else if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. )
						{
							/* remember that oscillations are happening  */
							conv.lgEdenOscl = TRUE;
							nLoopOscil = LoopEden;
							/* turn this back on 04 ep 26 */
							/* >>chng 04 sep 27, from * 0.5 to * 0.75 */
							damp = MAX(0.05,damp * 0.75 );
							if( PRT_EDEN_OSCIL ) 
								fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl TRUE, previous change %e new change %e\n",
								PreviousChange,dense.EdenTrue-dense.eden);

							/* changes in electron density are changing sign - be careful,
							* change_eden_rel2_tolerance multiplies the error tolerance on the electron density,
							* largest change allowed is related to this */
							change_eden_rel2_tolerance = 0.5;
							/* >>chng 04 aug 09, from 0.5 to 0.25 */
							change_eden_rel2_tolerance = 0.25;
							kase = 0;
						}
						else if( lgFailLastTry )
						{
							/* this is case where last evaluation of eden, after solution, 
							 * caused failure */
							change_eden_rel2_tolerance = 0.5;
							kase = 1;
						}
						/* >>chng 03 apr 22, add this branch 
						 * >>chng 04 feb 23, replace conv.lgEdenOscl with lgEden_ever_oscil */
						else if( LoopEden > 4 && !lgEden_ever_oscil &&
							fabs( (dense.EdenTrue-dense.eden)/dense.eden) > 3.*conv.EdenErrorAllowed )
						{
							/* we have gone a few times around this loop, the electron density is
							* not oscillating, and we have a long ways to go.  Open up the 
							* allowed change */
							change_eden_rel2_tolerance = 3.;
							kase = 2;
						}
						else if( conv.lgEdenOscl )
						{
							/* >>chng 04 aug 14, add this to stop oscil in primal.in */
							/* oscillatinos have occurred */
							change_eden_rel2_tolerance = 0.25;
							kase = 4;
						}
						else
						{
							/* stable so far . .. */
							/* change_eden_rel2_tolerance = 2.;*/
							/* had been 2, change to 1 stopped oscillations from developing in 
							* loc blr grid */
							change_eden_rel2_tolerance = 1.;
							kase = 3;
						}

						/* now remember this change for the next time through the loop  */
						PreviousChange = dense.EdenTrue - dense.eden ;

						/* old difference between assumed and trye */
						eden_assumed_minus_true_old = eden_assumed_minus_true;

						/* remember electron density before we update it */
						eden_assumed_old = dense.eden;

						/* is the correct electron density - is it too different? 
						* default on conv.EdenErrorAllowed is 0.01, dyanmics is 0.001 */
						EdenUpperLimit = dense.eden * (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
						EdenLowerLimit = dense.eden / (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);

#						define USE_NEW	TRUE
#						define PRT_NEW	FALSE
						if( PRT_NEW && USE_NEW )
							fprintf(ioQQQ,"DEBUG slope\t%li\t%li\t%.3f\t%.3f\t%.3f\t%.4e\t%.4e\t%.4f\t%.3f\t%.3e\n",
							nzone, 
							conv.nPres2Ioniz,
							slope , 
							slopeold , 
							slopenew ,
							dense.eden, 
							dense.EdenTrue, 
							(dense.eden-dense.EdenTrue)/SDIV(dense.eden) ,
							change_eden_rel2_tolerance ,
							EdenLowerLimit);
						if( lgEvalSlop_ls )
							slope = slope_ls;/**/
						if( USE_NEW && slope !=0. )
						{
							/* slope is d(ne_true) / d(ne_assumed) */
							/* >>chng 04 sep 26, add damp here too */
							eden_proposed = dense.eden + (dense.EdenTrue - dense.eden)/slope_ls * damp;
						}
						else
						{
							eden_proposed = dense.EdenTrue;
						}
#						if 0
						{
#include "grainvar.h"
						fprintf(ioQQQ,"DEBUG eden grn\t%e\t%e\t%e\t%e\n",
							dense.eden, eden_proposed,dense.EdenTrue, -gv.TotalEden);
						}
#						endif

						/* THIS IS IT !!! this is it !!! this is where eden changes. */
						/* get the new electron density */
						/*if( dense.EdenTrue > EdenUpperLimit )*/
						if( eden_proposed > EdenUpperLimit )
						{
							/* this branch, proposed eden too big */
							dense.eden = EdenUpperLimit;
						}
						/*else if( dense.EdenTrue < EdenLowerLimit )*/
						else if( eden_proposed < EdenLowerLimit )
						{
							/* proposed eden too small */
							dense.eden = EdenLowerLimit;
						}
						else
						{
							/* eden was within fac of the correct value, use it */
							/*dense.eden = dense.EdenTrue;*/
							dense.eden = eden_proposed;
						}

						if( trace.lgTrace && trace.lgNeBug )
						{
							fprintf( ioQQQ, 
								"     ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 
							dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue );
						}
					}
					/* we now have the proposed new electron density */
				}
				/* >>chng 00 oct 21, did not have this branch before so did not make small changes to
				* electron density (smaller than size that determined eden not converged */
				/* this branch electron density was pretty close, we just need to make a small update */
				else
				{
					/* this is test for whether will call ConvIoniz again - check how close old and correct
					* electron densities are before updating EdenTrue */
					if( fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue  > conv.EdenErrorAllowed/10. &&
						/* >>chng 03 nov 29, do not update if we are in not-ion limit of eden,
						 * this means that molecules and grains have most of the electrons */
						kase !=KASE_EDEN_NOT_ION )
					{

						/* now update eden to EdenTrue, since we are so close to is */
						/* >>chng 04 sep 20, from just update to taking into account slope */
						/*dense.eden = dense.EdenTrue;*/
						/* max is to avoid overshoots - we don't want a large correction this late */
						dense.eden = dense.eden + (dense.EdenTrue-dense.eden)/MAX2(1.,slope_ls)*damp ;

						if( trace.lgTrConvg>=3 )
						{
							fprintf( ioQQQ, 
								"   ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n",
								dense.EdenTrue , 
								dense.eden );
						}

						/* we have changed the density slightly, so now must reevaluate ionization with this new value */
						/* converge the current level of ioinization
						* but only do this if change was significant */
						/* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
						/* , this calls eden_sum which updates EdenTrue */
						if( ConvIoniz() )
						{
#							ifdef DEBUG_FUN
							fputs( " <->ConvEdenIoniz()\n", debug_fp );
#							endif
							return 1;
						}
					}
					else if( trace.lgTrConvg>=3 )
					{
						fprintf( ioQQQ, 
							"   ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n");
					}

					/* >>chng 01 mar 14, check whether electron density is no longer converged 
					* after reevaluating ionization */
					/* >>chng 04 sep 26, had div by min of eden or edentrue, use eden since
					 * EdenTrue can be negaive when not converged */
					if( 
						fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
						conv.EdenErrorAllowed )
					{
						/* difference in assumed and true electron density is larger
						* than tolerance, we are not yet converged, default is 0.01 */
						conv.lgConvEden = FALSE;
						/* >>chng 04 aug 04, do not set oscillations flag if not oscillating */
						/* make steps smaller by setting this flag
						conv.lgEdenOscl = TRUE; */
						lgFailLastTry = TRUE;
						if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set TRUE\n");
						/* this will stay true through this zone */
						lgEden_ever_oscil = TRUE;
					}
				}
				{
					/*@-redef@*/
					/* debug print statement for change in electron density */
					enum {DEBUG_LOC=FALSE};
					/*@+redef@*/
					if( DEBUG_LOC )
					{
						fprintf(ioQQQ,"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n", 
							nzone,dense.eden ,eden_assumed_old, (dense.eden-eden_assumed_old)/dense.eden, dense.EdenTrue);
					}
				}
				if( trace.lgTrace && trace.lgNeBug )
				{
					fprintf( ioQQQ, 
						"     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 
					dense.eden, eden_assumed_old, dense.eden/SDIV(eden_assumed_old),  dense.EdenTrue ,TorF( conv.lgConvEden));
				}

				if( trace.lgTrConvg>=3 )
				{
					fprintf( ioQQQ, "   ConvEdenIoniz3 loop %4ld ", 
					LoopEden );

					/* this is flag says whether or not eden/eden has converged */
					if( conv.lgConvEden )
					{
						fprintf( ioQQQ, " converged, eden rel error is %g ", 
							(dense.EdenTrue-dense.eden)/dense.EdenTrue );
					}
					else if( eden_assumed_old > SMALLFLOAT )
					{
						/*                            NB - use 8.4 fmt so chng sign not misaleign chngs*/
						fprintf( ioQQQ, " NOT Converged! corr:%8.4f prop:%9.5f ", 
						(dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 
							(dense.eden-eden_assumed_old)/eden_assumed_old );
					}

					fprintf( ioQQQ, " new ne=%.6e true=%.6e kase=%i slope=%.3e osc:%c", 
						dense.eden ,  dense.EdenTrue , kase , slope_ls,TorF(lgSlope_oscil));

					if( conv.lgEdenOscl )
					{
						fprintf( ioQQQ, " EDEN OSCIL1 damp %.4f\n", damp);
					}
					else
					{
						fprintf( ioQQQ, "\n");
					}
				}

				/* loop until converged, or we give up */
				++LoopEden;

				/* this is case where overshoots did occur, but no longer */
				/* >>chng 04 sep 23, introduce this reset */
				/* disable resetting loop unless slope is small */
				if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. )
				{
					conv.lgEdenOscl = FALSE;
					/* turn this back on 04 ep 26 */
					damp = 1.;
				}

				/* if last iteration through here and not converged, and no oscillations, 
				* and no ionization failures ,
				* then increase limit by 2x */
				if( LoopEden == (LOOPMAX-1) && !conv.lgEdenOscl && conv.nConvIonizFails==0)
					/* double the limit, but only one time, and only if no oscillations */
					/* >>chng 04 aug 04, from 2x to 4x, to follow eden changes through eden front */
					LoopLimit *= 4;

			}	while( !conv.lgConvEden && LoopEden < LoopLimit && !conv.lgAbort );
		}
		/* turned on with set eden solver new */
		else if( strcmp( conv.chSolverEden , "new" )== 0 )
		{
			/* will be used to look for bracketing in the electron density */
			double PreviousEdenError = 0.;
			double CurrentEdenError = 0.;
			double CurrentEden = 0.;
			double ProposedEden,
				FractionChangeEden = 0. ;

			/* new method of converging electron density */
			LoopEden = 0;
			conv.nConvIonizFails = 0;

			/* this will be increased by 2x if, at the end, no oscillations have occurred */
			LoopLimit = LOOPMAX;

			/* will be set true if sign of change, ever changes */
			conv.lgEdenOscl = FALSE;
			if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE1\n");

			/* this is relative change in electron density that we will permit - it will become
			* smaller each time we bracket the true electron density */
			change_eden_rel2_tolerance = 0.02;

			strcpy( conv.chConvIoniz, " NONE!!!!!" );

			/* this is ionization/electron density convergence loop
			* keep calling ionize until lgIonDone is true */
			do 
			{

				/* this flag will be set false below if electron density not within tolerance 
				* after ionization is recomputed */
				conv.lgConvEden = TRUE;

				if( trace.lgTrConvg>=3 )
					fprintf( ioQQQ, "   ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",dense.eden);

				/* converge the current level of ioinization, this calls eden_sum which updates EdenTrue */
				if( ConvIoniz() )
				{
#					ifdef DEBUG_FUN
					fputs( " <->ConvEdenIoniz()\n", debug_fp );
#					endif
					return 1;
				}

				/* the history of how we got here 
				EdenAssumed3 = EdenAssumed2,
				EdenObtained3 = EdenObtained2;
				EdenAssumed2 = EdenAssumed1,
				EdenObtained2 = EdenObtained1,*/

				/* remember the old electron density for possible printout */
				CurrentEden = dense.eden;

				/* this is the current error in the electron density */
				CurrentEdenError = dense.EdenTrue - dense.eden;

				/* conv.EdenErrorAllowed is allowed error
				* the default value of conv.EdenErrorAllowed is 0.01 in zerologic
				* is changed with the SET Eden Error command
				* eden_assumed_old was incoming value of eden
				* EdenTrue is correct value from eden_sum for current ionization
				* new eden will be set using these, trying to prevent oscillations */

				/* is error larger than the tolerance we are trying to beat?
				* first check is whether error is very small after very first check, since
				* ionization soln may not have settled down yet */
				if( 
					(LoopEden==0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed/2.) ||
					(LoopEden>0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed 
					&& FractionChangeEden < conv.EdenErrorAllowed / 2.) )
					break;

				/* diference is assumed and true electron density is larger
				* than tolerance, we are not yet converged, default is 0.01 */
				conv.lgConvEden = FALSE;
				strcpy( conv.chConvIoniz, "Ne big chg" );

				/* SEARCH is TRUE if this is only search for first solution */
				if( conv.lgSearch )
				{
					dense.eden = pow(10.,
						(log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
				}
				else
				{

					/* was the sign of the change in the electron density changing,
					* if so then we have bracked the target */
					if( PreviousEdenError*CurrentEdenError < 0. )
					{
						/* we have bracketed the correct electron density, 
						 * make changes smaller, also solve linear equation for zero error */
						slope = (PreviousEdenError - CurrentEdenError ) /
							( eden_assumed_old - CurrentEden );

						ProposedEden = eden_assumed_old - PreviousEdenError / slope;
						/* as sanity check, this must be between the two limits we examined,
						* since zero was between them */
						ASSERT( ProposedEden >= MIN2( eden_assumed_old , CurrentEden ) );
						ASSERT( ProposedEden <= MAX2( eden_assumed_old , CurrentEden ) );

						change_eden_rel2_tolerance *= 0.25;
						if( trace.lgTrConvg>=3 )
							fprintf( ioQQQ, 
							"   ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n",
							change_eden_rel2_tolerance,
							(dense.EdenTrue-dense.eden)/dense.EdenTrue, ProposedEden);
					}
					else
					{
						/* did not bracket the target, keep moving in its direction */

						/* the correct electron density - is it too different? */
						EdenUpperLimit = dense.eden * (1. + change_eden_rel2_tolerance);
						EdenLowerLimit = dense.eden / (1. + change_eden_rel2_tolerance);

						/* get the new electron density */
						if( dense.EdenTrue > EdenUpperLimit )
						{
							/* this branch, proposed eden too big */
							ProposedEden = EdenUpperLimit;
						}
						else if( dense.EdenTrue < EdenLowerLimit )
						{
							/* proposed eden too small */
							ProposedEden = EdenLowerLimit;
						}
						else
						{
							/* eden was within fac of the correct value, use it */
							ProposedEden = dense.EdenTrue;
						}
						if( trace.lgTrConvg>=3 )
							fprintf( ioQQQ, 
							"   ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n",
							change_eden_rel2_tolerance,
							(dense.EdenTrue-dense.eden)/dense.EdenTrue ,
							ProposedEden,
							dense.EdenTrue,slope);
					}

					/* now remember this change for the next time through the loop  */
					PreviousEdenError = CurrentEdenError ;
					eden_assumed_old = CurrentEden;
					FractionChangeEden = fabs( dense.eden - ProposedEden ) / dense.EdenTrue;
					dense.eden = ProposedEden;

					if( trace.lgTrace && trace.lgNeBug )
					{
						fprintf( ioQQQ, 
							"     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 
						dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue );
					}
				}
				/* loop until we give up  --  normal exit is break in center of loop */
				++LoopEden;
			}	while( LoopEden < LoopLimit &&  !conv.lgAbort );

			/* we now have the proposed new electron density */
			/* we have changed the density slightly, so now must reevaluate ionization with this new value */
			/* converge the current level of ioinization
			* but only do this if change was significant */
			if( trace.lgTrConvg>=3 )
				fprintf( ioQQQ, 
				"   ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n",
				(dense.EdenTrue - dense.eden)/dense.EdenTrue,dense.EdenTrue);
			dense.eden = dense.EdenTrue;
			/* , this calls eden_sum which updates EdenTrue */
			/* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
			if( ConvIoniz() )
			{
				if( trace.lgTrConvg>=3 )
					fprintf( ioQQQ, 
					"   ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n",
					(dense.EdenTrue - dense.eden)/dense.EdenTrue);
			}

			/* >>chng 01 mar 14, check whether electron density is no longer converged 
			* after reevaluating ionization */
			if( 
				fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue > 
				conv.EdenErrorAllowed )
			{
				/* diference is assumed and true electron density is larger
				* than tolerance, we are not yet converged, default is 0.01 */
				conv.lgConvEden = FALSE;
				if( trace.lgTrConvg>=3 )
					fprintf( ioQQQ, 
					"   ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n",
					(dense.eden-dense.EdenTrue)/dense.EdenTrue );
			}
			else
			{
				conv.lgConvEden = TRUE;
			}

			if( trace.lgTrace && trace.lgNeBug )
			{
				fprintf( ioQQQ, 
					"     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 
				dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue ,TorF( conv.lgConvEden));
				/* update density and temperature related variables
				tfidle(FALSE); */
			}

			if( trace.lgTrConvg>=3 )
			{
				fprintf( ioQQQ, "   ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e", 
				LoopEden, dense.eden ,  dense.EdenTrue );

				/* this is flag says whether or not eden/eden has converged */
				if( conv.lgConvEden )
				{
					fprintf( ioQQQ, " converged, eden rel error is %g ", 
						(dense.EdenTrue-dense.eden)/dense.EdenTrue );
				}
				else
				{
					fprintf( ioQQQ, " NOCONV corr:%7.3f prop:%7.3f ", 
					(dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 
						(dense.eden-eden_assumed_old)/eden_assumed_old );
				}
				if( conv.lgEdenOscl )
					fprintf( ioQQQ, " EDEN OSCIL2 \n");
			}

		}
		else
		{
			fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s \n" , conv.chSolverEden );
			ShowMe();
		}
		++LoopBig;
	}

	if( trace.lgTrConvg>=3 )
	{
		fprintf( ioQQQ, "   ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 
			TorF(conv.lgConvEden ),
			EdenEntry ,
			dense.eden ,
			(EdenEntry-dense.eden)/SDIV( dense.eden ) );
	}
	else if( trace.lgTrace )
	{
		fprintf( ioQQQ, " ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 
			fnzone,
			TorF(conv.lgConvEden ),
			EdenEntry ,
			dense.eden ,
			(EdenEntry-dense.eden)/SDIV( dense.eden ) );
	}

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

	return 0 ;
}

/* returns true if electron density is converged */
int lgConvEden(void)
{
	int lgRet;

	/* >>chng 04 sep 26, div by eden, not min of eden and edentrue since latter now poss neg */
	if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
		conv.EdenErrorAllowed )
	{
		lgRet = FALSE;
		conv.lgConvEden = FALSE;
	}
	else
	{
		lgRet = TRUE;
		conv.lgConvEden = TRUE;
	}
	return lgRet;
}
#undef	LOOPMAX
