/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */
#include <string.h>
#include "cddefines.h"
#include "tfidle.h"
#include "phycon.h"
#include "trace.h"
#include "elecden.h"
#include "converge.h"
/* limit to how many loops we will do */
#define	LOOPMAX	20

void ConvEdenIoniz(void)
{
	long int LoopEden ,
		/* this will be LOOPMAX at first, then doubled if no oscillations occur*/
		LoopLimit ,
		/* the number of failures passed up here from convIoniz */
		nConvIonizFails;

	static double EdenOld;
	double 
		/* upper and lower limits to the range of the change we want to consider */
		EdenUpperLimit,
		EdenLowerLimit ,
		factor ,
		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 */

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

	LoopEden = 0;
	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;

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

	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;

		/* remember the old electron density for possible printout */
		EdenOld = phycon.eden;

		/* converge the current level of ioinization */
		ConvIoniz();

		/* ESUM sets variable EdenTrue to true elec den but does not change eden */
		esum();

		if( !conv.lgConvIoniz )
		{
			++nConvIonizFails;
		}

		/* following block is to not let electron density change by
		 * too much
		 * ElecDen.EdenError is allowed error
		 * the default value of EdenError is 0.01 in zerologic
		 * is changed with the SET EdenError command
		 * EdenOld was incoming value of eden
		 * EdenTrue is correct value from esum 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(phycon.eden-ElecDen.EdenTrue)/MIN2(phycon.eden,ElecDen.EdenTrue) >= 
			ElecDen.EdenError )
		{
			/* diference is assumed and true electron density is larger
			 * than tolerance, we are not yet converged */
			conv.lgConvEden = FALSE;
			strcpy( conv.chConvIoniz, "Ne big chg" );

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

				/* was the sign of the change in the electron density changing,
				 * or has it ever changed? */
				if( conv.lgConvEden || PreviousChange*(ElecDen.EdenTrue-phycon.eden) < 0. )
				{
					/* remember that oscillations are happening  */
					conv.lgEdenOscl = TRUE;

					/* changes in electron density are changing sign - be careful,
					 * factor multiplies the error tolerance on the electron density,
					 * largest change allowed is related to this */
					factor = 0.5;
				}
				else
				{
					/* stable so far . .. */
					/* factor = 2.;*/
					/* had been 2, change to 1 stopped oscillations from developing in 
					 * loc blr grid */
					factor = 1.;
				}

				/* use larger factor if not oscillating and change in eden is large */
				if( !conv.lgEdenOscl && 
					fabs( (ElecDen.EdenTrue-phycon.eden)/phycon.eden) > 3.*ElecDen.EdenError)
				{
					/* use 3x larger factor */
					factor *= 3.;
				}

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

				/* is the correct electron density - is it too different? */
				EdenUpperLimit = phycon.eden * (1. + ElecDen.EdenError*factor);
				EdenLowerLimit = phycon.eden / (1. + ElecDen.EdenError*factor);

				/* get the new electron density */
				if( ElecDen.EdenTrue > EdenUpperLimit )
				{
					/* this branch, proposed eden too big */
					phycon.eden = (float)EdenUpperLimit;
				}
				else if( ElecDen.EdenTrue < EdenLowerLimit )
				{
					/* proposed eden too small */
					phycon.eden = (float)EdenLowerLimit;
				}
				else
				{
					/* eden was within fac of the correct value, use it */
					phycon.eden = (float)ElecDen.EdenTrue;
				}

				/* this is set by 'set eden' command */
				if( ElecDen.EdenSet > 0. )
				{
					phycon.eden = ElecDen.EdenSet;
				}

				if( trace.lgTrace && trace.lgNeBug )
				{
					fprintf( ioQQQ, 
						"     ConvEdenIoniz, chg ne to%10.3e EdenOld%10.3e ratio%10.3e EdenTrue=%10.3e\n", 
					  phycon.eden, EdenOld, phycon.eden/EdenOld,  ElecDen.EdenTrue );
				}
			}
			/* we now have the proposed new electron density */

			/* update density and temperature related variables */
			tfidle();
		}

		if( trace.lgTrConvg>2 )
		{
			fprintf( ioQQQ, "     ConvEdenIoniz%4ld new ne=%12.4e true=%12.4e", 
			  LoopEden, phycon.eden ,  ElecDen.EdenTrue );

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

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

		/* 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 && nConvIonizFails==0)
			/* double the limit, but only one time, and only if no oscillations */
			LoopLimit *= 2;

	}	while( !conv.lgConvEden && LoopEden < LoopLimit );

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

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