/* 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_inc increment optical depths once per zone, called after radius_increment */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "rfield.h"
#include "trace.h"
#include "dense.h"
#include "hyperfine.h"
#include "wind.h"
#include "prt.h"
#include "h2.h"
#include "hmi.h"
#include "opacity.h"
#include "radius.h"
#include "atomfeii.h"
#include "rt.h"

/*RT_tau_inc increment optical depths once per zone, called after radius_increment */
void RT_tau_inc(void)
{

	long int i, 
	  ipHi, 
	  nelem, 
	  ipLo,
	  ipISO;

	double factor;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " RT_tau_inc called.\n" );
	}

	/* call RT_line_all one last time in this zone, to get fine opacities defined */
	RT_line_all( FALSE , TRUE);

	opac.telec += (float)(radius.drad_x_fillfac*dense.eden*6.65e-25);
	opac.thmin += (float)(radius.drad_x_fillfac*hmi.Hmolec[ipMHm]*3.9e-17*
		(1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep));

	/* prevent maser runaway */
	rt.dTauMase = 0;
	rt.mas_species = 0;
	rt.mas_ion = 0;
	rt.mas_hi = 0;
	rt.mas_lo = 0;

	/* >>chng 03 apr 11, to iso loop */
	/* loop over all lines in iso sequences */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			/* this is the parent ion, for H lines, is 1, 
			 * for He is 1 for he-like and 2 for h-like */
			int ion = nelem+1-ipISO;
			if( ion <=dense.IonHigh[nelem] )
			{
				/* must convert to physical scale for opacities,
				* get ionic abundances first */
				if( dense.xIonDense[nelem][ion] > 1e-30 )
				{
					factor = dense.xIonDense[nelem][ion];
				}
				else
				{
					/* case where almost no parent ion - this will make
					 * very large line opacity, so background dest small */
					factor = 1.;
				}
				for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
				{
					for( ipLo=0; ipLo < ipHi; ipLo++ )
					{
						/* must temporarily make ipLnPopOpc physical */
						EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc *= factor;
						/* actually do the work */
						RT_line_one_tauinc(&EmisLines[ipISO][nelem][ipHi][ipLo] ,
							ipISO , nelem , ipHi , ipLo );
						/* go back to original units so that final correction ok */
						EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc /= factor;
					}
				}
				ipLo = 0;
				/* these are the extra lyman lines */
				for( ipHi=iso.quant_desig[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
				/*for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
				{
					/* must make ipLnPopOpc physical */
					iso.ExtraLymanLines[ipISO][nelem][ipHi].PopOpc = 
						/* just use this line to get population of ground state */
						EmisLines[ipISO][nelem][3][ipLo].PopOpc *factor;

					/* actually do the work */
					RT_line_one_tauinc(&iso.ExtraLymanLines[ipISO][nelem][ipHi] ,
							-1 ,ipISO , nelem , ipHi );

					/* reset to funny population of ground state */
					iso.ExtraLymanLines[ipISO][nelem][ipHi].PopOpc = 
						EmisLines[ipISO][nelem][3][ipLo].PopOpc;
				}
			}
		}
	}

	/* increment optical depths for all heavy element lines
	 * same routine does wind and static,
	 * does not start from 0 since first line is dummy */
	/*DumpLine( &TauLines[ipT63] );*/
	for( i=1; i <= nLevel1; i++ )
	{
		RT_line_one_tauinc(&TauLines[i],
			-2 , -2 , -2 , i );
	}

	/* all lines in cooling with g-bar */
	for( i=0; i < nWindLine; i++ )
	{
		/* do not include H-like or He-like in the level two lines since
		 * these are already counted in iso sequences */
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			RT_line_one_tauinc(&TauLine2[i] ,
				-3 , -3 , -3 , i );
		}
	}

	/* the block of inner shell lines */
	for( i=0; i < nUTA; i++ )
	{
		if( UTALines[i].Aul > 0. )
		{
			/* populations have not been set */
			UTALines[i].PopOpc = dense.xIonDense[UTALines[i].nelem-1][UTALines[i].IonStg-1];
			UTALines[i].PopLo = dense.xIonDense[UTALines[i].nelem-1][UTALines[i].IonStg-1];
			UTALines[i].PopHi = 0.;
			RT_line_one_tauinc(&UTALines[i], -4 , -4 , -4 , i );
		}
	}

	/* all hyper fine structure lines  */
	for( i=0; i < nHFLines; i++ )
	{
		/* remember current gas-phase abundances */
		float save = dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1];

		/* bail if no abundance */
		if( save<=0. ) continue;

		/* set gas-phase abundance to total times isotope ratio */
		dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1] *= hyperfine.HFLabundance[i] ;

		RT_line_one_tauinc(&HFLines[i] , -5 , -5 , -5 , i );

		/* put the correct gas-phase abundance back in the array */
		dense.xIonDense[HFLines[i].nelem-1][HFLines[i].IonStg-1] = save;
	}
#	if 0
	{
		/* this is an option to print out one of the two photon continua */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{	
			fprintf(ioQQQ,"isotope pops\t%.3e\t%.3e\t%.3e\t%.3e\n",
				/*HFLines[0].PopOpc , HFLines[0].opacity,
				TauLines[ipH21cm].PopOpc , TauLines[ipH21cm].opacity);*/
				HFLines[0].TauIn , HFLines[0].TauTot,
				TauLines[ipH21cm].TauIn , TauLines[ipH21cm].TauTot);
		}
	}
#	endif

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

	/* do large FeII atom if this is enabled */
	FeIITauInc();

	/* increment optical depth for the H2 molecule */
	H2_RT_tau_inc();

	/* following is for static atmosphere */
	if( wind.windv == 0. )
	{
		/* iron fe feii fe2 
		 * first are overlapping feii lines */
		for( i=0; i < NFEII; i++ )
		{
			/* optical depths for Feii dest of lya when large feii not used */
			fe2ovr_la.Fe2TauLte[i] += fe2ovr_la.feopc[i]*(float)radius.drad_x_fillfac;
		}
	}

	if( trace.lgTrace && trace.lgOptcBug )
	{
		fprintf( ioQQQ, " RT_tau_inc updated optical depths:\n" );
		prtmet();
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " RT_tau_inc returns.\n" );
	}

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

