/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*RT_stark compute stark broadening escape probabilities using Puetter formalism */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "dense.h"
#include "hydrogenic.h"
#include "phycon.h"
#include "rt.h"

void RT_stark(void)
{
	long int ipLo, 
	  ipHi;
	double aa , ah, 
	  stark, 
	  strkla;

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

	if( !rt.lgStarkON || dense.eden < 1e8 )
	{
		for( ipHi=0; ipHi < iso.numLevels[ipH_LIKE][ipHYDROGEN]; ipHi++ )
		{
			for( ipLo=0; ipLo < iso.numLevels[ipH_LIKE][ipHYDROGEN]; ipLo++ )
			{
				hydro.pestrk[ipHi][ipLo] = 0.;
				hydro.pestrk[ipLo][ipHi] = 0.;
			}
		}
	
#		ifdef DEBUG_FUN
		fputs( " <->RT_stark()\n", debug_fp );
#		endif
		return;
	}

	/* evaluate Stark escape probability from 
	 * >>ref Puetter Ap.J. 251, 446. */

	/* coefficients for Stark broadening escape probability
	 * needs factor of nu/nl ** 3 *xl to be Puetters AH */
	ah = 6.9e-6*1000./1e12/(phycon.sqrte*phycon.te10*phycon.te10*
	  phycon.te03*phycon.te01*phycon.te01)*dense.eden;

	/* coefficient for all lines except Ly alpha */
	stark = 0.264*pow(ah,0.4);

	/* coefficient for Ly alpha */
	strkla = 0.538*ah*4.*9.875*(phycon.sqrte/phycon.te10/phycon.te03);

	/* Lyman lines always have outer optical depths */
	/*ASSERT( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn> 0. );*/
	/* >>chng 02 mar 31, put in max, crashed on some first iteration 
	 * with negative optical depths,
	 * NB did not understand why neg optical depths started */
	aa = pow(SDIV(EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn),-0.75f);
	hydro.pestrk[ipH1s][ipH2p] = 
		(float)(strkla/2.*MAX2(1.,aa));

	/*TODO	2	- Stark is disabled for now since Lya escape causes density dependent
	 * feedback on the radiative transfer.  Would need to redo the escape
	 * probs every time the electron density is updated - see blr89.in for an 
	 * example */
	hydro.pestrk[ipH1s][ipH2p] = (float)MIN2(.01,hydro.pestrk[ipH1s][ipH2p]);
	hydro.pestrk[ipH1s][ipH2p] = 0.;
	hydro.pestrk[ipH2p][ipH1s] = 
		hydro.pestrk[ipH1s][ipH2p]*EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Aul;

	for( ipHi=3; ipHi < iso.numLevels[ipH_LIKE][ipHYDROGEN]; ipHi++ )
	{
		hydro.pestrk[0][ipHi] = (float)(stark*hydro.strkar[0][ipHi]/2.*pow(MAX2(1.,
		  EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipH1s].TauIn),-0.75));

		hydro.pestrk[0][ipHi] = (float)MIN2(.01,hydro.pestrk[0][ipHi]);
		hydro.pestrk[ipHi][0] = EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipH1s].Aul*
		  hydro.pestrk[0][ipHi];
	}

	/* all other lines */
	for( ipLo=ipH2s; ipLo < (iso.numLevels[ipH_LIKE][ipHYDROGEN] - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < iso.numLevels[ipH_LIKE][ipHYDROGEN]; ipHi++ )
		{
			aa = stark*hydro.strkar[ipLo][ipHi]*
			  pow(MAX2(1.,EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].TauIn),-0.75);
			hydro.pestrk[ipLo][ipHi] = 
				(float)MIN2(.01,aa);

			hydro.pestrk[ipHi][ipLo] = EmisLines[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Aul*
			  hydro.pestrk[ipLo][ipHi];
		}
	}

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

