/*highen do high energy radiation field - gas interaction, Compton scattering, etc */
#include "cddefines.h"
#include "physconst.h"
#include "opacpoint.h"
#include "trace.h"
#include "recoil.h"
#include "heating.h"
#include "radius.h"
#include "turhet.h"
#include "hextra.h"
#include "ionfracs.h"
#include "hmi.h"
#include "neutrn.h"
#include "ircoil.h"
#include "photrate.h"
#include "compton.h"
#include "rfield.h"
#include "heat.h"
#include "phycon.h"
#include "hgamsc.h"
#include "opac.h"
#include "crsnut.h"
#include "sexp.h"
#include "gammas.h"
#include "highen.h"

void highen( void )
{
	long int i;

	double CosRayDen, 
	  CosRayIonRate, 
	  CosRayHeatRate,
	  crsphi, 
	  energy, 
	  hcolde, 
	  heatin, 
	  sqthot;

	static double
	  HeatHRecoil, 
	  HeatHeRecoil;
	double hin;

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


	/**********************************************************************
	 *
	 * bound electron scattering of >2.3 kev photons if neutral
	 * comoff usually 1, 0 if "NO COMPTON EFFECT" command given
	 * lgCompRecoil usually t, f if "NO RECOIL IONIZATION" command given 
	 *
	 **********************************************************************/

	if( Compton.comoff > .0 && recoil.lgCompRecoil )
	{
		HeatHRecoil = 0.;
		hgamsc.BoundCompton = 0.;

		/* recoil starts at 194 Ryd = 2.6 keV */
		for( i=ircoil.ipRecoilIon-1; i < rfield.nflux; ++i)
		{
			crsphi = opac.OpacStack[i-1+OpacPoint.iopcom]*(rfield.flux[i] + 
			  rfield.outcon[i] + rfield.outlin[i]);

			/* direct hydrogen ionization due to compton scattering, does not include secondaries  */
			hgamsc.BoundCompton += (float)crsphi;

			/* recoil energy in Rydbergs
			 * heating modified for suprathermal secondaries below; ANU2=ANU**2 */
			energy = 2.66e-5*rfield.anu2[i] - 1.;

			/* heating is in rydbergs because efionz, exctef, heatef in ryd */
			HeatHRecoil += crsphi*energy;
		}
		/* net heating efficiency, now in ergs/sec/cm3 */
		HeatHRecoil *= EN1RYD;

		/* same thing for helium */
		HeatHeRecoil = 0.;
		for( i=ircoil.irche-1; i < rfield.nflux; ++i )
		{
			crsphi = opac.OpacStack[i-1+OpacPoint.iopcom]*(rfield.flux[i] + 
			  rfield.outcon[i] + rfield.outlin[i]);

			/* recoil energy in Rydbergs
			 * heating modified for suprathermal secondaries below; ANU2=ANU**2 */
			energy = 2.66e-5*rfield.anu2[i] - 1.8;
			HeatHeRecoil += crsphi*energy;
		}
		/* convert to ergs */
		HeatHeRecoil *= EN1RYD;
	}
	else
	{
		HeatHRecoil = 0.;
		HeatHeRecoil = 0.;
		hgamsc.BoundCompton = 0.;
	}

	/* bound compton off hydrogen */
	PhotRate.PhotoRate[0][0][13][0] = 0.;
	PhotRate.PhotoRate[1][0][13][0] = 0.;
	PhotRate.PhotoRate[2][0][13][0] = HeatHRecoil;
	/* bound compton off helium */
	PhotRate.PhotoRate[0][0][14][0] = 0.;
	PhotRate.PhotoRate[1][0][14][0] = 0.;
	/* second term is helium, two because two electrons */
	PhotRate.PhotoRate[2][0][14][0] = HeatHeRecoil*2.;

	/**********************************************************************
	 *
	 * heating and ionization by cosmic rays, other relativistic particles
	 * CRYDEN=density (1/CM**3), neutral rate assumes 15ev total
	 * energy loss, 13.6 into ionization, 1.4 into heating
	 * units erg/sec/cm**3
	 * iff not specified, CRTEMP is 2.6E9K
	 *
	 **********************************************************************/

	if( hextra.cryden > 0. )
	{
		/* this is current cosmic ray density */
		CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,hextra.crpowr);

		/* related to current temperature, when thermal */
		sqthot = sqrt(hextra.crtemp);

		/* rate hot electrons heat cold ones, Balbus and McKee 1982
		 * erg sec^-1 cm^-3,
		 * in sumheat we will multipy this rate by sum of neturals, but for this
		 * term we actually want eden, so mult by eden over sum of neut */
		hcolde = 5.5e-14/sqthot*CosRayDen*
			phycon.eden/(xIonFracs[1][0] + xIonFracs[1][1] + hmi.htwo);

		/* ionization rate; Balbus and McKee */
		CosRayIonRate = 1.22e-4/sqthot*
			log(2.72*pow(hextra.crtemp/1e8,0.097))*CosRayDen;

		/* option for thermal CRs, first is the usual (and default) relativistic case */
		if( hextra.crtemp > 2e9 )
		{
			/*  cosmic ray ionization rate s-1 cm-3; ext rel limit */
			CosRayIonRate *= 3.;

			/*  usual circumstance; relativistic cosmic rays, factor is 15 eV per ionization,
			 * 13.6 into ionization itself, rest into non-ionizing heat, but must later
			 * multiply by number of secondaries that occur */
			CosRayHeatRate = (CosRayIonRate*0.093*EN1RYD + hcolde);
		}
		else
		{
			/*  option for thermal cosmic rays */
			CosRayIonRate *= 10.;
			CosRayHeatRate = (hcolde + CosRayIonRate*0.093*EN1RYD);
		}
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, "     HIGHEN: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n", 
			  CosRayDen, CosRayIonRate, CosRayHeatRate, hextra.crtemp );
		}
	}
	else
	{
		CosRayIonRate = 0.;
		CosRayHeatRate = 0.;
	}
	/* will save ionization and heating into this cell,
	 * first will be directly added to secondary rate */
	PhotRate.PhotoRate[0][0][12][0] = CosRayIonRate;

	/* will assume all energy is low energy since already result of secondaries */
	PhotRate.PhotoRate[1][0][12][0] = CosRayHeatRate;
	PhotRate.PhotoRate[2][0][12][0] = 0.;

	/**********************************************************************
	 *
	 * add on extra heating due to turbulence, goes into [1] of [x][0][11][0]
	 *
	 **********************************************************************/

	if( turhet.TurbHeat != 0. )
	{
		/* turbulent heating only goes into the low-energy heat of this element */
		PhotRate.PhotoRate[1][0][11][0] = 
			(float)(turhet.TurbHeat*sexp(turhet.turrad/(float)(radius.depth)));
	}
	else
	{
		PhotRate.PhotoRate[1][0][11][0] = 0.;
	}


	/**********************************************************************
	 *
	 * option to add on fast neutron heating, goes into [0] & [2] of [x][0][11][0]
	 *
	 **********************************************************************/
	if( neutrn.lgNeutrnHeatOn )
	{
		/* CrsSecNeutron is 4E-26 cm^-2, cross sec for stopping rel neutrons
		 * this is heating due to fast neutrons, assumed to secondary ionize */
		PhotRate.PhotoRate[2][0][11][0] = neutrn.totneu*crsnut.CrsSecNeutron;
	}
	else
	{
		PhotRate.PhotoRate[2][0][11][0] = 0.;
	}
	/* neutrons assumed to only secondary ionize */
	PhotRate.PhotoRate[0][0][11][0] = 0.;


	/**********************************************************************
	 *
	 * pair production in elec field of nucleus 
	 *
	 **********************************************************************/
	PhotRate.PhotoRate[0][0][10][0] = GammaK(OpacPoint.ippr,rfield.nflux,OpacPoint.ioppr,1.);
	PhotRate.PhotoRate[1][0][10][0] = heat.HeatLowEnr;
	PhotRate.PhotoRate[2][0][10][0] = heat.HeatHiEnr;


	/**********************************************************************
	 *
	 * Compton energy exchange 
	 *
	 **********************************************************************/
	Compton.cmcool = 0.;
	Compton.cmheat = 0.;
	heatin = 0.;
	/* comoff usually 1, turns off Compton */
	if( Compton.comoff >0.01 )
	{
		for( i=0; i < rfield.nflux; i++ )
		{

			/* Compton cooling
			 * CSIGC is Tarter expression times ANU(I)*3.858E-25
			 * 6.338E-6 is k in inf mass Rydbergs, still needs factor of TE */
			Compton.comup[i] = (double)(rfield.flux[i]+rfield.outcon[i]+
			  rfield.outlin[i])*Compton.csigc[i]*(phycon.eden*4.e0*
			  6.338e-6*1e-15);

			Compton.cmcool += Compton.comup[i];

			/* Compton heating 
			 * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25
			 * CMHEAT is just spontaneous, HEATIN is just induced */
			Compton.comdn[i] = (double)(rfield.flux[i]+rfield.outcon[i]+
			  rfield.outlin[i])*Compton.csigh[i]*phycon.eden*1e-15;

			/* induced Compton heating */
			hin = (double)(rfield.flux[i]+rfield.outcon[i]+rfield.outlin[i])*
			  Compton.csigh[i]*rfield.occnum[i]*phycon.eden*1e-15;
			Compton.comdn[i] += hin;
			heatin += hin;

			/* following is total compton heating */
			Compton.cmheat += Compton.comdn[i];
		}

		/* remember how important induced compton heating is */
		if( Compton.cmheat > 0. )
		{
			Compton.cinrat = MAX2((float)(Compton.cinrat),
				heatin/MAX2(1e-37, (float)(Compton.cmheat)));
		}

		if( trace.lgTrace && trace.lgComBug )
		{
			fprintf( ioQQQ, "     HIGHEN: COOL num=%8.2e HEAT num=%8.2e HCOIL=%8.2e\n", 
			  Compton.cmcool, Compton.cmheat, HeatHRecoil );
		}
	}

	/* save compton heating rate into main heating save array, 
	 * no secondary ionizations from compton heating cooling */
	HeatingCom.heating[19][0] = Compton.cmheat;

	if( trace.lgTrace && trace.lgComBug )
	{
		fprintf( ioQQQ, 
			"     HIGHEN finds heating fracs= frac(compt)=%10.2e f(bcompt)="
			"%10.2e f(pair)%10.2e totHeat=%10.2e\n", 
		  HeatingCom.heating[19][0]/heat.htot, 
		  HeatHRecoil/heat.htot, 
		  HeatingCom.heating[21][0]/heat.htot,
		  heat.htot	);
	}

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