/*HeTran derive escape and destruction probabilities for He lines */
#include "cddefines.h"
#include "taulines.h"
#include "eshe1l.h"
#include "nhe1lvl.h"
#include "zonecnt.h"
#include "doppvel.h"
#include "escdst.h"
#include "touton.h"
#include "he.h"
#include "he1as.h"
#include "he1dmp.h"
#include "he1esc.h"
#include "phe1lv.h"
#include "he1tau.h"
#include "he3lines.h"
#include "he3tau.h"
#include "he1nxt.h"
#include "iphe1l.h"
#include "twoway.h"
#include "wind.h"
#include "rtescprob.h"
#include "eovrlp.h"
#include "he1.h"
#include "hetran.h"

void HeTran(void)
{
	long int i, 
	  j;

	static double DstPrb10830;
	double damp, 
	  esin, 
	  esout, 
	  tout;

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

	/* evaluate helium destruction probabilities
	 * this code was heavily modified by Paul T. O'Brien, UCL, to
	 * both inward and outward fractional escape
	 * the original code had only one-sided escapes */

	if( he1tauCOM.he1lim[0][1]*0.90 - he1nxtCOM.he1nxt[0][1] > 0. )
	{
		/* only evaluate if we have not overrun optical depth scale */
		he1escCOM.he1esc[0][1] = (float)eshe1l(&esin);
		he1escCOM.he1esc[1][0] = he1escCOM.he1esc[0][1]*he1as.He1EinA[0][1];
		/*he1fin[0][1] = twoway.fracin;*/
		he1tauCOM.he1con[0][1] = 
			(float)(he1jbr(2,1)*he1as.He1EinA[0][1]*esca0k2(he1nxtCOM.he1nxt[0][1]));

		damp = he1dmpCOM.he1dmp[1]/DoppVel.doppler[1]/0.75;

		/* >>chng 96 sep 4, added DestWght weighting
		 * >>chng 96 sep 4, following had extra div by two???? */
		phe1lv.he1dst[0][1] = 
			(float)( he1as.He1EinA[0][1]*
		  eovrlp(he.heii*phe1lv.he1n[0], he1tauCOM.he1opc[0][1],
		  iphe1lCom.iphe1l[0][1],DoppVel.doppler[1],he1escCOM.he1esc[0][1],
		  "  K2"));
		/* >>chng 98 July 24, from "INCO" to "  K2", 
		 * strongly affects 2.06 micron line*/
	}
	else
	{
		/* >>chng 97 sep 13, did not conserve energy, added following */
		he1tauCOM.he1con[0][1] = 
			(float)(he1jbr(2,1)*he1as.He1EinA[0][1]*esca0k2(he1nxtCOM.he1nxt[0][1]));
	}

	for( i=2; i < NHE1LVL; i++ )
	{
		damp = he1dmpCOM.he1dmp[i]/DoppVel.doppler[1]/1.808/(1. - 
		  1./POW2(i+1.));

		esin = esccom(he1nxtCOM.he1nxt[0][i],damp);
		tout = he1tauCOM.he1lim[0][i]*0.90 - he1nxtCOM.he1nxt[0][i];

		if( tout > 0. )
		{
			tout = he1tauCOM.he1lim[0][i] - he1nxtCOM.he1nxt[0][i];
			esout = esccom(tout,damp);

			he1escCOM.he1esc[0][i] = (float)(0.5*(esin + esout));
			/*he1fin[0][i] = (float)(esin/(esin + esout));*/

			/* other index is product of A and escape */
			he1escCOM.he1esc[i][0] = he1escCOM.he1esc[0][i]*he1as.He1EinA[0][i];

			/* >>chng 96 sep 4, added DestWght weighting */
			phe1lv.he1dst[0][i] = 
				(float)(he1as.He1EinA[0][i]*
			  eovrlp(he.heii* phe1lv.he1n[0],
			  he1tauCOM.he1opc[0][i],iphe1lCom.iphe1l[0][i],
			  DoppVel.doppler[1],he1escCOM.he1esc[0][i],escdst.chDstFun  ));
		}

		he1tauCOM.he1con[0][i] = 
			(float)(he1jbr(i+1,1)*esca0k2(he1nxtCOM.he1nxt[0][i])*
		  he1as.He1EinA[0][i]);
	}

	if( touton.lgTauOutOn )
	{
		for( j=1; j < (NHE1LVL - 1); j++ )
		{
			for( i=j + 1; i < NHE1LVL; i++ )
			{
				damp = he1dmpCOM.he1dmp[i]/DoppVel.doppler[1]/1.808/
				  (1./POW2(j+1.) - 1./POW2(i+1.));

				esin = esccom(he1nxtCOM.he1nxt[j][i],damp);

				tout = he1tauCOM.he1lim[j][i]*0.90 - he1nxtCOM.he1nxt[j][i];

				if( tout > 0. )
				{
					tout = he1tauCOM.he1lim[j][i] - he1nxtCOM.he1nxt[j][i];
					esout = esccom(tout,damp);

					/*he1fin[j][i] = (float)(esin/(esin + esout));*/
					he1tauCOM.he1con[j][i] = (float)(he1jbr(i+1,j+1)*
						esin*he1as.He1EinA[j][i]);

					he1escCOM.he1esc[j][i] = (float)(0.5*(esin + esout));
					he1escCOM.he1esc[i][j] = he1escCOM.he1esc[j][i]*
					  he1as.He1EinA[j][i];

					/* destruction of line by background opacity
					 * >>chng 96 sep 4, added DestWght weighting */
					phe1lv.he1dst[j][i] = 
						(float)(he1as.He1EinA[j][i]*
					  eovrlp(he.heii*phe1lv.he1n[j],he1tauCOM.he1opc[j][i],
					  iphe1lCom.iphe1l[j][i],DoppVel.doppler[1],
					  he1escCOM.he1esc[j][i],escdst.chDstFun));
				}
			}
		}

		esin = esccom(t206.TauIn,1e-3);
		esout = esccom(t206.TauTot,1e-3);
		t206.Pesc = (float)(0.5*(esin + esout));
		t206.FracInwd = (float)(esin/(esin + esout));
	}

	else
	{

		/* outward optical depths not defined */
		for( j=1; j < (NHE1LVL - 1); j++ )
		{
			for( i=j + 1; i < NHE1LVL; i++ )
			{
				damp = he1dmpCOM.he1dmp[i]/DoppVel.doppler[1]/1.808/
				  (1./POW2(j+1.) - 1./POW2(i+1.));

				esin = esccom(he1nxtCOM.he1nxt[j][i],damp);
				/*he1fin[j][i] = 1.;*/
				he1tauCOM.he1con[j][i] = (float)(he1jbr(i+1,j+1)*
					esin*he1as.He1EinA[j][i]);

				he1escCOM.he1esc[j][i] = (float)esin;

				he1escCOM.he1esc[i][j] = he1escCOM.he1esc[j][i]*
				  he1as.He1EinA[j][i];

				/* destruction of line by background opacity
				 * >>chng 96 sep 4, added DestWght weighting */
				phe1lv.he1dst[j][i] = 
					(float)(he1as.He1EinA[j][i]*
				  eovrlp(he.heii*phe1lv.he1n[j],he1tauCOM.he1opc[j][i],iphe1lCom.iphe1l[j][i],
				  DoppVel.doppler[1],he1escCOM.he1esc[j][i],escdst.chDstFun
				  ));
			}
		}

		t206.Pesc = (float)esccom(t206.TauIn,1e-3);
		t206.FracInwd = 1.;
	}

	if( ZoneCnt.nzone <= 1 )
	{
		DstPrb10830 = eovrlp(he3lines.p2s,8.80e-7,he3lines.ipHe3l[0],
		  DoppVel.doppler[1],he3tau[0].Pesc,escdst.chDstFun );
	}
	else
	{
		DstPrb10830 = eovrlp(he3lines.p2s,8.80e-7,
		  he3lines.ipHe3l[0],DoppVel.doppler[1],he3tau[0].Pesc,
		  escdst.chDstFun);
	}

	he3tau[IPT10830-1].Pdest = (float)DstPrb10830;

	if( wind.windv == 0. )
	{
		for( i=0; i < NHE3TAU; i++ )
		{
			/* always use complete redis for helium lines */
			he3tau[i].Pesc = 
				(float)escsub(he3tau[i].TauIn,
			  he3tau[i].TauTot,1e-4);

			he3tau[i].FracInwd = twoway.fracin;
		}
	}

	else
	{
		/* this is wind solution
		 * ESCCOM is one-sided complete distribution (OK for wind) */
		for( i=0; i < NHE3TAU; i++ )
		{
			he3tau[i].Pesc = 
				(float)esccom(he3tau[i].TauIn,  1e-4);
			he3tau[i].FracInwd = 0.5;
		}
	}

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