/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*RT_tau_reset after first iteration, updates the optical depths, mirroring this
 * routine but with the previous iteration's variables */
#include "cddefines.h"
#include "taulines.h"
#include "trace.h"
#include "iso.h"
#include "rfield.h"
#include "opacity.h"
#include "h2.h"
#include "mole.h"
#include "geometry.h"
#include "dense.h"
#include "atomfeii.h"
#include "colden.h"
#include "rt.h"

/* ====================================================================== */
/*RT_tau_reset update total optical depth scale, 
 * called after iteration is complete */
void RT_tau_reset(void)
{
	long int i, 
	  ipHi, 
	  ipISO,
	  nelem, 
	  ipLo;

	double WeightNew;

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

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

	/* pumping of CaH */
	opac.tpcah[1] = opac.tpcah[0];
	opac.tpcah[0] = opac.taumin;

	/* save column densities of H species */
	for( i=0; i<NCOLD; ++i )
	{
		colden.colden_old[i] = colden.colden[i];
	}
	for( i=0; i < NUM_HEAVY_MOLEC; i++ )
	{
		co.hevcol_old[i] = co.hevcol[i];
	}

	/* ======================================================================== */
	/* must take average of old and new optical depths - were old in place? */
	if( iteration <= 1 )
	{
		/* this is first pass */
		WeightNew = 1.;
	}
	else
	{
		WeightNew = 0.75;
	}

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

	opac.telec = opac.taumin;
	opac.thmin = opac.taumin;

	/* do optical depths in extra Lyman lines */
	for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* the main set of all lines in the model atoms */
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				for( ipLo=0; ipLo < (iso.numLevels_max[ipISO][nelem] - 1); ipLo++ )
				{
					for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
					{
						/*RT_line_one_tau_reset computes average of old and new optical depths 
						* for new scale at end of iter */
 						RT_line_one_tau_reset(&EmisLines[ipISO][nelem][ipHi][ipLo],WeightNew);
					}
				}
				/* the extra Lyman lines */
				for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
				{
					iso.ExtraLymanLines[ipISO][nelem][ipHi].TauIn = opac.taumin;
				}
			}
		}
	}

	/* >>>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( opac.lgCaseB )
	{
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				float f;
				/* La may be case B, tlamin set to 1e9 by default with case b command */
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn = opac.tlamin;
				/* >>>chng 99 nov 22, did not reset TauCon */
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauCon = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
				EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauTot = 
				  2.f*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
				f = opac.tlamin/EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].opacity;

				for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
				{
					EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn = 
						f*EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].opacity;
					/* reset line optical depth to continuum source */
					EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauCon = EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauIn;
					EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].TauTot = 
					  2.f*EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].TauIn;
				}
			}
		}

		/* now do helium like sequence - different since collapsed levels 
		 * all go to ground */
		for( nelem=1; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				double Aprev;
				float ratio;
				/* La may be case B, tlamin set to 1e9 by default with case b command */
				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn = opac.tlamin;

				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;

				EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauTot = 
				  2.f*EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].TauIn;

				ratio = opac.tlamin/EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].opacity;

				/* this will be the trans prob of the previous lyman line, will use this to 
				 * find the next one up in the series */
				Aprev = EmisLines[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Aul;

				i = ipHe2p1P+1;
				/* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess
				 * which are which - this will work for any number of levels */
				while( i < iso.numLevels_max[ipHE_LIKE][nelem] )
				/*while( i < N_HE_LYMAN && ipHeLyman[i] < iso.numLevels_max[ipHE_LIKE][nelem] )*/
				{
					/* >>chng 02 mar 19 use proper test for resonance collapsed lines */
					/*if( EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul> Aprev/10. )*/
					if( EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul> Aprev/10. ||
						iso.quant_desig[ipHE_LIKE][nelem][i].s < 0 )
					{
						Aprev = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul;
						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn = 
							ratio*EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].opacity;
						/* reset line optical depth to continuum source */
						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauCon = EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn;
						EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauTot = 
						  2.f*EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn;
						/*fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",nelem, i, 
							EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul, EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn );*/
					}
					++ i;
				}
			}
		}
	}

	/* all heavy element lines */
	for( i=1; i <= nLevel1; i++ )
	{
		RT_line_one_tau_reset(&TauLines[i],WeightNew);
		/* >>chng 96 jul 06 following sanity check added */
		ASSERT( fabs(TauLines[i].TauIn) <= fabs(TauLines[i].TauTot) );
	}

	/* zero out fine opacity fine grid fine mesh array */
	memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(float) );

	/* all level 2 heavy element lines */
	for( i=0; i < nWindLine; i++ )
	{
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			RT_line_one_tau_reset(&TauLine2[i],WeightNew);
		}
	}

	/* all hyperfine structure lines */
	for( i=0; i < nHFLines; i++ )
	{
		RT_line_one_tau_reset(&HFLines[i],WeightNew);
	}

	/* inner shell lines */
	for( i=0; i < nUTA; i++ )
	{
		if( UTALines[i].Aul > 0. )
		{
			RT_line_one_tau_reset(&UTALines[i],WeightNew);
		}
	}

	/* co carbon monoxide lines */
	for( i=0; i < nCORotate; i++ )
	{
		RT_line_one_tau_reset(&C12O16Rotate[i],WeightNew);
		RT_line_one_tau_reset(&C13O16Rotate[i],WeightNew);
	}

	/* the large H2 molecule */
	H2_RT_tau_reset();

	/* large FeII atom */
	FeII_RT_tau_reset();

	if( opac.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
			 * TauAbsGeo(i,2) = 2.*TauAbsFace(i) */
			opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
			/* TauScatGeo(i,2) = 2.*TauScatFace(i) */
			opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
			opac.TauScatFace[i] = opac.taumin;
			opac.TauAbsFace[i] = opac.taumin;
		}
	}
	else if( geometry.lgSphere )
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			/* [1] is total optical depth from previous iteration,
			 * [0] is optical depth at current position */
			opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
			opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
			opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
			opac.TauScatGeo[0][i] = opac.TauScatFace[i];
			opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
			opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
			/* TauAbsFace and TauScatFace are abs and sct optical depth to ill face */
			opac.TauScatFace[i] = opac.taumin;
			opac.TauAbsFace[i] = opac.taumin;
		}
	}
	else
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
			opac.TauTotalGeo[0][i] = opac.taumin;
			opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
			opac.TauAbsGeo[0][i] = opac.taumin;
			opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
			opac.TauScatGeo[0][i] = opac.taumin;
			opac.TauScatFace[i] = opac.taumin;
			opac.TauAbsFace[i] = opac.taumin;
		}
	}

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

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

