/*HydroOTS evaluate model hydrogen atom ots rates */
#include "cddefines.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "ionrange.h"
#include "opac.h"
#include "pbowen.h"
#include "elmton.h"
#include "phycon.h"
#include "ionfracs.h"
#include "trace.h"
#include "nh.h"
#include "rfield.h"
#include "heots.h"
/* following three needed for some debug statements */
#include "radius.h"
#include "zonecnt.h"
#include "e2tau.h"
#include "addotscon.h"
#include "addotslin.h"
/* the Bowen HeII yield */
#define BOWEN 0.3f

void HydroOTS(void)
{
	long int ipHi, 
	  ipLo, 
	  ipZ, 
	  n;
	float bwnfac ,
	  ots660 ,
	  otsnew;
	

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

	/**************************************************************************
	 *
	 * the bowen HeII - OIII fluorescense problem
	 *
	 **************************************************************************/

	ipZ = 1;
	if( elmton.lgElmtOn[ipZ] )
	{
		/* conversion per unit atom to OIII, at end of sub we divide by it,
		 * to fix lines back to proper esc/dest probs */
		bwnfac = BOWEN * (1.f - HydroLines[ipZ][IP2P][IP1S].Pesc );

		/* escaping part that will be added to the outward beam in metdif,
		 * this uses ipLnEsc instead of ipLnDesp */
		heots.esc660 = HydroLines[ipZ][IP2P][IP1S].Aul*
		  HydroLines[ipZ][IP2P][IP1S].PopHi*xIonFracs[ipZ+2][ipZ]*
		  HydroLines[ipZ][IP2P][IP1S].Pesc *BOWEN*2.0f * 0.8f ;

		/* the last factor accounts for fact that two photons are produced,
		 * and the branching ratio */
		ots660 = HydroLines[ipZ][IP2P][IP1S].Aul*
		  HydroLines[ipZ][IP2P][IP1S].PopHi*xIonFracs[ipZ+2][ipZ]*
		  HydroLines[ipZ][IP2P][IP1S].Pdest *BOWEN*2.0f * 0.8f ;

		/* now add this to the ots field */
		AddOTSLin(ots660 , pbowen.ip660 );

		/* decrease the destruction prob by the amount we will add elsewhere,
		 * ok since dest probs updated on every iteration*/
		HydroLines[ipZ][IP2P][IP1S].Pdest *= bwnfac;
	}

	else
	{
		bwnfac = 1.;
	}

	/* make ots fields due to hydrogen lines and continua */

	/* loop over all elements */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		/* do if this stage exists */
		if( (IonRange.IonHigh[ipZ] == ipZ + 2)  )
		{

			/**************************************************************************
			 *
			 * generate hydrogenic line ots rates 
			 *
			 **************************************************************************/

			/* now loop over all possible levels, but cannot include two photon
			 * since there is no pointer to this continuum */
			for( ipHi=IP2P; ipHi <= nhlevel; ipHi++ )
			{
				for( ipLo=IP1S; ipLo < ipHi; ipLo++ )
				{
					/* ots rates, the destp prob was set in hydropesc */
					HydroLines[ipZ][ipHi][ipLo].ots = 
						HydroLines[ipZ][ipHi][ipLo].Aul*
						HydroLines[ipZ][ipHi][ipLo].PopHi*xIonFracs[ipZ+2][ipZ]*
						HydroLines[ipZ][ipHi][ipLo].Pdest;

					assert( HydroLines[ipZ][ipHi][ipLo].ots >= 0. );

					/* finally dump the ots rate into the stack */
					AddOTSLin(HydroLines[ipZ][ipHi][ipLo].ots,
						HydroLines[ipZ][ipHi][ipLo].ipCont );
				}
			}
			{
				/* debugging code for Lya problems */
				enum {DEBUG = FALSE};
				if( DEBUG )
				{
					if( ipZ==0 ) 
					{
						fprintf(ioQQQ,"%li hlyaots\t%g\t%g\t%g\t%g\t%g\n",
							ZoneCnt.nzone, radius.depth ,
							HydroLines[0][IP2P][IP1S].Pdest,
							HydroLines[0][IP2P][IP1S].ots , 
							opac.opac[HydroLines[0][IP2P][IP1S].ipCont-1] ,
							rfield.otslinNew[HydroLines[0][IP2P][IP1S].ipCont-1] );
					}
				}
			}

			/**************************************************************************
			 *
			 * ots recombination bound-free b-f continua continuum
			 *
			 **************************************************************************/

			/* put in OTS continuum */
			for( n=IP1S; n <= nhlevel; n++ )
			{
				otsnew = (float)(hydro.hrec[ipZ][n][ipRecRad]*
					(1. - hydro.hrec[ipZ][n][ipRecEsc])*phycon.eden*
					xIonFracs[ipZ+2][ipZ]);
				assert( otsnew >= 0. );

				/*otsnew = (rfield.otscon[nh.ipHn[ipZ][n]-1]*
					(1.f - damp) + 
					damp*hydro.hrec[ipZ][n][ipRecRad]*
					(1.f - hydro.hrec[ipZ][n][ipRecEsc])*phycon.eden*
					xIonFracs[ipZ+2][ipZ]/
					opac.opac[nh.ipHn[ipZ][n]-1])/2.f;*/

				/* pointer used in this routine is decremented by one there */
				AddOTSCon(otsnew,nh.ipHn[ipZ][n]);
				/* debugging code for rec continua */
				{
					enum {DEBUG = FALSE};
					if( DEBUG )
					{
						if( ipZ==1 && n==IP1S) 
						{
							fprintf(ioQQQ,"%li h con ots\t%g\t%g\t%g\t%g\t%g\t%g\n",
								ZoneCnt.nzone, radius.depth ,
								otsnew,
								hydro.hrec[ipZ][n][ipRecEsc] , 
								opac.opac[nh.ipHn[ipZ][n]-1] ,
								rfield.otsconNew[nh.ipHn[ipZ][n]-1] ,
								e2tau.e2TauAbs[nh.ipHn[ipZ][n]-1] );
						}
					}
				}
			}
		}
	}

	ipZ = 1;
	if( elmton.lgElmtOn[ipZ] && bwnfac > 0. )
	{
		/* increase the destruction prob by the amount we wdecreased it above */
		HydroLines[ipZ][IP2P][IP1S].Pdest /= bwnfac;
	}

	if( trace.lgTrace )
	{
		fprintf(ioQQQ,"     HydroOTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n",
			HydroLines[0][IP2P][IP1S].Pdest , HydroLines[0][IP2P][IP1S].ots , 
			rfield.otslinNew[HydroLines[0][IP2P][IP1S].ipCont-1]  ,
			opac.opac[HydroLines[0][IP2P][IP1S].ipCont-1]
			);
	}

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