/*p8446 drive the solution of OI level populations, Ly-beta pumping */
#include "cddefines.h"
#include "taulines.h"
#include "fluors.h"
#include "h.h"
#include "doppvel.h"
#include "bowoi.h"
#include "trace.h"
#include "popoi.h"
#include "touton.h"
#include "ionfracs.h"
#include "rtescprob.h"
#include "addotslin.h"
#include "hydrogenic.h"
#include "oilevl.h"
#include "p8446.h"

void p8446(double *coloi)
{
	long int i;
	double esin, 
	  eslb, 
	  esoi, 
	  flb, 
	  foi, 
	  opaclb, 
	  opacoi, 
	  tin, 
	  tout, 
	  xlb, 
	  xoi;
	static double esab;
	static double aoi = 7.24e7;
	static double alb = 5.575e7;

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

	/* A's from Pradhan; OI pump line; Ly beta, 8446 */

	/* Notation;
	 * PMPH31 net rate hydrogen level 3 depopulated to 1
	 * PMPO15 and PMPO51 are similar, but for oxygen
	 * ESCH31 is emission in 1025 transition */

	/* called by LINES before calc really start, protect here
	 * also used for cases where OI not present */
	if( xIonFracs[1][7] <= 0. )
	{
		for( i=0; i < 6; i++ )
		{
			popoiCom.popoi[i] = 0.;
		}
		/* return zero */
		*coloi = 0.;
		bowoi.pmpo15 = 0.;
		bowoi.pmpo51 = 0.;
		bowoi.pmph31 = HydroLines[0][3][IP1S].Aul*
		  HydroLines[0][3][IP1S].Pesc;
		bowoi.esch31 = HydroLines[0][3][IP1S].Aul*
		  HydroLines[0][3][IP1S].Pesc;

		/* all trace output turned on with "trace ly beta command' */
		if( trace.lgTr8446 && trace.lgTrace )
		{
			fprintf( ioQQQ, 
				"       P8446 called for first time, finds A*escape prob 3-1 =%10.2e\n", 
			  bowoi.esch31 );
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->p8446()\n", debug_fp );
#		endif
		return;
	}

	if( touton.lgTauOutOn )
	{
		/* two sided escape prob */
		tin = HydroLines[0][3][IP1S].TauIn + 
			TauLines[ipTO1025].TauIn;
		esin = esccom(tin,1e-4);
		tout = (HydroLines[0][3][IP1S].TauTot + 
		  TauLines[ipTO1025].TauTot)*0.9 - 
		  HydroLines[0][3][IP1S].TauIn - 
		  TauLines[ipTO1025].TauIn;

		if( trace.lgTr8446 && trace.lgTrace )
		{
			fprintf( ioQQQ, "       P8446 tin, tout=%10.2e%10.2e\n", 
			  tin, tout );
		}

		if( tout > 0. )
		{
			tout = HydroLines[0][3][IP1S].TauTot + 
			  TauLines[ipTO1025].TauTot - 
			  HydroLines[0][3][IP1S].TauIn - 
			  TauLines[ipTO1025].TauIn;

			/* do not update esab if we overran optical depth scale */
			esab = 0.5*(esin + esccom(tout,1e-4));
		}
	}
	else
	{
		/* one-sided escape probability */
		esab = esccom(HydroLines[0][3][IP1S].TauIn+
		  TauLines[ipTO1025].TauIn,1e-4);
	}

	esoi = TauLines[ipTO1025].Pesc;
	eslb = HydroLines[0][3][IP1S].Pesc;

	/* all trace output turned on with "trace ly beta command' */
	if( trace.lgTr8446 && trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"       P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n", 
		  DoppVel.doppler[0], DoppVel.doppler[7], eslb, esoi, esab );
	}

	/* find relative opacities for two lines */
	opacoi = 2.92e-9*xIonFracs[1][7]*0.5556/DoppVel.doppler[7];
	opaclb = 1.22e-8*h.hii*hydro.hn[0][IP1S]/DoppVel.doppler[0];

	/* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */
	xoi = opacoi/(opacoi + opaclb);
	xlb = opaclb/(opacoi + opaclb);

	/* find relative line-widths, assume same rest freq */
	foi = MIN2(DoppVel.doppler[0],DoppVel.doppler[7])/DoppVel.doppler[7];
	flb = MIN2(DoppVel.doppler[0],DoppVel.doppler[7])/DoppVel.doppler[0]*
	  (1. - TauLines[ipTO1025].Pesc);

	if( trace.lgTr8446 && trace.lgTrace )
	{
		fprintf( ioQQQ, 
			"       P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n", 
		  opaclb, opacoi, xlb, xoi, flb, foi );
	}

	/* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate
	 * lgFluor set false with no induced command, usually true */
	if( fluors.lgFluor )
	{
		bowoi.pmpo15 = (float)(flb*alb*(hydro.hn[0][3]*h.hii)*xoi*(1. - esab)/
		  xIonFracs[1][7]);
		/* net decay rate from upper level */
		bowoi.pmpo51 = (float)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - 
		  esab)*foi));
	}
	else
	{
		bowoi.pmpo15 = 0.;
		bowoi.pmpo51 = 0.;
	}

	/* find level populations for OI */
	oilevl(xIonFracs[1][7],coloi);

	/* escape term from n=3 of H; followed by pump term to n=3 */
	bowoi.pmph31 = (float)(alb*(1. - (1. - flb)*(1. - eslb) - xlb*(1. - esab)*
	  flb));

	/*pmph13 = xlb*(1. - esab)*foi*aoi*popoiCom.popoi[4];*/
	bowoi.pmph31 = HydroLines[0][3][IP1S].Aul*
		HydroLines[0][3][IP1S].Pesc;

	/* following is actual emission rate, used to predict Ly-beta in lines.for */
	bowoi.esch31 = (float)(alb*(eslb*(1. - flb) + esab*flb));

	/* all trace output turned on with "trace ly beta command' */
	if( trace.lgTr8446 && trace.lgTrace )
	{
		fprintf( ioQQQ, "       P8446 PMPH31=%10.2e EscP3-1%10.2e ESCH31=%10.2e\n", 
		  bowoi.pmph31, bowoi.pmph31/alb, bowoi.esch31/alb );
	}

	/* continuum pumping due to J=1, 0 sub states.
	 * neglect J=2 since L-Beta very optically thick */

	/* lower level populations */
	TauLines[ipT1304].PopLo = popoiCom.popoi[0];
	TauLines[ipTO1025].PopLo = popoiCom.popoi[0];
	TauLines[ipT1039].PopLo = popoiCom.popoi[0];
	TauLines[ipT8446].PopLo = popoiCom.popoi[1];
	TauLines[ipT4368].PopLo = popoiCom.popoi[1];
	TauLines[ipTOI13].PopLo = popoiCom.popoi[2];
	TauLines[ipTOI11].PopLo = popoiCom.popoi[2];
	TauLines[ipTOI29].PopLo = popoiCom.popoi[3];
	TauLines[ipTOI46].PopLo = popoiCom.popoi[4];

	/* upper level populations */
	TauLines[ipT1304].PopHi = popoiCom.popoi[1];
	TauLines[ipTO1025].PopHi = popoiCom.popoi[4];
	TauLines[ipT1039].PopHi = popoiCom.popoi[3];
	TauLines[ipT8446].PopHi = popoiCom.popoi[2];
	TauLines[ipT4368].PopHi = popoiCom.popoi[5];
	TauLines[ipTOI13].PopHi = popoiCom.popoi[3];
	TauLines[ipTOI11].PopHi = popoiCom.popoi[4];
	TauLines[ipTOI29].PopHi = popoiCom.popoi[5];
	TauLines[ipTOI46].PopHi = popoiCom.popoi[5];

	/* opacity factor
	 * t1304(ipLnPopOpc) = (popoi(1)-popoi(2)*3.0)
	 * 666 error! following needed to get badbugs/bug8.in to work */
	TauLines[ipT1304].PopOpc = popoiCom.popoi[0];
	TauLines[ipTO1025].PopOpc = (float)(popoiCom.popoi[0] - popoiCom.popoi[4]*0.6);
	TauLines[ipT1039].PopOpc = (float)(popoiCom.popoi[0] - popoiCom.popoi[3]*3.0);

	/* t8446(ipLnPopOpc) = (popoi(2)-popoi(3)*.33)
	 * 666 error! following needed to get badbugs/bug5.in to work */
	TauLines[ipT8446].PopOpc = 
		(float)(MAX2(0.,popoiCom.popoi[1]-popoiCom.popoi[2]*.33));
	TauLines[ipT4368].PopOpc = (float)(popoiCom.popoi[1] - popoiCom.popoi[5]*.33);
	TauLines[ipTOI13].PopOpc = (float)(popoiCom.popoi[2] - popoiCom.popoi[3]*3.0);
	TauLines[ipTOI11].PopOpc = (float)(popoiCom.popoi[2] - popoiCom.popoi[4]*0.6);
	TauLines[ipTOI29].PopOpc = (float)(popoiCom.popoi[3] - popoiCom.popoi[5]*.33);
	TauLines[ipTOI46].PopOpc = (float)(popoiCom.popoi[4] - popoiCom.popoi[5]*1.67);

	AddOTSLin(popoiCom.popoi[1]*TauLines[ipT1304].Pdest*TauLines[ipT1304].Aul,
	  TauLines[ipT1304].ipCont);

	/* call AddOTSLin( popoi(4)*t1039(ipLnDesP)*t1039(ipLnAul),
	 *  1 INT(t1039(ipLnIpCont)) ,0.5) */
	AddOTSLin(popoiCom.popoi[2]*TauLines[ipT8446].Pdest*TauLines[ipT8446].Aul,
	  TauLines[ipT8446].ipCont);

	AddOTSLin(popoiCom.popoi[5]*TauLines[ipT4368].Pdest*TauLines[ipT4368].Aul,
	  TauLines[ipT4368].ipCont);

	AddOTSLin(popoiCom.popoi[3]*TauLines[ipTOI13].Pdest*TauLines[ipTOI13].Aul,
	  TauLines[ipTOI13].ipCont);

	AddOTSLin(popoiCom.popoi[4]*TauLines[ipTOI11].Pdest*TauLines[ipTOI11].Aul,
	  TauLines[ipTOI11].ipCont);

	AddOTSLin(popoiCom.popoi[5]*TauLines[ipTOI29].Pdest*TauLines[ipTOI29].Aul,
	  TauLines[ipTOI29].ipCont);

	AddOTSLin(popoiCom.popoi[5]*TauLines[ipTOI46].Pdest*TauLines[ipTOI46].Aul,
	  TauLines[ipTOI46].ipCont);

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