/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*ion_photo fill array PhotoRate with photoionization rates for heavy elements */
#include "cddefines.h"
#include "yield.h"
#include "heavy.h"
#include "opacity.h"
#include "dense.h"
#include "thermal.h"
#include "conv.h"
#include "grainvar.h"
#include "elementnames.h"
#include "gammas.h"
#include "ionbal.h"
#include "ca.h"
#include "mole.h"
#include "phycon.h"
#include "hmi.h"
#include "rfield.h"
#include "atoms.h"
#include "iso.h"
#include "oxy.h"
#include "atmdat.h"
#include "fe.h"

void ion_photo(
	/* nlem is atomic number on C scale, 0 for H */
	long int nelem , 
	/* debugging flag to turn on print */
	bool lgPrintIt )
{
	long int ion, 
	  iphi, 
	  iplow, 
	  ipop, 
	  limit_hi, 
	  limit_lo,
	  ns;

	DEBUG_ENTRY( "ion_photo()" );

	/* IonLow(nelem) and IonHigh(nelem) are bounds for evaluation*/

	/*begin sanity checks */
	ASSERT( nelem < LIMELM );
	ASSERT( dense.IonLow[nelem] >= 0 );
	ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
	ASSERT( dense.IonHigh[nelem] <= nelem + 1);
	/*end sanity checks */

	/* NB - in following, nelem is on c scale, so is 29 for Zn */

	/* min since iso-like atom rates produced in iso_photo.
	 * IonHigh and IonLow range from 0 to nelem+1, but for photo rates
	 * we want 0 to nelem since cannot destroy fully stripped ion.  
	 * since iso-seq done elsewhere, want to actually do IonHigh-xxx.*/
	/* >>chng 00 dec 07, logic on limit_hi now precisely identical to ion_solver */
	/* >>chng 02 mar 31, change limit_hi to < in loop (had been <=) and
	 * also coded for arbitrary number of iso sequences */
	/*limit_hi = MIN2(nelem,dense.IonHigh[nelem]);
	limit_hi = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
	limit_hi = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );

	/* >>chng 03 sep 26, always do atom itself since may be needed for molecules */
	limit_hi = MAX2( 1 , limit_hi );

	/* when grains are present want to use atoms as lower bound to number of stages of ionization,
	 * since atomic rates needed for species within grains */
	if( !conv.nPres2Ioniz && gv.lgDustOn() )
	{
		limit_lo = 0;
	}
	else
	{
		limit_lo = dense.IonLow[nelem];
	}

	/* >>chng 01 dec 11, lower bound now limit_lo */
	/* loop over all ions for this element */
	/* >>chng 02 mar 31, now ion < limit_hi not <= */
	for( ion=limit_lo; ion < limit_hi; ion++ )
	{
		/* loop over all shells for this ion */
		for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
		{
			/* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
			if( (ns==(Heavy.nsShells[nelem][ion]-1) || opac.lgRedoStatic) )
			{
				/* option to redo the rates only on occasion */
				iplow = opac.ipElement[nelem][ion][ns][0];
				iphi = opac.ipElement[nelem][ion][ns][1];
				ipop = opac.ipElement[nelem][ion][ns][2];

				t_phoHeat photoHeat;

				/* compute the photoionization rate, ionbal.lgPhotoIoniz_On is 1, set 0
				 * with "no photoionization" command */
				ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
					GammaK(iplow,iphi,
					ipop,t_yield::Inst().elec_eject_frac(nelem,ion,ns,0),
					&photoHeat )*ionbal.lgPhotoIoniz_On;

				/* these three lines must be kept parallel with the lines
				 * in GammaK ion*/

				/* the heating rate */
				ionbal.PhotoRate_Shell[nelem][ion][ns][1] = photoHeat.HeatLowEnr*ionbal.lgPhotoIoniz_On;
				ionbal.PhotoRate_Shell[nelem][ion][ns][2] = photoHeat.HeatHiEnr*ionbal.lgPhotoIoniz_On;
			}
		}

		/* add on compton recoil ionization for atoms to outer shell */
		/* >>chng 02 mar 24, moved here from ion_solver */
		/* this is the outer shell */
		ns = (Heavy.nsShells[nelem][ion]-1);
		/* this must be moved to photoionize and have code parallel to iso_photo code */
		ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ionbal.CompRecoilIonRate[nelem][ion];
		/* add the heat as secondary-ionization capable heating */
		ionbal.PhotoRate_Shell[nelem][ion][ns][2] += ionbal.CompRecoilHeatRate[nelem][ion];
	}

	/* option to print information about these rates for this element */
	if( lgPrintIt )
	{
		/* option to print rates for particular shell */
		ns = 5;
		ion = 1;
		GammaPrt(
		  opac.ipElement[nelem][ion][ns][0], 
		  opac.ipElement[nelem][ion][ns][1], 
		  opac.ipElement[nelem][ion][ns][2], 
		  ioQQQ, /* io unit we will write to */
		  ionbal.PhotoRate_Shell[nelem][ion][ns][0], 
		  0.05);

		/* outer loop is from K to most number of shells present in atom */
		for( ns=0; ns < Heavy.nsShells[nelem][0]; ns++ )
		{
			fprintf( ioQQQ, "\n %s", elementnames.chElementNameShort[nelem] );
			fprintf( ioQQQ, " %s" , Heavy.chShell[ns]);
			/* MB hydrogenic photo rate may not be included in beow */
			for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
			{
				if( Heavy.nsShells[nelem][ion] > ns )
				{
					fprintf( ioQQQ, " %8.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
				}
				else
				{
					break;
				}
			}
		}
		fprintf(ioQQQ,"\n");
	}
	/* >>chng 11 may 06, moved from ion_calci.cpp */
	/* Ly-alpha photoionization of Ca+
	 * valence shell is reevaluated by ion_photo on every call, so this does not double count */
	if( nelem == ipCALCIUM )
	{
		long ns = 6, ion = 1;
		ionbal.PhotoRate_Shell[nelem][ion][ns][0] += ca.dstCala;
	}
	if( nelem == ipCARBON )
	{
		/* >>chng 05 aug 10, add Leiden hack here to get same C0 photo rate
		 * as UMIST - negates difference in grain opacities */
		if(mole_global.lgLeidenHack)
		{
			int nelem=ipCARBON , ion=0 , ns=2;
			ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
				(HMRATE((1e-10)*3.0,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*
				exp(-(3.0*rfield.extin_mag_V_point))/1.66));
			/* heating rates */
			ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
			ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
		}
	}
	if( nelem == ipNITROGEN )
	{
		/* photoionization from 2D of NI, is atomic nitrogen present? */
		if( dense.xIonDense[ipNITROGEN][0] > 0. )
		{
			t_phoHeat photoHeat;
			// photo rate, population atoms.p2nit evaluated in cooling
			atoms.d5200r = (realnum)GammaK(opac.in1[0],opac.in1[1],opac.in1[2],1.,&photoHeat);
			/* valence shell photoionization, followed by heating; [0][6] => atomic nitrogen */
			ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][0]*
			  (1. - atoms.p2nit) + atoms.p2nit*atoms.d5200r;
			ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1] = ionbal.PhotoRate_Shell[ipNITROGEN][0][2][1]*
			  (1. - atoms.p2nit) + photoHeat.HeatNet*atoms.p2nit;
		}
		else
		{
			atoms.p2nit = 0.;
			atoms.d5200r = 0.;
		}
	}
	if ( nelem == ipMAGNESIUM )
	{
		if( dense.IonLow[ipMAGNESIUM] <= 1 )
		{
			t_phoHeat dummy;
			/* photoionization from excited upper state of 2798 */
			realnum rmg2l = (realnum)GammaK(opac.ipmgex,
					iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,opac.ipOpMgEx,1., &dummy );
					ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] += rmg2l*atoms.popmg2;

			if( nzone <= 1 )
			{
				atoms.xMg2Max = 0.;
			}
			else if( ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0] > 1e-30 )
			{
				/* remember max relative photoionization rate for possible comment */
				atoms.xMg2Max = (realnum)(MAX2(atoms.xMg2Max,rmg2l*atoms.popmg2/
					 ionbal.PhotoRate_Shell[ipMAGNESIUM][1][3][0]));
			}
		}
		else
		{
			atoms.xMg2Max = 0.;
		}
	}

	if ( nelem == ipOXYGEN )
	{

		t_phoHeat dummy;
		/* oxygen, atomic number 8 */
		if( !dense.lgElmtOn[ipOXYGEN] )
		{
			oxy.poiii2Max = 0.;
			oxy.poiii3Max = 0.;
			oxy.r4363Max = 0.;
			oxy.r5007Max = 0.;
			oxy.poiii2 = 0.;
			oxy.AugerO3 = 0.;
			oxy.s3727 = 0.;
			oxy.s7325 = 0.;
			thermal.heating[7][9] = 0.;
			oxy.poimax = 0.;
			return;
		}
		else
		{
			double aeff;

			/* photoionization from O++ 1D
			 *
			 * estimate gamma function by assuming no frequency dependence
			 * betwen 1D and O++3P edge */
			/* destroy upper level of OIII 5007*/
			oxy.d5007r = (realnum)(GammaK(opac.ipo3exc[0],opac.ipo3exc[1],
			  opac.ipo3exc[2] , 1., &dummy ));

			/* destroy upper level of OIII 4363*/
			oxy.d4363 = (realnum)(GammaK(opac.ipo3exc3[0],opac.ipo3exc3[1],
			  opac.ipo3exc3[2] , 1., &dummy ));

			/* destroy upper level of OI 6300*/
			oxy.d6300 = (realnum)(GammaK(opac.ipo1exc[0],opac.ipo1exc[1],
			  opac.ipo1exc[2] , 1., &dummy ));

			/* A21 = 0.0263 */
			aeff = 0.0263 + oxy.d5007r;

			/* 1. as last arg makes this the relative population */
			oxy.poiii2 = (realnum)(atom_pop2(2.5,9.,5.,aeff,2.88e4,1.)/aeff);
			{
				enum {DEBUG_LOC=false};
				if( DEBUG_LOC )
				{
					fprintf(ioQQQ,"pop rel  %.1e rate %.1e  grnd rate %.1e\n",
						oxy.poiii2 , oxy.d5007r ,ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] );
				}
			}

			/* photoionization from excited states */
			if( nzone > 0 )
			{
				/* neutral oxygen destruction */
				ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]*
				  (1. - oxy.poiexc) + oxy.d6300*oxy.poiexc;

				/* doubly ionized oxygen destruction */
				ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] = ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]*
				  (1. - oxy.poiii2 - oxy.poiii3) + oxy.d5007r*oxy.poiii2 +
				  oxy.d4363*oxy.poiii3;

				if( ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > 1e-30 && dense.IonLow[ipOXYGEN] <= 2 )
				{
					if( (oxy.d5007r*oxy.poiii2 + oxy.d4363*oxy.poiii3)/
					  ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0] > (oxy.r4363Max +
					  oxy.r5007Max) )
					{
						oxy.poiii2Max = (realnum)(oxy.d5007r*oxy.poiii2/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
						oxy.poiii3Max = (realnum)(oxy.d4363*oxy.poiii3/ionbal.PhotoRate_Shell[ipOXYGEN][2][2][0]);
					}
					oxy.r4363Max = (realnum)(MAX2(oxy.r4363Max,oxy.d4363));
					oxy.r5007Max = (realnum)(MAX2(oxy.r5007Max,oxy.d5007r));
				}

				/* ct into excited states */
				if( dense.IonLow[ipOXYGEN] <= 0 && (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0] +
				  atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1]) > 1e-30 )
				{
					oxy.poimax = (realnum)(MAX2(oxy.poimax,oxy.d6300*oxy.poiexc/
					  (ionbal.PhotoRate_Shell[ipOXYGEN][0][2][0]+
					  atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]* dense.xIonDense[ipHYDROGEN][1])));
				}
			}
			else
			{
				oxy.poiii2Max = 0.;
				oxy.poiii3Max = 0.;
				oxy.r4363Max = 0.;
				oxy.r5007Max = 0.;
				oxy.poimax = 0.;
			}
		}
		long int iup;
		/* save atomic oxygen photodistruction rate for 3727 creation */
		if( dense.IonLow[ipOXYGEN] == 0 && oxy.i2d < rfield.nflux )
		{
			oxy.s3727 = (realnum)(GammaK(oxy.i2d,oxy.i2p,opac.iopo2d , 1., &dummy ));

			iup = MIN2(iso_sp[ipH_LIKE][1].fb[0].ipIsoLevNIonCon,rfield.nflux);
			oxy.s7325 = (realnum)(GammaK(oxy.i2d,iup,opac.iopo2d , 1., &dummy ));

			oxy.s7325 -= oxy.s3727;
			oxy.s3727 = oxy.s3727 + oxy.s7325;

			/* ratio of cross sections */
			oxy.s7325 *= 0.66f;
		}
		else
		{
			oxy.s3727 = 0.;
			oxy.s7325 = 0.;
		}

		oxy.AugerO3 = (realnum)ionbal.PhotoRate_Shell[ipOXYGEN][0][0][0];

		oxy.s3727 *= dense.xIonDense[ipOXYGEN][0];
		oxy.s7325 *= dense.xIonDense[ipOXYGEN][0];
		oxy.AugerO3 *= dense.xIonDense[ipOXYGEN][0];
	}
	if( nelem == ipIRON )
	{
		if( !dense.lgElmtOn[ipIRON] )
		{
			fe.fekcld = 0.;
			fe.fekhot = 0.;
			fe.fegrain = 0.;
		}
		else
		{
			const int NDIM = ipIRON+1;

			static const double fyield[NDIM+1] = {.34,.34,.35,.35,.36,.37,.37,.38,.39,.40,
			  .41,.42,.43,.44,.45,.46,.47,.47,.48,.48,.49,.49,.11,.75,0.,0.,0.};

			long int i, limit, limit2;
			/* now find total Auger yield of K-alphas
			 * "cold" iron has M-shell electrons, up to Fe 18 */
			fe.fekcld = 0.;
			limit = MIN2(18,dense.IonHigh[ipIRON]);

			for( i=dense.IonLow[ipIRON]; i < limit; i++ )
			{
				ASSERT( i < NDIM + 1 );
				fe.fekcld +=
					(realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
				  fyield[i]);
			}

			/* same sum for hot iron */
			fe.fekhot = 0.;
			limit = MAX2(18,dense.IonLow[ipIRON]);

			limit2 = MIN2(ipIRON+1,dense.IonHigh[ipIRON]);
			ASSERT( limit2 <= LIMELM + 1 );

			for( i=limit; i < limit2; i++ )
			{
				ASSERT( i < NDIM + 1 );
				fe.fekhot +=
					(realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*dense.xIonDense[ipIRON][i]*
				  fyield[i]);
			}

			/* Fe Ka from grains - Fe in grains assumed to be atomic
			 * gv.elmSumAbund[ipIRON] is number density of iron added over all grain species */
			i = 0;
			/* fyield is 0.34 for atomic fe */
			fe.fegrain = ( gv.lgWD01 ) ? 0.f : (realnum)(ionbal.PhotoRate_Shell[ipIRON][i][0][0]*fyield[i]*
						 gv.elmSumAbund[ipIRON]);
		}
	}

	return;
}