/*ConvPresTempEdenIoniz drive pressure, ionization, and thermal balance */
#include <string.h>
#include "cddefines.h"
#include "phycon.h"
#include "elecden.h"
#include "tcool.h"
#include "heat.h"
#include "pressure.h"
#include "trace.h"
#include "converge.h"

void ConvPresTempEdenIoniz(void)
{
	long int loop,
		LoopMax=20;
	double Error;
	float OldChange;
	int lgPresOscil;

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

	/* this will count number of times we call ionize in this zone */
	conv.nPres2Ioniz = 0;
	loop = 0;
	conv.lgConvPres = FALSE;

	/* this will be the limit, which we will increase if no oscillations occur */
	LoopMax = 20;
	OldChange = 0.;

	/* this will be flag to check for pressure oscillations */
	lgPresOscil = FALSE;

	while( (!conv.lgConvPres && loop < LoopMax) )
	{
		/* change current densities of all constituents if necessary, 
		 * PresChng evaluates lgPresOK, true if pressure is now ok */
		PresChng();

		/* heating cooling balance while doing ionization,
		 * this is where the heavy lifting is done*/
		ConvTempEdenIoniz();

		/* check whether pressure is oscillating */
		if( (pressure.pnow - pressure.pcorec ) * OldChange < 0. )
		{
			/* the sign of the change in pressure has changed, so things
			 * are oscillating.  This would be a problem */
			lgPresOscil = TRUE;
		}
		OldChange = pressure.pnow - pressure.pcorec ;

		/* only print convergence trace at this level if NOT constant density ("CDEN") */
		if( trace.lgTrConvg>0 && strcmp(pressure.chCPres,"CDEN") )
		{
			fprintf( ioQQQ, " ConvPresTempEdenIoniz%4ld Curr Pres:%10.2e Corr Pres:%10.2e Hden:%10.2e\n", 
			  loop, pressure.pnow, pressure.pcorec, phycon.hden );
		}

		/* increment loop counter */
		++loop;

		/* if we hit limit of loop, but no oscillations have happened, then we are
		 * making progress, and can keep going */
		if( loop == LoopMax && !lgPresOscil )
		{
			LoopMax = MIN2( 100 , LoopMax*2 );
		}
	}

	if( !conv.lgConvPres )
	{
		/* announce the pressure failure */
		ConvFail("pres");
	}

	/* say something if we did not converge */
	if( !conv.lgConvEden )
	{
		ConvFail("eden") ;
	}

	if( !conv.lgConvTemp )
	{
		/* announce the temperature failure */
		ConvFail("temp");
	}

	/* say something if we did not converge */
	if( !conv.lgConvIoniz )
	{
		ConvFail("ioni") ;
	}

	/* remember mean and largest errors on electron density */
	Error = fabs(phycon.eden - ElecDen.EdenTrue)/ElecDen.EdenTrue;
	conv.BigEdenError = (float)MAX2(Error , conv.BigEdenError );
	conv.AverEdenError += (float)Error;

	/* remember mean and largest errors between heating and cooling */
	Error = fabs(tcool.ctot - heat.htot) / tcool.ctot;
	conv.BigHeatCoolError = (float)MAX2(Error , conv.BigHeatCoolError );
	conv.AverHeatCoolError += (float)Error;

	/* remember mean and largest pressure errors */
	Error = fabs(pressure.pnow - pressure.pcorec) / pressure.pcorec;
	conv.BigPressError = (float)MAX2(Error , conv.BigPressError );
	conv.AverPressError += (float)Error;

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