/* 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_tauinc increment optical depths for all heavy element lines, zone by zone,
 * mainly called by RT_tau_inc, but also by FeII */
#include "cddefines.h"
#include "doppvel.h"
#include "geometry.h"
#include "rfield.h"
#include "radius.h"
#include "wind.h"
#include "rt.h"

void RT_line_one_tauinc(EmLine * t ,
	/* following four are flags to generate info if some species has extreme maser */
	long int maser_flag_species,
	long int maser_flag_ion,
	long int maser_flag_hi,
	long int maser_flag_lo )

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

	/* routine increments optical depths for static or exanding atmosphere */

	/* use static solution for both d-critical and static, for d-critical speeds are only
	 * roughly sonic and thermal for H itself, Lya is dominant, and flow across HII region
	 * has nearly constant speed.  so this seems like better bet.
	 * for d-critical velocity gradients are too small to use a purely local calulation
	 * of the total optical depth - the depth would be overestimated */
	/*TODO	1	dynamics; this test assumes that neg vel are subsonic, so that sobolev length
	 * would overestimate the optical depth, since ion is at most present over computed
	 * slab, and possibly more.  */
	if( wind.windv <= 0. )
	{
		/* static and negative velocity solutions */
		float dTau_total;
		float opacity = t->opacity / DoppVel.doppler[ t->nelem -1] ;
		/* this is optical depth within this one line itself */
		float dTau_line = (float)( t->PopOpc * opacity * radius.drad_x_fillfac );

		/* if fine opacities fine grid fine mesh are calculated then use them to get optical depth 
		 * to continuum source */
		if( t->ipFine>0 && t->ipFine<rfield.nfine && rfield.lgOpacityFine )
		{
			/* pick up fine opacity at center of line */
			dTau_total = rfield.fine_opac[t->ipFine] * (float)radius.drad_x_fillfac;

			/* sanity check - total fine opacity should be greater than single ine opacity */
			ASSERT( fabs(dTau_line) < 0.01 || 
				rfield.fine_opac[t->ipFine]*1.1 >= opacity*t->PopOpc  );
		}
		else
		{
			dTau_total = dTau_line;
		}
		/*if( rfield.lgOpacityFine )
			fprintf(ioQQQ,"DEBUG lgOpacityFine oops\n");*/
		/* >>chng 03 aug 20, only use self abs for this optical depth */
		t->TauIn += dTau_line;

		/* >>chng 03 feb 12, add correction for lines that are not this line */
		/* this is optical depth to continuum source, includes line overlap,
		 * is used for continuum pumping */
		t->TauCon += dTau_total;

		/* keep track of any masers */
		if( dTau_line < rt.dTauMase )
		{
			rt.dTauMase = dTau_line;
			rt.mas_species = maser_flag_species;
			rt.mas_ion = maser_flag_ion;
			rt.mas_hi = maser_flag_hi;
			rt.mas_lo = maser_flag_lo;
			/*if( dTau_line < -0.1 )
				fprintf(ioQQQ,
				" maser hittt, iteraton %li zone %li species %li ion %li hi, lo %li %li val %.2e\n",
				iteration, nzone, 
				rt.maser_flag_species , rt.maser_flag_ion , rt.maser_flag_hi , rt.maser_flag_lo , dTau_line );*/
		}

		/* this will cause comment to be printed in prtComments */
		if( rt.dTauMase < -1. )
		{
			rt.lgMaserCapHit = TRUE;
		}
	}

	else
	{
		double absAccelTot = fabs(wind.AccelTot);
		float oldtau = t->TauIn;
		double hold;
		/* wind solution */

		/* this is effective length scale for sobolev or lvg approximation.
		 * This is equation 3 of
		 * >>refer	RT	wind	Castor, J.I., Abbott, D.C., & Klein, R.I., 1975, ApJ, 195, 157
		 *
		 * >>chng 01 mar 24, change to this formalism from previous, constant velocity
		 * approximation 
		 */
		/*t->TauIn = wind.wndeff * t->PopOpc * t->opacity;*/
		/* do not let the physical scale exceed the current depth */
		hold = SDIV(absAccelTot);
		t->TauIn = (float)(MIN2( radius.depth, DoppVel.doppler[ t->nelem -1]/ hold ) *
			geometry.FillFac * t->PopOpc * t->opacity);
		if( t->TauIn < 0. && t->TauIn < oldtau )
		{
			/* >>chng 04 dec 26,
			 * we have a maser, dont let it run away - it can do this in this case wind no
			 * memory of optical depth scale up to here - let it change by, at most,
			 * -1 from previous optical depth */
			t->TauIn = MAX2( oldtau-1.f, t->TauIn) ;
		}

		t->TauCon = t->TauIn;

		t->TauTot = t->TauIn;
	}

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