/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*highen do high energy radiation field - gas interaction, Compton scattering, etc */
#include "cddefines.h"
#include "physconst.h"
#include "trace.h"
#include "heavy.h"
#include "radius.h"
#include "hextra.h"
#include "hmi.h"
#include "thermal.h"
#include "dense.h"
#include "neutrn.h"
#include "doppvel.h"
#include "ionbal.h"
#include "rfield.h"
#include "opacity.h"
#include "gammas.h"
#include "highen.h"

void highen( void )
{
	long int i,
		ion,
		nelem;

	double CosRayDen, 
	  crsphi, 
	  heat_cold_electrons, 
	  heatin, 
	  sqthot;

	double hin;

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


	/**********************************************************************
	 *
	 *						COMPTON RECOIL IONIZATION
	 *
	 * 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 
	 *
	 **********************************************************************/

	/* comoff turned off with no compton command,
	 * lgCompRecoil - no recoil ionization */
	if( rfield.comoff > .0 && ionbal.lgCompRecoil )
	{
		for( nelem=0; nelem<LIMELM; ++nelem )
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				ionbal.CompRecoilIonRate[nelem][ion] = 0.;
				ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
				if( dense.xIonDense[nelem][ion] > 0. )
				{
					/* recoil ionization starts at 194 Ryd = 2.6 keV */
					/* this is the ionization potential of the valence shell */
					/* >>chng 02 may 27, lower limit is now 1 beyond actual threshold
					 * since recoil energy at threshold was very small, sometimes negative */
					for( i=ionbal.ipCompRecoil[nelem][ion]; i < rfield.nflux; ++i)
					{
						double recoil_energy;
						crsphi = opac.OpacStack[i-1+opac.iopcom] * rfield.SummedCon[i] *
							/* this is number of electrons in valence shell of this ion */
							ionbal.nCompRecoilElec[nelem-ion];

						/* direct hydrogen ionization due to compton scattering, 
						 * does not yet include secondaries,
						 * last term accounts for number of valence electrons that contribute */
						ionbal.CompRecoilIonRate[nelem][ion] += crsphi ;

						/* recoil energy in Rydbergs
						 * heating modified for suprathermal secondaries below; ANU2=ANU**2 */
						/* >>chng 02 mar 27, use general formula for recoil energy */
						/*energy = 2.66e-5*rfield.anu2[i] - 1.;*/
						recoil_energy  = rfield.anu2[i] / 
							( EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT) * EN1RYD* EN1RYD - 
							Heavy.Valence_IP_Ryd[nelem][ion];

						/* heating is in rydbergs because efionz, exctef, heatef in ryd */
						ionbal.CompRecoilHeatRate[nelem][ion] += crsphi*recoil_energy;
					}
					/* net heating rate, per atom, convert ryd/sec/cm3 to ergs/sec/atom */
					ionbal.CompRecoilHeatRate[nelem][ion] *= EN1RYD;

					ASSERT( ionbal.CompRecoilHeatRate[nelem][ion] >= 0.);
					ASSERT( ionbal.CompRecoilIonRate[nelem][ion] >= 0.);
				}
			}
		}
	}
	else
	{
		for( nelem=0; nelem<LIMELM; ++nelem )
		{
			for( ion=0; ion<nelem+1; ++ion )
			{
				ionbal.CompRecoilIonRate[nelem][ion] = 0.;
				ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
			}
		}
	}

	/**********************************************************************
	 *
	 *                          COSMIC RAYS
	 *
	 * 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, as determined from original density times
		 * possible dependence on radius */
		CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,(double)hextra.crpowr);

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

		/* rate hot electrons heat cold ones, Balbus and McKee 1982
		 * units 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 */
		heat_cold_electrons = 5.5e-14/sqthot*CosRayDen*
			/*dense.eden/(dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHELIUM][0] + hmi.Hmolec[ipMH2g]);*/
			dense.eden/(dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHELIUM][0] + hmi.H2_total);

		/* ionization rate; Balbus and McKee */
		ionbal.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 )
		{
			/* usual circumstance; relativistic cosmic rays, 
			 * cosmic ray ionization rate s-1 cm-3; ext rel limit */
			ionbal.CosRayIonRate *= 3.;

		}
		else
		{
			/*  option for thermal cosmic rays */
			ionbal.CosRayIonRate *= 10.;
		}
		/* u factor is 15 eV per ionization,
		 * 13.6 into ionization itself, rest into non-ionizing heat, 
		 * (15-13.6)/13.6=0.103, heating in Ryd, cnvrt to ergs eith EN1RYD
		 * will multiply by number of secondaries that occur when used in heatsum*/
		/* >>chng 04 jan 27, from 0.093 to 2.574 as per following */
		/* cr heating from Table 1 of
		 *>>refer	cr	heating	Wolfire et al.1995, ApJ, 443, 152
		 * For every ionization due to cosmic rays, ~35eV of heat is added 
		 * to the system.  This manifests itself in the ionbal.CosRayHeatRate term
		 * by the 2.574*EN1RYD term, which is just the energy in ergs in 35 eV.
		 * Change made by Nick Abel and Gargi Shaw, 04 Jan 27.   In heatsum 
		 * we  multiply by the number of secondaries that occur */
		/* >>chng 04 mar 24, div last term by ratio of sec to prim since do not want
		 * to multiply by sec to prim in final result, in heatsum */
		ionbal.CosRayHeatRate = ionbal.CosRayIonRate*2.574*EN1RYD + heat_cold_electrons;

		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, ionbal.CosRayIonRate, ionbal.CosRayHeatRate, hextra.crtemp );
		}
	}
	else
	{
		ionbal.CosRayIonRate = 0.;
		ionbal.CosRayHeatRate = 0.;
	}
	/*fprintf(ioQQQ,"DEBUG cr %.2f %.3e %.3e %.3e\n",
		fnzone,
		hextra.cryden ,
		ionbal.CosRayIonRate ,
		ionbal.CosRayHeatRate );*/

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

	/* TurbHeat added with hextra command, DispScale added with turbulence dissipative */
	if( (hextra.TurbHeat+DoppVel.DispScale) != 0. )
	{
		/* turbulent heating only goes into the low-energy heat of this element */
		/* >>>>chng 00 apr 28, functional form of radius dependence had bee turrad/depth
		 * and so went to infinity at the illuminated face.  Changed to current form as
		 * per Ivan Hubeny comment */
		if( hextra.turrad > 0. )
		{
			/* if turrad is >0 then vary heating with depth */
			ionbal.ExtraHeatRate = 
				hextra.TurbHeat*sexp(radius.depth /hextra.turrad);

			/* >>chng 00 aug 16, added option for heating from back side */
			if( hextra.turback != 0. )
			{
				ionbal.ExtraHeatRate += 
					hextra.TurbHeat*sexp((hextra.turback - radius.depth) /hextra.turrad);
			}
		}
		/* this is turbulence dissipate command */
		else if( DoppVel.DispScale > 0. )
		{
			double turb = DoppVel.TurbVel * sexp( radius.depth / DoppVel.DispScale );
			/* if turrad is >0 then vary heating with depth */
			/* >>chng 01 may 10, add extra factor of length over 1e13 cm */
			ionbal.ExtraHeatRate = 3.45e-28 / 2.82824 * turb * turb * turb
				/ (DoppVel.DispScale/1e13) ;
		}
		else
		{
			/* constant extra heating */
			ionbal.ExtraHeatRate = (float)hextra.TurbHeat;
		}

	}
	else
	{
		ionbal.ExtraHeatRate = 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 */
		/* neutrons assumed to only secondary ionize */
		ionbal.xNeutronHeatRate = neutrn.totneu*neutrn.CrsSecNeutron;
	}
	else
	{
		ionbal.xNeutronHeatRate = 0.;
	}


	/**********************************************************************
	 *
	 * pair production in elec field of nucleus 
	 *
	 **********************************************************************/
	ionbal.PairProducPhotoRate[0] = GammaK(opac.ippr,rfield.nflux,opac.ioppr,1.);
	ionbal.PairProducPhotoRate[1] = thermal.HeatLowEnr;
	ionbal.PairProducPhotoRate[2] = thermal.HeatHiEnr;

	/**********************************************************************
	 *
	 * Compton energy exchange 
	 *
	 **********************************************************************/
	rfield.cmcool = 0.;
	rfield.cmheat = 0.;
	heatin = 0.;
	/* comoff usually 1, turns off Compton */
	if( rfield.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 */
			rfield.comup[i] = (double)(rfield.flux[i]+rfield.ConInterOut[i]+
			  rfield.outlin[i]+ rfield.outlin_noplot[i])*rfield.csigc[i]*(dense.eden*4.e0*
			  6.338e-6*1e-15);

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

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

			/* induced Compton heating */
			hin = (double)(rfield.flux[i]+rfield.ConInterOut[i]+rfield.outlin[i]+rfield.outlin_noplot[i])*
			  rfield.csigh[i]*rfield.OccNumbIncidCont[i]*dense.eden*1e-15;
			rfield.comdn[i] += hin;
			heatin += hin;

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

		/* remember how important induced compton heating is */
		if( rfield.cmheat > 0. )
		{
			rfield.cinrat = MAX2(rfield.cinrat,heatin/rfield.cmheat);
		}

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

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

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

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

