/*HeMetalOTS compute diffuse fields due to helium atom, ion, triplets, metal recombination,
 * called by ioize */
#include "cddefines.h"
#include "showme.h"
#include "nhe3lvl.h"
#include "nhe1lvl.h"
#include "exptau.h"
#include "opac.h"
#include "he3n.h"
#include "phe1lv.h"
#include "he.h"
#include "phycon.h"
#include "tepowers.h"
#include "doppvel.h"
#include "pressure.h"
#include "heots.h"
#include "hhe2phtots.h"
#include "he1tau.h"
#include "nhe1.h"
#include "he3lines.h"
#include "d10830.h"
#include "iphe1l.h"
#include "zonecnt.h"
#include "ip626.h"
#include "he3as.h"
#include "sphere.h"
#include "addotscon.h"
#include "addotslin.h"
#include "ionfracs.h"
#include "diffheav.h"
#include "heavy.h"
#include "widlin.h"

void HeMetalOTS(void)
{
	long int i, 
		ip,
		ipla,
	  j, 
	  nelem,
	  n;
	float dam584, 
		diffuse ,
		difflya,
		esc,
	  escfac, 
	  ots,
	  outward;
	float totline,
	  otsfac;
	float otshere; 

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

	/* He I recombination continuum */
	for( n=0; n < 6; n++ )
	{
		otshere = (float)(phe1lv.he1rec[0][n]*(1. - phe1lv.he1rec[2][n])*
		  phycon.eden*he.heii) ;
		if( otshere < 0. )
		{
			fprintf( ioQQQ, " HeMetalOTS found negative OTS He1 CON, n=%4ld val=%10.1e\n", 
			  n, otshere );
			ShowMe();
			puts( "[Stop in hediff]" );
			exit(1);
		}
		AddOTSCon(otshere,nhe1Com.nhe1[n]);
	}

	/* HeI recombination continua at 1.81 Ryd */
	heots.esc504 = (float)(phe1lv.he1rec[0][0]*phe1lv.he1rec[2][0]*phycon.eden*
	  he.heii);

	/* HeI Ly-alpha He I ly a
	 * escaping part of HeI Ly-alpha lya ly a */
	heots.esc584 = 
		he.heii*phe1lv.he12p*phe1lv.he1mis[0][1]*
		exptau.ExpmTau[iphe1lCom.iphe1l[0][1]-1];

	/* do not want to add much outward HeI Lya since it is so important
	 * use inward looking optical depth to decide outward flux */
	otshere = (float)(phe1lv.he12p*he.heii*phe1lv.he1dst[0][1] + he.heii*phe1lv.he12p*
	  phe1lv.he1mis[0][1]*(1.f - exptau.ExpmTau[iphe1lCom.iphe1l[0][1]-1]));

	if( otshere < 0. )
	{
		fprintf( ioQQQ, " HeMetalOTS found negative OTS He1 alhpa,%4i val=%10.1e\n", 
		  1, otshere );
		ShowMe();
		puts( "[Stop in hediff]" );
		exit(1);
	}
	AddOTSLin(otshere,iphe1lCom.iphe1l[0][1] );

	/* all other lines of HeI - tests show this is very small effect */
	for( j=2; j < 7; j++ )
	{
		for( i=0; i <= (j - 1); i++ )
		{
			otshere = phe1lv.he1n[j]*he.heii*phe1lv.he1dst[i][j];
			if( otshere < 0. )
			{
				fprintf( ioQQQ, " HeMetalOTS found negative OTS He1l%4ld%4ld val=%10.1e\n", 
				  j, i, otshere );
				ShowMe();
				puts( "[Stop in hediff]" );
				exit(1);
			}
			AddOTSLin(otshere,iphe1lCom.iphe1l[i][j]);
		}
	}

	/* HeI 2tripS to ground - one photon decay
	 * this is a major source of ionization for parishii */
	totline = he.heii*he3nCom.he3n[0]*he3as.a2sg;
	if( sphere.lgSphere )
	{
		/* all goes outward or ots */
		escfac = 0.;
	}
	else
	{
		escfac = 0.5f*exptau.ExpmTau[ip626.ipHe23sGrnd-1];
	}

	/* number of phots escaping from illuminated face of cloud */
	heots.esc626 = (float)(totline*escfac);

	/* otsfac = 1. - ExpZone(ipHe23sGrnd)
	 * when optically thin, this goes to zero, when thick, to 100% */
	otsfac = 1.f - 2.f*escfac;

	/* >>chng 96 mar 7, above had been 1-escfac, so 50% ots when thin to cont */
	AddOTSLin(totline*otsfac,ip626.ipHe23sGrnd);

	/* following is zero when thick to cont, 0.5 when thin */
	outward = MAX2(0.f,1.f-otsfac-escfac);

	/* outward escaping line */
	heots.out626 = (float)(totline*outward);

	/* HeI 10830 OTS */
	AddOTSLin(he3lines.p2p*d10830.DR10830,he3lines.ipHe3l[0]);

	/* HeI Ly-a 584 radiation pressure */
	if( (phe1lv.he1n[0] > 1e-20 && he1tauCOM.he1tau[0][1] > 1.1) && 
	  opac.opac[iphe1lCom.iphe1l[0][1]-1] > 0. )
	{
		dam584 = 0.0921f/TePowers.sqrte;
		pressure.plines[5] = (float)(9.54e-5*phe1lv.he12p/6./(phe1lv.he1n[0]/2.)*
		  widlin(he1tauCOM.he1tau[0][1],he1tauCOM.he1lim[0][1],dam584,
		  DoppVel.doppler[1]));
	}
	else
	{
		pressure.plines[5] = 0.;
	}

	/*>>chng 99 feb 17 this code had been in metdif, moved here in cleanup */
	/* all locally destroyed two photon continua of hydrogenic H, HeI, HeII */
	if( hhe2phtots.H12phtOTS < 0. )
	{
		fprintf( ioQQQ, " metdif found negative OTS CON, H12phtOTS val=%10.1e\n", 
		  hhe2phtots.H12phtOTS );
		ShowMe();
		puts( "[Stop in metdif]" );
		exit(1);
	}
	AddOTSCon(hhe2phtots.H12phtOTS,hhe2phtots.ipH12pht);
	{
		enum {DEBUG=FALSE};
		if( DEBUG )
		{
#			include "taulines.h"
			fprintf(ioQQQ," hemetalots expmtau %e %e\n" , 
				exptau.ExpmTau[HydroLines[0][IP2P][IP1S].ipCont - 30 ],
				hhe2phtots.H12phtOTS);
		}
	}

	if( hhe2phtots.He12phtOTS < 0. )
	{
		fprintf( ioQQQ, " metdif found negative OTS CON, He12phtOTS val=%10.1e\n", 
		  hhe2phtots.He12phtOTS );
		ShowMe();
		puts( "[Stop in metdif]" );
		exit(1);
	}
	AddOTSCon(hhe2phtots.He12phtOTS,hhe2phtots.ipHe12pht);

	if( hhe2phtots.He22phtOTS < 0. )
	{
		fprintf( ioQQQ, " metdif found negative OTS CON, He22phtOTS val=%10.1e\n", 
		  hhe2phtots.He22phtOTS );
		ShowMe();
		puts( "[Stop in metdif]" );
		exit(1);
	}
	AddOTSCon(hhe2phtots.He22phtOTS,hhe2phtots.ipHe22pht);


	/* add recombination continua for elements heavier than He */
	/* nelem = 2 is Li */
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		/* do not include hydrogenic species in following */
		for( i=0; i < nelem; i++ )
		{
			if( xIonFracs[i+2][nelem] > 0. )
			{
				ip = Heavy.ipHeavy[i][nelem];

				/* this is extinction to face */
				esc = exptau.ExpmTau[ip-1];

				/* DiffusHeavy is rad rec to ground, with electron density,
				 * it is evaluated in routine MakeRecomb */
				diffuse = DiffHeav.DiffusHeavy[i][nelem]*xIonFracs[i+2][nelem];

				if( opac.opac[ip-1] > 0. )
				{
					ots = diffuse*MIN2(0.f,1.f-esc);
				}
				else
				{
					ots = 0.;
				}

				if( ots < 0. )
				{
					fprintf( ioQQQ, " metdif found negative OTS heav, n=%4ld%4ld val=%10.1e nzone=%4ld\n", 
					  nelem, i, ots, ZoneCnt.nzone );
					ShowMe();
					puts( "[Stop in metdif]" );
					exit(1);
				}

				AddOTSCon(ots,ip);

				/* now do the recombination Lya */
				ipla = Heavy.ipLyHeavy[i][nelem];
				esc = exptau.ExpmTau[ipla-1];

				/* xLyaHeavy is set to a fraction of rad rec in MakeRecomb */
				difflya = DiffHeav.xLyaHeavy[i][nelem]*xIonFracs[i+2][nelem];

				/*begin sanity check */
				if( difflya < 0. )
				{
					fprintf( ioQQQ, " METDIF Negative Lya flux for element, ion=%4ld%4ld\n", 
					  nelem, i );
					fprintf( ioQQQ, " METDIF xLyaHeavy, xIonFracs=%10.2e%10.2e\n", 
					  DiffHeav.xLyaHeavy[i][nelem], xIonFracs[i+2][nelem] );

					ShowMe();
					puts( "[Stop in metdif]" );
					exit(1);
				}
				/*end sanity check */

				ots = difflya*MIN2(0.f,1.f-esc);

				if( ots < 0. )
				{
					fprintf( ioQQQ, " metdif found negative OTS line heav%4ld%4ld val=%10.1e\n", 
					  nelem, i, ots );
					fprintf( ioQQQ, " 1-esc =%10.2e\n", 1. - 
					  esc );
					ShowMe();
					puts( "[Stop in metdif]" );
					exit(1);
				}
				AddOTSLin(ots,ipla);
			}
		}
	}


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