/* 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_static do line radiative transfer for static geometry,
 * evaluates escape and destruction probability, 
 * called by HydroPEsc, RT_line_all  */
/*RT_line_wind do line radiative transfer for wind geometry */
#include "cddefines.h"
#include "rfield.h"
#include "doppvel.h"
#include "dense.h"
#include "opacity.h"
#include "lines_service.h"
#include "conv.h"
#include "radius.h"
#include "rt.h"

/*RT_line_static do line radiative transfer for static geometry,
 * called by RT_line_all */
static void RT_line_static(
	/* the em line we will work on  */
	EmLine * t , 
	/* when this is true do both esc and dest probs,
	 * when falst, only dest probs */
	int lgDoEsc ,
	/* this is option to not include line self shielding across this zone.
	 * this can cause pump to depend on zone thickness, and leads to unstable
	 * feedback in some models with the large H2 molecule, due to Solomon
	 * process depending on zone thickness and level populations. */
	int lgShield_this_zone );

/*RT_line_one do rt for emission line structure - calls RT_line_static or RT_line_wind */
void RT_line_one(
	/* the em line we will work on  */
	EmLine * t , 
	/* when this is true do both esc and dest probs,
	 * when falst, only dest probs */
	int lgDoEsc ,
	/* if true then we want to update the fine opacity scale */
	int lgDoFine_opac_update ,
	/* this is option to not include line self shielding across this zone.
	 * this can cause pump to depend on zone thickness, and leads to unstable
	 * feedback in some models with the large H2 molecule, due to Solomon
	 * process depending on zone thickness and level populations. */
	int lgShield_this_zone )
{

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

	/* skip line transfer if requested with 'no line transfer' command, but never skip Lya */
	if( !rfield.lgDoLineTrans )
	{
		if( /*&& conv.nPres2Ioniz */ (t->iRedisFun != ipLY_A) )
			return;
	}

	/* always call single routine */
	RT_line_static(
		/* the line we will work on  */
		t , 
		/* when this is true do both esc and dest probs,
			* when falst, only dest probs */
		lgDoEsc ,
		lgShield_this_zone);

	/* >>chng 03 feb 12, add this array */
	/* define fine opacity fine grid fine mesh */
	/* rfield.lgOpacityFine flag set flast with no fine opacities command */
	if( lgDoFine_opac_update && t->ipFine >= 0 && t->PopOpc > 0. &&
		t->ipFine<rfield.nfine && rfield.lgOpacityFine )
	{
		long int i , ip_lo , ip_hi;

		/* number of fine opacity cells corresponding to one doppler width for current
		 * temperature, velocity field, and nuclear mass,
		 * rfield.fine_opac_velocity_width is width per cell, cm/s */
		float cells_wide_1x = DoppVel.doppler[t->nelem-1]/rfield.fine_opac_velocity_width;

		/* line center opacity - want float since will add to fine opacity array,
		 * which is float */
		float opac_line =  (float)t->PopOpc * t->opacity / DoppVel.doppler[t->nelem-1];

		/* this is effective optical depth to this point. Do not work with line if this product
		 * is substantially less than unity */
		double dTauEffec = opac_line*radius.depth_x_fillfac;

		/* >>chng 04 mar 29, use dTauEffec to choose whether line is important */
		/*if( t->TauIn < 0.0001 )*/
		/*if( dTauEffec > 0.1 )*/
		/* >>chng 05 feb 27, keep all opacities */
		if( dTauEffec > SMALLFLOAT )
		{
			/* >>chng 04 mar 29, rewrite to use symmetry between 
			 * two halves of line profile */
			/* width of optically thick line, do 4x with exponential core,
			 * must be at least one cell */
			long int nCells_core = (long)(cells_wide_1x*4.+1.5f);
			float opacity;

			/* core is symmetric - make sure both upper and lower bounds
			 * are within continuum energy bin */
			if( t->ipFine - nCells_core < 1 )
				nCells_core = t->ipFine - 1;
			if( t->ipFine + nCells_core > rfield.nfine )
				nCells_core = t->ipFine -rfield.nfine - 1;

			/* want core to be at least one cell wide */
			nCells_core = MAX2( 1 , nCells_core );

			/* line center itself, must not double count here */
			rfield.fine_opac[t->ipFine] += opac_line;
			for( i=1; i<nCells_core; ++i )
			{
				/* square of distance from line center in units of doppler width */
				double xsq = POW2(i/cells_wide_1x);

				/* the line opacity at that point, includes doppler core and damping wings */
				opacity = opac_line*(float)(
					exp( -xsq ) + t->damp/MAX2(1., (1.772f*xsq)) );

				rfield.fine_opac[t->ipFine+i] += opacity;
				rfield.fine_opac[t->ipFine-i] += opacity;
			}

			/* include damping wings if optical depth is large */
			/* >>chng 04 mar 29, rewrite to use symmetry in two halves of line */
			if( dTauEffec*t->damp/9. > 0.1 )
			{
				/* do damping wings if very optically thick
				 * this is number of cells to extend the damping wings - cells_wide is one dop width 
				 * 56.419 is 1 / (0.01 * sqrt pi) ) */
				/* tests with th85orion, stopping at half -h2 point, showed that 0.01 was
				 * needed for outer edge, given the defin of dTauEffec */
				float x = (float)sqrt( dTauEffec * t->damp *56.419 ) * cells_wide_1x;
				long int nCells_damp;
				/* >>chng 04 mar 24, add test on size of x, which can exceed
				* limits of long in extreme optical depths */
				if( x<LONG_MAX )
				{
					nCells_damp = (long)x;

					if( t->ipFine-nCells_damp < 1 )
						nCells_damp = t->ipFine-1 ;

					if( t->ipFine+nCells_damp > rfield.nfine )
						nCells_damp = rfield.nfine - t->ipFine-1;
				}
				else
				{
					/* x was too big, just set to extreme range, which is
					 * number of cells to nearest boundary */
					nCells_damp = MIN2( rfield.nfine-t->ipFine , t->ipFine )-1;
				}

				ASSERT( nCells_core>0 );
				for( i=nCells_core; i<nCells_damp; ++i )
				{
					/* distance from line center in doppler units */
					float xsq = POW2(i/cells_wide_1x);
					/* term in denominator is x^2 */
					opacity = opac_line*t->damp/(1.772f*xsq );
					rfield.fine_opac[t->ipFine+i] += opacity;
					rfield.fine_opac[t->ipFine-i] += opacity;
				}
			}
		}

		else if( dTauEffec> 0.001 )
		{
			/* optically thin line - use square topped line profile */
			ip_lo = t->ipFine - (long)(cells_wide_1x+0.5f);
			ip_hi = t->ipFine + (long)(cells_wide_1x+0.5f);
			ip_lo = MAX2(1, ip_lo );
			/* make sure at least one cell wide */
			ip_hi = MAX2( ip_hi , t->ipFine+1);
			ip_hi = MIN2( rfield.nfine , ip_hi );
			for( i=ip_lo; i<ip_hi; ++i )
			{
				rfield.fine_opac[i] += opac_line;
			}
		}
	}

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

	return;
}

/*RT_line_static do line radiative transfer for static geometry */
static void RT_line_static(
	/* the em line we will work on  */
	EmLine * t , 
	/* when this is true do both esc and dest probs,
	 * when false, only dest probs */
	int lgDoEsc ,
	/* this is option to not include line self shielding across this zone.
	 * this can cause pump to depend on zone thickness, and leads to unstable
	 * feedback in some models with the large H2 molecule, due to Solomon
	 * process depending on zone thickness and level populations. */
	int lgShield_this_zone )
{
	int lgGoodTau;
	long int nelem;
	int nRedis = -1;

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

	/* this is the routine that computes much of the radiative transfer
	 * physics for lines for static case.  */

	/*ASSERT( t->IonStg < LIMELM+4 );*/
	/*ASSERT( t->nelem-1 < LIMELM+4 );*/

	/* >>chng 01 aug 18, do not evaluate if no population */
	/* molecules are evaluated here, these must be valid */
	/* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
	 * that rates are evaluated at start of a new iteration */
	if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0.  && !lgDoEsc )
	{
		/* zero population, return after setting everything with side effects */
		t->Pesc = 1.f;

		/* inward escaping fraction */
		t->FracInwd = 0.5;

		/* pumping rate */
		t->pump = 0.;

		/* destruction probability */
		t->Pdest = 0.;
		t->Pelec_esc = 0.;

		/* damping constant */
		t->damp = 0.;

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

		return;
	}
	/* check that we don't have a runaway maser */
	else if( t->TauIn < -30. )
	{
		fprintf( ioQQQ, "PROBLEM RT_line_static called with large negative optical depth, zone %.2f, setting lgAbort true.\n",
			fnzone );
		DumpLine(t);
		/* >>>chng 00 apr 23, return busted true instead of exit here, so that code will
		 * not hang under MPI */
		lgAbort = TRUE;

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

	/* this checks if we have overrun the optical depth scale
	 * possible for line to mase, reason for sec test
	 * >>chng 97 jan 08, on iterations where the first iteration had NO species,
	 * but it does exist on later iterations, added test for existence of ion */
	/* >>chng 03 may 14, go to function visible to whole code */
	/* test on lgTauOutOn is to say tau ok if this is first iteration 
	 * and outward optical depths turned off */
	if( !opac.lgTauOutOn || lgTauGood( t ) )
	{
		lgGoodTau = TRUE;
	}
	else
	{
		/* do not do anything if we have overrun the scale, */
		lgGoodTau = FALSE;
	}

	/* atomic number of this element on the C scale*/
	nelem = t->nelem-1;

	/* damping constant for the line, save it */
	t->damp = t->dampXvel / DoppVel.doppler[t->nelem-1];
	/*if( t->damp <0. )
	{
		fprintf(ioQQQ,"DEBUG hit it \n");
	}*/
	ASSERT( t->damp >= 0. );

	/* static solution - which type of line will determine
	 * which redistribution function */
	/* iRedisFun == 1 - alpha resonance line, partial redistribution,
	 * ipPRD == 1 */
	if( t->iRedisFun == ipPRD )
	{
		/* incomplete redistribution with wings */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = esc_PRD( t->TauIn , t->TauTot , t->damp );

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}

#			if 0
			/* update pumping even if over run optical depth scale */
			if( rfield.lgInducProcess )
			{
				/* the inward-looking escape probability, TauCon is optical depth to
				 * the continuum source, including all overlapping lines*/
				rt.wayin = (float)esc_PRD_1side(t->TauCon,t->damp);
			}
#			endif
		}
		nRedis = ipDEST_INCOM;
	}

	/* complete redistribution without wings - t->ipLnRedis is ipCRD == -1 */
	else if( t->iRedisFun == ipCRD )
	{
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				/* >>chng 01 mar -6, escsub will call any of several esc prob routines,
				 * depending of how core is set.  We always want core-only for this option,
				 * so call  esca0k2(tau) directly */
				t->Pesc = esc_CRDcore( t->TauIn , t->TauTot );

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}

#			if 0
			/* continuum pumping rate
			 * the "no induced" command causes ContPump to set 0 */
			if( rfield.lgInducProcess )
			{
				/* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp )
				 * this is to make it more like C90 was
				 * complete redistribution with doppler core only */
				rt.wayin = (float)esca0k2(t->TauCon);
			}
#			endif
		}
		nRedis = ipDEST_K2;
	}

	/* chng 01 feb 27, add CRD with wings, = 2 */
	else if( t->iRedisFun == ipCRDW )
	{
		/* complete redistribution with damping wings */
		if( lgDoEsc )
		{
			if( lgGoodTau )
			{
				t->Pesc = esc_CRDwing( t->TauIn , t->TauTot , t->damp );

				/* inward escaping fraction */
				t->FracInwd = rt.fracin;
			}
#			if 0
			/* continuum pumping rate
			 * the "no induced" command causes ContPump to set 0 */
			if( rfield.lgInducProcess )
			{
				/* wayin = esc_CRDwing_1side( t(ipLnTauCon),damp )
				 * this is to make it more like C90 was
				 * complete redistribution with doppler core only */
				/* >>chng 01 may 29, had been esc_PRD_1side, as per above comment,
				 * but this should include damping wings.  */
				/*rt.wayin = (float)esc_PRD_1side(t->TauCon,damp);*/
				rt.wayin = (float)esc_CRDwing_1side(t->TauCon,t->damp);
			}
#			endif
		}

		nRedis = ipDEST_K2;
	}

	/* chng 03 may 14, special Lya case*/
	else if( t->iRedisFun == ipLY_A )
	{
		/* incomplete redistribution with wings, for special case of Lya
		 * uses fits to Hummer & Kunasz numerical results 
		 * this routine is different because escape and dest probs
		 * are evaluated together, so no test of lgDoEsc */
		if( lgGoodTau )
		{
			double dest , esin;

			/* this will always evaluate escape prob, no matter what lgDoEsc is.
			* The destruction prob comes back as dest */
			t->Pesc = RTesc_lya(&esin , &dest , t->PopOpc , nelem);

			/* this is current destruction rate */
			t->Pdest = dest;

			/*  this is fraction of line which is inward escaping */
			t->FracInwd = rt.fracin;
		}
#		if 0
		/* update pumping even if over run optical depth scale */
		if( rfield.lgInducProcess )
		{
			/* the inward-looking escape probability */
			rt.wayin = (float)esc_PRD_1side(t->TauCon,t->damp);
		}
#		endif
	}
	else
	{
		fprintf( ioQQQ, " RT_line_static called with redistribution function zero\n" );
		ShowMe();
		puts( "[Stop in RT_line_static]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* pumping by incident and diffuse continua */
	/* >>chng 03 may 15, this block had been repeated all four times above, move down
	 * here to use common code, only rt.wayin is defined above */
	/* >>chng 03 may 17, do not evaluate when drad is zero, before first
	 * zone thickness is established */
	if( (lgDoEsc  || (t->iRedisFun == ipLY_A) ) && rfield.lgInducProcess )
	{
		float dTau;
		double pumpsave = t->pump;
		float shield_continuum;

		/* >>chng 04 apr 18, break out line continuum shielding into this function */
		shield_continuum = RT_continuum_shield_fcn( t );

		/* continuum upward pumping rate, A gu/gl abs prob occnum
		 * the "no induced" command causes ContPump to set 0 */
		/*t->pump = t->Aul * t->gHi / t->gLo * rt.wayin **/
		/* >>chng 05 feb 16, add outward diffuse emission */
		/*t->pump = t->Aul * t->gHi / t->gLo * shield_continuum *
			rfield.OccNumbIncidCont[t->ipCont-1];*/
		t->pump = t->Aul * t->gHi / t->gLo * shield_continuum *(
			rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1] );

		/* >>chng 99 dec 02, add correction for line optical depth across zone */
		/* correction for attenuation within line across zone */
		/* >>chng 99 mar 18, from 0 to 1e-8 in following to avoid underflow, 
		 * this caused an overestimate of pumping rates */
		/* >>chng 00 dec 19, nova.in oscillated since correction factor was
		 * very large across zones that were already very optically thick, dr changed,
		 * causing feedback between ionization and thickness
		 * the intention of this correction was to only apply when on the 
		 * verge of optically thick,
		 * change lower bound on product t->dTau * radius.drad_x_fillfac to 1e-3 to avoid
		 * constantly getting unity, also check on relative change optical depth/optical depth*/
		/* this can drive oscillations in very neutral gas, as blr92.in, blr89.in
		 * or nova.in, when higher lyman lines have produce 
		 * t->dTau * radius.drad_x_fillfac about unity, where small changes in
		 * pump can change ionization, dTau, and hence pump */
		/*if( t->dTau * radius.drad_x_fillfac > 1e-8 )*/
		/* >>chng 02 jun 14, nova.in again - upper limit had been 10 - but
		 * there were major jumps in heating rate at this point - increase 10 -> 1e10
		 * not clear why there should be an upper limit at all, other than to save time? */
		/* update dTau since affects pump rate */
		/* >>chng 03 jun 18, include option to not shield across this zone */
		/* >>chng 03 jul 19, add continuum opacity */
		dTau =  (float)((t->PopOpc * t->opacity / DoppVel.doppler[nelem] + opac.opacity_abs[t->ipCont-1])*
			 radius.drad_x_fillfac_mean);
		/* >>chng 04 mar 04, do not want to include maser action across this zone */
		dTau = MIN2(1.f, dTau);
		/*if( 0 && lgShield_this_zone && dTau * radius.drad_x_fillfac > 1e-3 && dTau * radius.drad_x_fillfac < 1e10 )*/
		/* >>chng 04 mar 04, use mean of past few zone thicknesses rather than
		 * current zone thickness since this can induce wild oscillations, see
		 * nova.in test case */
		if( lgShield_this_zone && dTau > 1e-3 && dTau < 1e10 )
		{
			/* during very first search dTau is -1e38*/
			/* >>chng 04 mar 04, use mean, see above comment */
			/*t->pump *= log(1. + dTau * radius.drad_x_fillfac ) / (dTau * radius.drad_x_fillfac);*/
			t->pump *= log(1. + dTau ) / dTau;
		}

		/* >>chng 00 Oct 03, add pumping by local diffuse continuum,
		 * rfield.DiffPumpOn is unity unless process disabled by setting to 0 in  parsedont.c
		 * with no DIFFuse LINE PUMPing command */
		if( dTau*rfield.DiffPumpOn > SMALLFLOAT )
		{
			/* >>chng 02 mar 14, added brackets around scale*rfield.OccNumbDiffCont so
			 * that smallest and largest number are multiplied first to prevent overflow, PvH */
			double scale;
			if( opac.lgTauOutOn )
			{
				if( lgGoodTau )
				{
					/* >>chng 03 jul 19, add check for tau = 1 in continuum */
					/* scale for tau = 1 in the continuum */
					double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
					/* inward looking depth, dTau is opacity in line, cm^-1,
					* so multiplying by this scale is same as dividing by opacity in 
					* source function */
					double scale_in =  MIN3( 1./ dTau , radius.depth    , scale_con);
					double scale_out = MIN3( 1./ dTau , radius.Depth2Go , scale_con );
					t->pump += t->Aul * t->gHi / t->gLo * 0.5 *((scale_in+scale_out)*
						rfield.OccNumbDiffCont[t->ipCont-1]);
				}
			}
			else 
			{
				/* >>chng 03 jul 19, add check for tau = 1 in continuum */
				/* scale for tau = 1 in the continuum */
				double scale_con = 1. / SDIV(opac.opacity_abs[t->ipCont-1]+opac.opacity_sct[t->ipCont-1]);
				/* first iteration, only inward looking depths really known */
				scale = 1./ dTau;
				scale = MIN3( scale , radius.depth , scale_con );
				t->pump += t->Aul * t->gHi / t->gLo * (scale*
					rfield.OccNumbDiffCont[t->ipCont-1]);
			}
		}
		/* this stop oscillations from developing in very deep blr models 
		 * like blr89 and blr98, but must not be applied too early since
		 * line optical depths change greatly in first few zones */
		if( nzone>50 )
			t->pump = (pumpsave + t->pump) / 2.;
	}

	/* option to kill induced processes */
	else if( !rfield.lgInducProcess )
	{
		t->pump = 0.;
	}

	/* escape following scattering off an electron */
	/* >>chng 02 may 08, put elec scat here, had been only for h-like species,
	 * in hydropesc */
	/* this is turned off with the no scattering escape command */
	if( rt.lgElecScatEscape )
	{
		/* in the EmLine structure PopOpc is the pop relative to the ion,
		 * but has been converted to a physical population before this routine
		 * was called.  No need to convert here */
		double opac_line = t->PopOpc * t->opacity/DoppVel.doppler[nelem];

		if( opac_line > SMALLFLOAT )
		{
			/* the old escape probability */
			double eesc_save = t->Pelec_esc;
			/* the opacity in electron scattering */
			double es = dense.eden*6.65e-25;
			/* >>chng 04 may 12, use total scattering
			double es = opac.opacity_sct[t->ipCont-1]; */
			/* the optical depth in electron scattering*/
			double taues = es * radius.depth; 
			/* this is equation 5 of 
			 *>>refer	line	desp	Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/
			/* >>chng 03 may 07, add correction for total optical depth to electrons */
			double opacity_ratio = es/(es+opac_line) * (taues/(1.+taues));
			/* >>chng 04 may 12, just do line to total opac 
			double opacity_ratio = es/(es+opac_line);*/
			t->Pelec_esc = opacity_ratio * MAX2(0.,1.-t->Pesc-t->Pdest);
			/*KLUDGE >>chng 03 feb 02 
			 * take mean of old and new if electron scattering escape */
			if( !conv.lgSearch )
				t->Pelec_esc = t->Pelec_esc*0.75 + eesc_save*0.25;  
		}
		else
			t->Pelec_esc = 0.;
	}

	/* >>>chng 99 dec 18, did not have the test for lgGoodTau, and so clobbered
	 * dest prob when overrun */
	/* only do this if not Lya special case, since dest prob already done */
	if( lgGoodTau && t->iRedisFun!=ipLY_A)
	{
		double DestSave = t->Pdest;
		t->Pdest = RT_DestProb(
			/* the abundance of the species */
			t->PopOpc,
			/* line center opacity in funny units (needs a vel) */
			t->opacity,
			/* array index for line in continuum array,
			 * on f not c scale */
			t->ipCont,
			/* line width in velocity units */
			DoppVel.doppler[nelem],
			/* escape probability */
			t->Pesc,
			/* redistribution function */
			nRedis);

		/* >>chng 03 may 10, this had been 3/4 new + 1/4 old, but had to be combined
		 ^ with half and half when this routine was called, to converge very high Z high U
		 * quasar models.  bring all these factors together here */
		/*t->Pdest = 0.75f*t->Pdest + 0.25f * DestSave;*/
#		define OLDFAC	(0.625)
		if( nzone>1 )
			t->Pdest = (1.-OLDFAC)*t->Pdest + OLDFAC * DestSave;
	}

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

	return;
}

#if 0
/*RT_line_wind do line radiative transfer for wind geometry */
static void RT_line_wind(EmLine * t , int lgDoEsc);
/*RT_line_wind do line radiative transfer for wind geometry */
static void RT_line_wind(EmLine * t , int lgDoEsc)
{
	float velocity;

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

	/* 
	 * this is the routine that computes much of the radiative transfer
	 * physics for lines in wind case.  Lines diveded into two types,
	 * high quality (level1) and lower quality (atom_level2) g-bar lines
	 */

	/* >>chng 01 aug 18, do not evaluation if no population */
	/* >>chng 01 oct 25, add test for lgDoEsc, this is to make sure
	 * that rates are evaluated at start of a new iteration */
	if( dense.xIonDense[ t->nelem-1][t->IonStg-1] == 0.  && !lgDoEsc )
	{
		/* zero population, return after setting everything with side effects */
		t->Pesc = 1.f;

		/* inward escaping fraction */
		t->FracInwd = 0.5;

		/* pumping rate */
		t->pump = 0.;

		/* destruction probability */
		t->Pdest = 0.;
		t->Pelec_esc = 0.;

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

	/*confirm pops and cs ok */
	ASSERT( t->ColOvTot <= 1. && t->ColOvTot >= 0. );

	/* t->nelem is atomic number of species */
	ASSERT( t->nelem >= 1 && t->ipCont > 0 );

	/* line width, thermal + turbulence */
	velocity = DoppVel.doppler[t->nelem-1];

	t->damp = 0.;

	/* use only one-sided escape probabilities
	 * get escape probability, this gets fractions too */
	if( lgDoEsc )
	{
		t->Pesc = esc_CRDwing_1side(t->TauIn,t->damp);

		/* inward fraction */
		t->FracInwd = 0.5;

		/* continuum pumping */
		if( rfield.lgInducProcess )
		{
			/* >>chng 05 feb 16, include diffuse outward emission */
			/*t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*
			  rfield.OccNumbIncidCont[t->ipCont-1];*/
			t->pump = t->Aul*t->gHi/t->gLo*t->Pesc*(
			  rfield.OccNumbIncidCont[t->ipCont-1] + rfield.OccNumbContEmitOut[t->ipCont-1]);
		}
		else
		{
			t->pump = 0.;
		}
	}

	/* get destruction probability */
	t->Pdest = 
		RT_DestProb(
		t->PopOpc,
	    /* its line absorption cross section */
		t->opacity,
		/* index for line in continuum array, on f not c scale */
		t->ipCont,
	    velocity,
		t->Pesc, 
		ipDEST_K2);
	/*TODO	2	- consider using ipDEST_INCOM instead of K2, it is much faster */

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

	return;
}
#endif


