/*update update total optical depth scale, called after iteration is complete */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "nhe1lvl.h"
#include "deptau.h"
#include "hydrogenic.h"
#include "caseb.h"
#include "double.h"
#include "xry.h"
#include "he3tau.h"
#include "itercnt.h"
#include "he1nxt.h"
#include "taumin.h"
#include "he1tau.h"
#include "rfield.h"
#include "showme.h"
#include "phe1lv.h"
#include "trace.h"
#include "stopcalc.h"
#include "touton.h"
#include "elmton.h"
#include "sphere.h"
#include "opac.h"
#include "otsmin.h"
#include "pmpcah.h"
#include "nhe1.h"
#include "autoit.h"
#include "taucon.h"
#include "colden.h"
#include "taddlya.h"
#include "punconv.h"
#include "noconverge.h"
#include "zonecnt.h"
#include "phycon.h"
#include "pop371.h"
#include "receff.h"
#include "tav.h"
#include "update.h"

void update(void)
{
	int lgConverged;

	long int i, 
	  ipHi, 
	  ipZ, 
	  ipLo;

	double AverFac, 
	  TauOff, 
	  f1, 
	  f2, 
	  fhe1;

	static double OldEfrac=0.;

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

	/* save column density */
	coldenCom.ColDenSav[IterCnt.iter-1] = coldenCom.colden[IPCOLUMNDENSITY-1];

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " UPDATE estimating new optical depths\n" );
		if( trace.lgHBug && trace.lgHBugFull )
		{
			fprintf( ioQQQ, " New Hydrogen outward optical depths:\n" );
			for( ipHi=1; ipHi <= nhlevel; ipHi++ )
			{
				fprintf( ioQQQ, "%3ld", ipHi );
				for( ipLo=0; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ, "%10.2e", 
						HydroLines[trace.ipZTrace][ipHi][ipLo].TauIn );
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	/* pumping of CaH */
	AverFac = log(MAX2(1e-37,pmpcah.tpcah[0]));
	f1 = sqrt(MAX2(1.,AverFac));
	f2 = pow(6.1e-5*pmpcah.tpcah[0],0.5);
	pmpcah.conca = (float)MAX2(f1,f2 );

	AverFac = log(MAX2(1e-37,HydroLines[0][6][2].TauIn));
	AverFac = sqrt(MAX2(1.,AverFac));
	f1 = pow(1.7e-6*MAX2(1e-37,HydroLines[0][6][2].TauIn),0.5);
	pmpcah.conhe = (float)MAX2(AverFac, f1);

	pmpcah.tpcah[1] = pmpcah.tpcah[0];
	pmpcah.tpcah[0] = tauminCom.taumin;

	/* this is an option to keep iterating until it converges
	 * iterate to convergence option
	 * autocv is percentage difference in optical depths allowed,
	 * =0.20 in block data */
	lgConverged = TRUE;
	strcpy( NoConverge.chNotConverged, "Converged!" );

	if( IterCnt.iter > 1 && autoit.lgAutoIt )
	{
		for( ipZ=0; ipZ < LIMELM; ipZ++ )
		{
			if( elmton.lgElmtOn[ipZ] )
			{
				/* check both H-alpha and Ly-alpha - only if balmer lines thick */
				/* following checks if Ha optical deth significant */
				if( HydroLines[ipZ][3][2].TauIn > 0.5 )
				{
					if( fabs(HydroLines[ipZ][IP2P][IP1S].TauTot/
					  (HydroLines[ipZ][IP2P][IP1S].TauIn*double_.DoubleTau)-1.) > 
					  autoit.autocv )
					{
						/* not converged to within AUTOCV, normally 15 percent */
						lgConverged = FALSE;
						strcpy( NoConverge.chNotConverged, "H Lya     " );
						if( PunConv.lgPunConv )
						{
							fprintf( PunConv.ipPunConv, " H Lya\n" );
						}
					}

					if( fabs(HydroLines[ipZ][3][2].TauTot/
					  (HydroLines[ipZ][3][2].TauIn*double_.DoubleTau)-1.) > 
					  autoit.autocv )
					{
						/* not converged to within AUTOCV, normally 15 percent */
						lgConverged = FALSE;
						strcpy( NoConverge.chNotConverged, "H B a     " );

						if( PunConv.lgPunConv )
						{
							fprintf( PunConv.ipPunConv, " H B a\n" );
						}
					}
				}
			}
		}

		if( (he1tauCOM.he1tau[0][1] - TAddLya.TAddHeI) > 0.5 )
		{
			if( fabs(he1tauCOM.he1tau[0][1]*double_.DoubleTau/he1tauCOM.he1lim[0][1]-1.) > 
				autoit.autocv )
			{
				lgConverged = FALSE;
				strcpy( NoConverge.chNotConverged, "He1 Lya   " );

				if( PunConv.lgPunConv )
				{
					fprintf( PunConv.ipPunConv, " He1 Lya\n" );
				}

			}
		}

		if( PunConv.lgPunConv && lgConverged )
		{
			fprintf( PunConv.ipPunConv, " converged\n" );
		}

		/* lower limit to number of iterations if converged */
		if( lgConverged )
			IterCnt.itermx = MIN2(IterCnt.itermx,IterCnt.iter);

		/* >>chng 96 dec 20, moved following to within if on lgAutoIt
		 * this is test for stopping on first zone */
		if( phycon.te < StopCalc.tend && ZoneCnt.nzone == 1 )
		{
			lgConverged = TRUE;
			strcpy( NoConverge.chNotConverged, "          " );
			IterCnt.itermx = MIN2(IterCnt.itermx,IterCnt.iter);
		}
	}

	/* must take average of old and new optical depths - were old in place?
	 * */
	if( IterCnt.iter <= 1 )
	{
		/* this is first pass */
		f1 = 1.;
	}
	else
	{
		/* >>chng 96 june 6, kill logic for oscillations, do not use
		 * optical depths this was too conservative.  instead check whether
		 * i-front breaking out by looking at electron fractions */
		if( (OldEfrac - 0.98)*(phycon.ElecFrac - 0.98) < 0. )
		{
			f1 = 0.2;
		}
		else
		{
			/* no oscillations, save to speed up convergence */
			f1 = 0.75;
		}
	}

	/* will use this to check on electron density oscillations */
	OldEfrac = phycon.ElecFrac;

	AverFac = f1;
	f2 = 1. - f1;

	/* static is false unless sphere static set */
	if( sphere.lgStatic )
	{
		for( ipHi=1; ipHi < NHE1LVL; ipHi++ )
		{
			/* DoubleTau is usually 1, set to 2 with the DoubleTau option to
			 * simulate two-sided photoionization */
			if( IterCnt.iter == 1 )
			{
				TauOff = TAddLya.TAddHeI*he1tauCOM.he1opc[0][ipHi]/he1tauCOM.he1opc[0][1];
			}
			else
			{
				TauOff = 0.;
			}

			he1tauCOM.he1lim[0][ipHi] = (float)pow(10.,log10(MAX2(tauminCom.taumin,
			  (he1tauCOM.he1tau[0][ipHi]-TauOff)*double_.DoubleTau))*
			  f1 + log10(MAX2(tauminCom.taumin,he1tauCOM.he1lim[0][ipHi]))* f2);

			if( he1tauCOM.he1lim[0][ipHi] <= 0. )
			{
				/* this probably due to underflow above */
				he1tauCOM.he1lim[0][ipHi] = (float)(TauOff/1e6);
			}
			he1nxtCOM.he1nxt[0][ipHi] = (float)(he1tauCOM.he1lim[0][ipHi]/2.);
			he1tauCOM.he1tau[0][ipHi] = (float)(he1tauCOM.he1lim[0][ipHi]/2.);
		}
	}
	else
	{
		for( ipHi=1; ipHi < NHE1LVL; ipHi++ )
		{
			/* DoubleTau is usually 1, set to 2 with the DoubleTau option to
			 * simulate two-sided photoionization */
			if( IterCnt.iter == 1 )
			{
				TauOff = TAddLya.TAddHeI*he1tauCOM.he1opc[0][ipHi]/he1tauCOM.he1opc[0][1];
			}
			else
			{
				TauOff = 0.;
			}
			/* he1lim(ipHi,1) = (he1tau(ipHi,1)-TauOff)*f1*DoubleTau + he1lim(ipHi,1)*f2 */
			he1tauCOM.he1lim[0][ipHi] = (float)pow(10.,log10(MAX2(tauminCom.taumin,
			  (he1tauCOM.he1tau[0][ipHi]-TauOff)*double_.DoubleTau))*
			  f1 + log10(MAX2(tauminCom.taumin,he1tauCOM.he1lim[0][ipHi]))*
			  f2);
			if( he1tauCOM.he1lim[0][ipHi] <= 0. )
			{
				/* this probably due to underflow above */
				he1tauCOM.he1lim[0][ipHi] = (float)(TauOff/1e6);
			}
			he1nxtCOM.he1nxt[0][ipHi] = 
				(float)MIN2(tauminCom.taumin,he1tauCOM.he1lim[0][ipHi]/2.);
			he1tauCOM.he1tau[0][ipHi] = 
				(float)MIN2(tauminCom.taumin,he1tauCOM.he1lim[0][ipHi]/2.);
		}
	}

	if( sphere.lgStatic )
	{
		for( ipLo=1; ipLo < (NHE1LVL - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < NHE1LVL; ipHi++ )
			{
				/* he1lim(ipHi,ipLo) = he1tau(ipHi,ipLo)*f1*DoubleTau + he1lim(ipHi,ipLo) * f2 */
				he1tauCOM.he1lim[ipLo][ipHi] = (float)pow(10.,log10(MAX2(tauminCom.taumin,
				  (he1tauCOM.he1tau[ipLo][ipHi])*double_.DoubleTau))*f1 + 
				  log10(MAX2(tauminCom.taumin,he1tauCOM.he1lim[ipLo][ipHi]))*f2);

				he1nxtCOM.he1nxt[ipLo][ipHi] = (float)(he1tauCOM.he1lim[ipLo][ipHi]/2.);
				he1tauCOM.he1tau[ipLo][ipHi] = (float)(he1tauCOM.he1lim[ipLo][ipHi]/2.);
			}
		}
	}
	else
	{
		for( ipLo=1; ipLo < (NHE1LVL - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < NHE1LVL; ipHi++ )
			{
				/* he1lim(ipHi,ipLo) = he1tau(ipHi,ipLo)*f1*DoubleTau + he1lim(ipHi,ipLo) * f2 */
				he1tauCOM.he1lim[ipLo][ipHi] = (float)pow(10.,log10(MAX2(tauminCom.taumin,
				  (he1tauCOM.he1tau[ipLo][ipHi])*double_.DoubleTau))*f1 + 
				  log10(MAX2(tauminCom.taumin,he1tauCOM.he1lim[ipLo][ipHi]))*f2);
				he1nxtCOM.he1nxt[ipLo][ipHi] = 
					(float)MIN2(tauminCom.taumin,he1tauCOM.he1lim[ipLo][ipHi]/2.);
				he1tauCOM.he1tau[ipLo][ipHi] = 
					(float)MIN2(tauminCom.taumin,he1tauCOM.he1lim[ipLo][ipHi]/ 2.);
			}
		}
	}

	/* force hydrogen outward optical depths to be on */
	touton.lgTauOutOn = TRUE;

	/* possibly fix La optical depths
	 * tlamin is set to large number when case b set */
	if( (tauminCom.tlamin > 1e-5 && (!sphere.lgStatic)) || caseb.lgCaseb )
	{
		fhe1 = tauminCom.tlamin/he1tauCOM.he1opc[0][1];
		he1tauCOM.he1tau[0][1] = tauminCom.tlamin;
		he1tauCOM.he1lim[0][1] = tauminCom.tlamin*2.f;
		he1nxtCOM.he1nxt[0][1] = tauminCom.tlamin;
		for( ipHi=2; ipHi < NHE1LVL; ipHi++ )
		{
			he1tauCOM.he1tau[0][ipHi] = (float)(fhe1*he1tauCOM.he1opc[0][ipHi]);
			he1tauCOM.he1lim[0][ipHi] = he1tauCOM.he1tau[0][ipHi]*2.f;
			he1nxtCOM.he1nxt[0][ipHi] = he1tauCOM.he1tau[0][ipHi];
		}
	}

	/* following sets (1)=inward optical depth to TAUMIN,
	 * (2)=outward = (log) mean of old optical depths,
	 * sets (3) to escap prob for first zone */
	for( ipHi=0; ipHi < NHE3TAU; ipHi++ )
	{
		tav(&he3tau[ipHi],AverFac);
	}

	tav(&t206,AverFac);

	if( sphere.lgSphere )
	{
		phe1lv.he1rec[1][0] = (float)MAX2(0.,otsminCom.otsmin);
		phe1lv.he1rec[1][1] = 1.;
	}
	else
	{
		phe1lv.he1rec[1][0] = (float)receff(nhe1Com.nhe1[0]);
	}

	taucon.telec = tauminCom.taumin;
	taucon.thmin = tauminCom.taumin;

	/* hydrogenic species */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		if( elmton.lgElmtOn[ipZ] )
		{
			for( ipLo=0; ipLo <= (nhlevel - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
				{
					/*tav computes average of old and new optical depths 
					 * for new scale at end of iter */
					tav(&HydroLines[ipZ][ipHi][ipLo],0.5);
				}
			}
		}
	}

	/* >>>chng 99 nov 11 did not have case b for hydrogenic species on second and
	 * higher iterations */
	/* option to clobber these taus for lyman lines, if case b is set */
	if( caseb.lgCaseb )
	{
		for( ipZ=0; ipZ < LIMELM; ipZ++ )
		{
			if( elmton.lgElmtOn[ipZ] )
			{
				double f;
				/* La may be case B, tlamin set to 1e9 by default with case b command */
				HydroLines[ipZ][IP2P][IP1S].TauIn = tauminCom.tlamin;
				/* >>>chng 99 nov 22, did not reset TauCon */
				HydroLines[ipZ][IP2P][IP1S].TauCon = HydroLines[ipZ][IP2P][IP1S].TauIn;
				HydroLines[ipZ][IP2P][IP1S].TauTot = 
				  (float)(2.*HydroLines[ipZ][IP2P][IP1S].TauIn);
				f = tauminCom.tlamin/HydroLines[ipZ][IP2P][IP1S].opacity;

				for( ipHi=3; ipHi <= nhlevel; ipHi++ )
				{
					HydroLines[ipZ][ipHi][IP1S].TauIn = 
						(float)(f*HydroLines[ipZ][ipHi][IP1S].opacity);
					/* reset line optical depth to continuum source */
					HydroLines[ipZ][ipHi][IP1S].TauCon = HydroLines[ipZ][ipHi][IP1S].TauIn;
					HydroLines[ipZ][ipHi][IP1S].TauTot = 
					  (float)(2.*HydroLines[ipZ][IP2P][IP1S].TauIn);
				}
			}
		}
	}

	/* all heavy element lines */
	for( i=1; i <= nLevel1; i++ )
	{
		tav(&TauLines[i],AverFac);
		/* >>chng 96 july 6 following sanity check added
		 *begin sanity check */
		if( fabs(TauLines[i].TauIn) > 
			fabs(TauLines[i].TauTot) )
		{
			fprintf( ioQQQ, " UPDATE finds TauIn>TauOut, line=%5ld%10.2e%10.2e\n", 
			  i, TauLines[i].TauIn, TauLines[i].TauTot );
			ShowMe();
			puts( "[Stop in update]" );
			exit(1);
		}
		/*end sanity check */
	}

	/* all heavy element lines not in cooling */
	for( i=0; i < nWindLine; i++ )
	{
		tav(&TauLine2[i],AverFac);
	}

	/* Verner's FeII atom */
	if( FeII.lgFeIION )
	{
		FeIITauAver();
	}

	if( caseb.lgCaseb )
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			/* DEPABS and SCT are abs and sct optical depth for depth only
			 * we will not change total optical depths, just reset inner to half
			 * tauabs(i,2) = 2.*depabs(i) */
			opac.tauabs[0][i] = (float)(opac.tauabs[1][i]/2.);
			/* tausct(i,2) = 2.*depsct(i) */
			opac.tausct[0][i] = (float)(opac.tausct[1][i]/2.);
			deptau.depsct[i] = tauminCom.taumin;
			deptau.depabs[i] = tauminCom.taumin;
		}
	}
	else if( sphere.lgSphere )
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			/* DEPABS and SCT are abs and sct optical depth for depth only */
			opac.tauabs[1][i] = (float)(2.*deptau.depabs[i]);
			opac.tauabs[0][i] = deptau.depabs[i];
			opac.tausct[1][i] = (float)(2.*deptau.depsct[i]);
			opac.tausct[0][i] = deptau.depsct[i];
			opac.TauTotal[1][i] = opac.tausct[1][i] + opac.tauabs[1][i];
			opac.TauTotal[0][i] = opac.tausct[0][i] + opac.tauabs[0][i];
			deptau.depsct[i] = tauminCom.taumin;
			deptau.depabs[i] = tauminCom.taumin;
		}
	}
	else
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			opac.TauTotal[1][i] = opac.TauTotal[0][i];
			opac.TauTotal[0][i] = tauminCom.taumin;
			opac.tauabs[1][i] = opac.tauabs[0][i];
			opac.tauabs[0][i] = tauminCom.taumin;
			opac.tausct[1][i] = opac.tausct[0][i];
			opac.tausct[0][i] = tauminCom.taumin;
			deptau.depsct[i] = tauminCom.taumin;
			deptau.depabs[i] = tauminCom.taumin;
		}
	}

	/* this is optical depth at x-ray point defining effective optical depth */
	xry.tauxry = opac.tauabs[0][xry.ipxry-1];

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