/*strk compute stark broadening escape probabilities using Puetter formalism */
#include "cddefines.h"
#include "taulines.h"
#include "cdminmax.h"/* side effects, variable args */
#include "nhe1lvl.h"
#include "tepowers.h"
#include "strkon.h"
#include "he1as.h"
#include "hydrogenic.h"
#include "he1str.h"
#include "he1tau.h"
#include "he1nxt.h"
#include "phycon.h"
#include "strk.h"

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

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

	if( !strkon.lgStarkON || phycon.eden < 1e8 )
	{
		for( ipHi=0; ipHi <= nhlevel; ipHi++ )
		{
			for( ipLo=0; ipLo <= nhlevel; ipLo++ )
			{
				hydro.pestrk[ipHi][ipLo] = 0.;
				hydro.pestrk[ipLo][ipHi] = 0.;
			}
		}
		for( ipHi=0; ipHi < NHE1LVL; ipHi++ )
		{
			for( ipLo=0; ipLo < NHE1LVL; ipLo++ )
			{
				he1str.peshe1[ipHi][ipLo] = 0.;
			}
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->strk()\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/(TePowers.sqrte*TePowers.te10*TePowers.te10*
	  TePowers.te03*TePowers.te01*TePowers.te01)*phycon.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*(TePowers.sqrte/TePowers.te10/TePowers.te03);

	/* Lyman lines always have outer optical depths */
	hydro.pestrk[IP1S][IP2P] = 
		(float)(strkla/2.*fmax(1.,
		pow(HydroLines[0][IP2P][IP1S].TauIn,-0.75)));

	hydro.pestrk[IP1S][IP2P] = (float)MIN2(.01,hydro.pestrk[IP1S][IP2P]);

	hydro.pestrk[IP2P][IP1S] = 
		hydro.pestrk[1][2]*HydroLines[0][IP2P][IP1S].Aul;

	he1str.peshe1[0][1] = (float)(strkla/2.*(pow(MAX2(1.,he1nxtCOM.he1nxt[0][1]),-0.75) + 
	  pow(vfmax(1.,he1tauCOM.he1lim[0][1]/10.,he1tauCOM.he1lim[0][1]-
	  he1nxtCOM.he1nxt[0][1],FEND),-0.75)));

	he1str.peshe1[0][1] = (float)MIN2(.01,he1str.peshe1[0][1]);

	he1str.peshe1[1][0] = he1str.peshe1[0][1]*he1as.He1EinA[0][1];

	for( ipHi=3; ipHi <= nhlevel; ipHi++ )
	{
		hydro.pestrk[1][ipHi] = (float)(stark*hydro.strkar[1][ipHi]/2.*pow(MAX2(1.,
		  HydroLines[0][ipHi][IP1S].TauIn),-0.75));

		hydro.pestrk[1][ipHi] = (float)MIN2(.01,hydro.pestrk[1][ipHi]);
		hydro.pestrk[ipHi][1] = HydroLines[0][ipHi][IP1S].Aul*
		  hydro.pestrk[1][ipHi];
	}

	for( ipHi=2; ipHi < NHE1LVL; ipHi++ )
	{
		he1str.peshe1[0][ipHi] = (float)(strkla/2.*(pow(MAX2(1.,he1nxtCOM.he1nxt[0][ipHi]),-0.75) + 
		  pow(vfmax(1.,he1tauCOM.he1lim[0][ipHi]/10.,he1tauCOM.he1lim[0][ipHi]-
		  he1nxtCOM.he1nxt[0][ipHi],FEND),-0.75)));

		he1str.peshe1[0][ipHi] = (float)MIN2(.01,he1str.peshe1[0][ipHi]);
		he1str.peshe1[ipHi][0] = he1str.peshe1[0][ipHi]*he1as.He1EinA[0][ipHi];
	}

	/* all other lines */
	for( ipLo=IP2S; ipLo <= (nhlevel - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
		{
			hydro.pestrk[ipLo][ipHi] = 
				(float)fmin(.01,stark*hydro.strkar[ipLo][ipHi]*
			  pow(MAX2(1.,HydroLines[0][ipHi][ipLo].TauIn),-0.75));

			hydro.pestrk[ipHi][ipLo] = HydroLines[0][ipHi][ipLo].Aul*
			  hydro.pestrk[ipLo][ipHi];
		}
	}

	for( ipLo=1; ipLo < (NHE1LVL - 1); ipLo++ )
	{
		for( ipHi=ipLo + 1; ipHi < NHE1LVL; ipHi++ )
		{
			he1str.peshe1[ipLo][ipHi] = (float)fmin(.01,stark*hydro.strkar[ipLo+1][ipHi+1]*
			  pow(MAX2(1.,he1nxtCOM.he1nxt[ipLo][ipHi]),-0.6));
			he1str.peshe1[ipHi][ipLo] = he1as.He1EinA[ipLo][ipHi]*he1str.peshe1[ipLo][ipHi];
		}
	}


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