/* 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_line_one_tau_reset computes average of old and new optical depths for new scale at end of iter,
 * called by update, also FeIILevelPops */
#include "cddefines.h"
#include "opacity.h"
#include "geometry.h"
#include "lines_service.h"
#include "prt.h"
#include "rt.h"

void RT_line_one_tau_reset(EmLine * t, 
  double WeightNew)
{
	char chSave[11];

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

	/* option to print masing lines, set with print maser */
	if( prt.lgPrtMaser && 
		( t->TauTot < -0.01 || t->TauIn <  -0.01) )
	{
		strcpy( chSave, chLineLbl(t) );
		fprintf( ioQQQ, " Masing line:%10.10s t(in, out)=%10.2e%10.2e\n", 
		  chSave, t->TauIn, t->TauTot );
	}

	/* iteration is 1 when this routine is called at end of first iteration
	 * lgStatic is false by default, set true with command sphere SPHERE STATIC */
	if( iteration == 1 && (!geometry.lgStatic) )
	{
		/* end of first iteration, first estimate of total optical depth 
		 * is now the inward optical depth - set T(ipLnTauTot) to T(1)
		 * DoubleTau normally 1, set to 2 by DoubleTau command in order
		 * to simulate two-sided photoionization */
		t->TauTot = t->TauIn*rt.DoubleTau;
		t->TauIn = MIN2(opac.taumin,t->TauTot/2.f);
	}
	else if( geometry.lgSphere && geometry.lgStatic )
	{
		/* static sphere, both sides interact */
		if( iteration == 1 )
		{
			/* on first iteration TauIn was 0 for first zone */
			t->TauTot = 2.f*t->TauIn;
		}
		else
		{
			t->TauTot = (float)(t->TauIn*WeightNew + t->TauTot*
			  (1. - WeightNew));
		}
		t->TauIn = t->TauTot/2.f;
	}
	else
	{
		/* take some sort of mean of old and new limiting optical depths */
		if( t->TauIn > 0. && t->TauTot > 0. )
		{
			/* this is case where no maser occurs */
#			if 0
			t->TauTot = (float)pow(10.,log10(t->TauTot)*(1. - WeightNew) + 
				log10(t->TauIn*rt.DoubleTau)*WeightNew);
#			endif
			/* >>chng 03 nov 06, try taking linear average */
			t->TauTot = (float)(t->TauTot*(1. - WeightNew) + 
				t->TauIn*rt.DoubleTau*WeightNew);
			t->TauIn = MIN2(opac.taumin,t->TauTot/ 2.f);
		}
		else
		{
			/* case where maser did occr */
			t->TauTot = (t->TauIn*rt.DoubleTau + t->TauTot)/2.f;
			t->TauIn = MIN2(opac.taumin,t->TauTot/2.f);
		}
	}

	/* this is escape prob */
	t->Pesc = 0.5f*(1.f + 1.f/MAX2(1.f,t->TauTot));

	/* this is fraction inward */
	t->FracInwd = MIN2(1.f,1.5f-t->Pesc);

	/* this is destruction probability
	 * >>chng 96 sep 4, was not here, needed to add since now taking
	 * mean of old and new dest prob */
	t->Pdest = 0.;
	t->Pelec_esc = 0.;

	/* optical depth to the continuum source */
	t->TauCon = opac.taumin;

	/* >>chng 01 sep 01, zero out some pops and energies */
	t->PopLo = 0.;
	/* >>chng 97 jul 21, added following zero
	 * population of upper level */
	t->PopHi = 0.;
	/* population of lower level with correction for stim emission */
	t->PopOpc = 0.;
	/* following two heat exchange excitation, deexcitation */
	t->cool = 0.;
	t->heat = 0.;
	/* intensity of line */
	t->xIntensity = 0.;
	/* number of photons emitted in line */
	t->phots = 0.;

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

