/*ionize main routine to drive ionization solution for all species,
 * it is only called by ConvIoniz */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "toler.h"
#include "tseton.h"
#include "trace.h"
#include "h.h"
#include "heat.h"
#include "phycon.h"
#include "elmton.h"
#include "he.h"
#include "tcool.h"
#include "secondaries.h"
#include "highen.h"
#include "helium.h"
#include "elecden.h"
#include "makechartran.h"
#include "sumheat.h"
#include "rt.h"
#include "ionheavy.h"
#include "sumcontinuum.h"
#include "addopac.h"
#include "chksumcon.h"
#include "freeht.h"
#include "iiibod.h"
#include "grain.h"
#include "dielsupres.h"
#include "lgconverg.h"
#include "hemetalots.h"
#include "ionrange.h"
#include "ionfracs.h"
#include "opac.h"
#include "coolr.h"
#include "prtotsrate.h"
/* this contains declaration for cddrive */
#include "comole.h"
#include "trimstages.h"
#include "converge.h"

void ionize( 
	/* this is zero the first time ionize is called by convIoniz,
	 * counts number of call thereafter */
	 long loopi )
{
	double HeatOld;
	static double SecondOld;
	long int ipZ,
		ipOTSchange;
	static double SumOTS=0. , OldSumOTS[2] ;
	double RatioOTSInc;

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

	if( loopi==0 )
	{
		/* these will be used to look for oscillating ots rates */
		OldSumOTS[0] = 0.;
		OldSumOTS[1] = 0.;
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"   IONIZE called.  TE:%9.3e  HI:%8.2e  HII:%8.2e  Ne:%9.3e  CSUP:%8.2e\n", 
		  phycon.te, h.hi, h.hii, phycon.eden, Secondaries.csupra);
	}

	/* this routine is in sumheat, and zeroes out the heating array */
	ZeroHeat();

	/* if this is very first call, say not converged, since no meaningful way to
	 * check on changes in quantities.  this counter is false only on very first
	 * call, when equal to zero */
	if( !conv.nTotalIoniz )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "first call" );
	}

	/* this will be flag to check whether any ionization stages
	 * were trimed down */
	conv.lgIonStageTrimed = FALSE;

	/* must redo photoionizaiton rates for all species on second and later tries */
	/* always reevaluate the rates when . . . */
	if( 
		/* opac.lgOpacStatic (usually TRUE), is set FALSE with no static opacity command */
		!opac.lgOpacStatic || 
		/* we are in search mode */
		conv.lgSearch || 
		/* this is the first call to this zone */
		conv.nPres2Ioniz == 0 )
	{
		/* we need to redo ALL photoionization rates */
		opac.lgRedoStatic = TRUE;
		/* this adjusts the lowest and highest stages of ionization we will consider,
		 * only safe to call when lgRedoStatic is true since this could lower the 
		 * lowest stage of ionization, which needs all its photo rates */

		/* conv.nTotalIoniz is only 0 (FALSE) on the very first call to ionize,
		 * when we do not know what the ionization distribution is, since not yet done */
		if( conv.nTotalIoniz )
		{
			/* we will never trim hydrogen itself (ipZ==0) down.  
			 * Criteria for helium are different
			 * from rest since He+ CT with CO is so important - keep He+ going
			 * far longer than we would expect */
			if( elmton.lgElmtOn[ipHE] )
			{
				if( IonRange.IonHigh[ipHE] == 3 )
				{
					/* fully stripped helium is present - how abundant is it? 
					 * this limit is the one used in C90, probably not carefully thought out */
					if( xIonFracs[3][ipHE]/xIonFracs[0][ipHE] < 1e-6F )
					{
						/* decrement the counter to the highest stage, and zero its abundance */
						--IonRange.IonHigh[ipHE];
						xIonFracs[3][ipHE] = 0.;
						he.heiii = 0.;
					}
				}
				else if( IonRange.IonHigh[ipHE] == 2 )
				{
					/* no stripped ion, but check on single ion.  this is the limit that
					 * was present in c90, small for CO destruction */
					if( xIonFracs[2][ipHE]/xIonFracs[0][ipHE] < 1e-12F )
					{
						/* decrement the counter to the highest stage, and zero its abundance */
						--IonRange.IonHigh[ipHE];
						xIonFracs[2][ipHE] = 0.;
						he.heii = 0.;
					}
				}
			}

			/* TrimStages only used for Li and heavier, not H or He */
			for( ipZ=2; ipZ<LIMELM; ++ipZ )
			{
				if( elmton.lgElmtOn[ipZ] )
				{
					/* TrimStages will set conv.lgIonStageTrimed TRUE is any ion has its
					 * lowest stage of ionization dropped */
					TrimStages(ipZ);
				}
			}
		}

		/* now check if anything trimed down */
		if( conv.lgIonStageTrimed )
		{
			/* something was trimed down, so say that ionization not yet stable */
			/* say that ionization has not converged, secondaries changed by too much */
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "IonTrimmed" );
		}
	}
	else
	{
		opac.lgRedoStatic = FALSE;
	}

	/* do all hydrogenic species */
	Hydrogenic();

	/* evaluate Compton heating, bound E compton, cosmic rays */
	highen();

	/* evaluate free-free heating rate */
	freeht();

	/* find corrections for three-body rec - collisional ionization */
	iiibod();

	/* helium ionization balance */
	helium();

	/* find grain temperature, charge, and gas heating rate */
	grain();

	/* deduce dielctronic supression factors */
	DielSupres();

	/* fill in master charge transfer array, and zero out some arrays that track effects */
	MakeCharTran();

	/* do the ionization balance */
	IonLithi();
	IonBeryl();
	IonBoron();
	IonCarbo();
	IonOxyge();
	/* do carbon monoxide after oxygen */
	codriv();
	IonNitro();
	IonFluor();
	IonNeon();
	IonSodiu();
	IonMagne();
	IonAlumi();
	IonSilic();
	IonPhosi();
	IonSulph();
	IonChlor();
	IonArgon();
	IonPotas();
	IonCalci();
	IonScand();
	IonTitan();
	IonVanad();
	IonChrom();
	IonManga();
	IonIron();
	IonCobal();
	IonNicke();
	IonCoppe();
	IonZinc();

	/* lgConverg is a function to check whether ionization has converged */
	/* check whether oxygen ionization changed by more than relative error
	 * given by second number */
	if( !lgConverg(1, 0.05 ) )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "Helium Ion" );
	}

	if( !lgConverg(7, 0.15 ) )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "Oxygen Ion" );
	}

	/* check whether iron ionization changed by much */
	if( !lgConverg(25 , 0.15 ) )
	{
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "Iron Ioniz" );
	}

	/* get total cooling, this can add heat */
	coolr(&tcool.ctot );
		
	/* get total heating rate */
	HeatOld = heat.htot;

	/* SumHeat will update ElecFrac, 
	 * secondary electron ionization and excitation efficiencies,
	 * and sum up current secondary rates - remember old values before we enter */
	SecondOld = Secondaries.csupra;

	SumHeat();

	/* test whether we can just set the rate to the new one, or whether we should
	 * take average and do it again.  Secondaries.sec2total was set in hydrogenic, and
	 * is ratio of secondary to total hydrogen destruction rates */
	if( (Secondaries.sec2total>0.001) && 
		fabs( 1. - SecondOld/MAX2(SMALLFLOAT,Secondaries.csupra) ) > 0.1 && 
		SecondOld > 0. &&  Secondaries.csupra > 0.)
	{
		/* say that ionization has not converged, secondaries changed by too much */
		conv.lgConvIoniz = FALSE;
		strcpy( conv.chConvIoniz, "SecIonRate" );
		conv.BadConvIoniz[0] = SecondOld;
		conv.BadConvIoniz[1] = Secondaries.csupra;
	}

	if( HeatOld > 0. && !tseton.lgTSetOn )
	{
		/* check if heating has converged - toler is final match */
		if( fabs(1.-heat.htot/HeatOld) > tolerCom.toler*.5 )
		{
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "Big d Heat" );
			conv.BadConvIoniz[0] = HeatOld;
			conv.BadConvIoniz[1] = heat.htot;
		}
	}

	/* evaluate current opacities, AddOpac is called only here during actual calculation,
	 * although also called in ConvInitTemp if upper limit to continuum increased */
	AddOpac();

	/* conv.nTotalIoniz was set to zero in zero and reset, and incremented below,
	 * and so counts the number of times this routine has been called during current iteration.  
	 * the very first time we are called this is zero, false, so ots routines not called, 
	 * since both opacities and dest prob are junk */
	if( conv.nTotalIoniz > 0 )
	{
		/* do hydrogenic OTS rates for H lines and all continua since
		 * we now have ionization balance of all species.  Note that this is not
		 * entirely self-consistnet, since dest probs here are not the same as
		 * the ones used in the model atoms.  Problems??  if near convergence
		 * then should be nearly identical */
		/* line dest rates are evaluated in HydroPesc */
		HydroOTS();
		HeMetalOTS();
	}

	/* remember old ots rates */
	OldSumOTS[0] = OldSumOTS[1];
	OldSumOTS[1] = SumOTS;

	/* now update several components of the continuum, this only happens after
	 * we have gone through entire soln for this zone at least one time. 
	 * this var is 0 on very first call to this zone 
	 * done due to wild oscillation in ots on first call */
	SumContinuum(&SumOTS , &RatioOTSInc , &ipOTSchange );

	/* now check whether the ots rates changed */
	if( SumOTS> 0. )
	{
		/* the ots rate must be converged to the error in the electron density */
		/*fprintf(ioQQQ,"NewOldOTS= %g %g\n",1. - OldSumOTS / SumOTS , RatioOTSInc);*/
		/*if( fabs(1.-OldSumOTS/SumOTS) > 0.10 )*/
		/* how fine should this be converged?? orginally had above, 10%, but take
		 * smaller ratio?? */
		if( fabs(1.-OldSumOTS[1]/SumOTS) > ElecDen.EdenError )
		{
			/* this branch, ionization not converged due to large change in ots rates.
			 * check whether ots rates are oscillating, if loopi > 1 so we have enough info*/
			if( loopi > 1 )
			{
				/* here we have three values, are they changing sign? */
				if( (OldSumOTS[0]-OldSumOTS[1]) * ( OldSumOTS[1] - SumOTS ) < 0. )
				{
					/* ots rates are oscillating, if so then use small fraction of new dest rates 
					 * when updating Lyman line dest probs in HydroPesc*/
					conv.lgOscilOTS = TRUE;
				}
				/*else
				{
					 not oscillating, take mean of old and new 
					 chck effects of not setting false within zone conv.lgOscilOTS = FALSE;
				}*/
			}
			/*else
			{
				 not enough information at present iteration 
				conv.lgOscilOTS = FALSE;
			}*/
			/*fprintf(ioQQQ," ots change, oscil=%li\n", conv.lgOscilOTS);*/

			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "OTSRatChng" );
			conv.BadConvIoniz[0] = OldSumOTS[1];
			conv.BadConvIoniz[1] = SumOTS;
			/* uncomment this to print contributors to ots rate */
			/*PrtOtsRate(SumOTS*0.1 , 'l' );*/
		}
	}

	/* this is only a sanity check that the summed continua accurately reflect
	 * all of the individual components.  Only include this when NDEBUG is not set,
	 * we are in not debug compile */
#	if !defined(NDEBUG)
	ChkSumCon(0);
#	endif

	/* option to print OTS continuum with TRACE OTS */
	if( trace.lgTrace && trace.lgOTSBug )
	{
		/* find ots rates here, so we only print fraction of it,
		 * SumOTS is both line and continuum contributing to ots, and is mult by opacity */
		SumContinuum(&SumOTS , &RatioOTSInc , &ipOTSchange);
		/* number is weakest rate to print */
		PrtOtsRate(SumOTS*0.05 , 'b' );
	}

	/* now reevaluate only destruction probabilities */
	RTMake(FALSE);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "   IONIZE retrn, HI:%10.3e  HII:%10.3e  EDEN:%10.3e HTOT:%10.3e lgIonDone=%3c\n", 
		  h.hi, h.hii, phycon.eden, heat.htot, TorF(conv.lgConvIoniz) );
	}
	
	/* this counts number of times we are called by PresChng, in this zone */
	++conv.nPres2Ioniz;

	/* this counts how many times we have been called in this model,
	 * located here at bottom of routine so that number is false on very first
	 * call, since set to 0 in zero */
	++conv.nTotalIoniz;

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