/* 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 */
#include "cddefines.h"
#include "atmdat.h"
#include "conv.h"
#include "dense.h"
#include "heavy.h"
#include "helike_cs.h"
#include "hydroeinsta.h"
#include "hydrogenic.h"
#include "hydro_vs_rates.h"
#include "ionbal.h"
#include "iso.h"
#include "opacity.h"
#include "phycon.h"
#include "physconst.h"
#include "rfield.h"
#include "secondaries.h"
#include "trace.h"
#include "taulines.h"

/* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};

void iso_collisional_ionization( long ipISO, long nelem )
{
	ASSERT( ipISO < NISO );

	DEBUG_ENTRY( "iso_collisional_ionization()" );
	
	t_iso_sp* sp = &iso_sp[ipISO][nelem];

	/* collisional ionization from ground */
	sp->fb[0].ColIoniz = iso_ctrl.lgColl_ionize[ipISO] *
		t_ADfA::Inst().coll_ion_wrapper( nelem, nelem-ipISO, phycon.te );

	iso_put_error(ipISO,nelem,sp->numLevels_max,0,IPCOLLIS,0.20f,0.20f);
	
	for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
	{
		if( nelem == ipISO )
		{
			/* use routine from Vriens and Smeets (1981). */
			/* >>refer	iso neutral	col.ion.	Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */
			sp->fb[ipHi].ColIoniz = hydro_vs_ioniz( sp->fb[ipHi].xIsoLevNIonRyd, phycon.te );
		}
		else
		{
			/* ions */
			/* use hydrogenic ionization rates for ions
			 * >>refer	iso ions	col.ion.	Allen 1973, Astro. Quan. for low Te.
			 * >>refer	iso ions	col.ion.	Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
			 * */
			sp->fb[ipHi].ColIoniz = 
				Hion_coll_ioniz_ratecoef( ipISO, nelem, N_(ipHi), sp->fb[ipHi].xIsoLevNIonRyd, phycon.te  );
		}
		
		// iso_ctrl.lgColl_ionize is option to turn off collisions, "atom XX-like collis off" command
		sp->fb[ipHi].ColIoniz *= iso_ctrl.lgColl_ionize[ipISO];
 	
		iso_put_error(ipISO,nelem,sp->numLevels_max,ipHi,IPCOLLIS,0.20f,0.20f);
	}
	
	/* Here we arbitrarily scale the highest level ionization to account for the fact
	 * that, if the atom is not full size, this level should be interacting with higher
	 * levels and not just the continuum.  We did add on collisional excitation terms instead
	 * but this caused a problem at low temperatures because the collisional ionization was 
	 * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann
	 * factor.  The result was a collisional recombination that had residual exponentials of
	 * the form exp(x/kT), which blew up at small T.	*/
	if( 0 && !sp->lgLevelsLowered )
	{
		sp->fb[sp->numLevels_max-1].ColIoniz *= 100.;
	}
	
	return;
}

void iso_suprathermal( long ipISO, long nelem )
{
	DEBUG_ENTRY( "iso_suprathermal()" );

	/* check that we were called with valid parameters */
	ASSERT( ipISO < NISO );
	ASSERT( nelem >= ipISO );
	ASSERT( nelem < LIMELM );

	t_iso_sp* sp = &iso_sp[ipISO][nelem];

	/***********************************************************************
	 *                                                                     *
	 * get secondary excitation by suprathermal electrons                  *
	 *                                                                     *
	 ***********************************************************************/

	for( long i=1; i < sp->numLevels_max; i++ )
	{
		if( sp->trans(i,0).ipCont() > 0 )
		{
			/* get secondaries for all iso lines by scaling LyA 
			 * excitation by ratio of cross section (oscillator strength/energy) 
			 * Born approximation or plane-wave approximation based on
			 *>>refer	HI	excitation	Shemansky, D.E., et al., 1985, ApJ, 296, 774 */
			sp->trans(i,0).Coll().rate_lu_nontherm_set() = secondaries.x12tot *
				(sp->trans(i,0).Emis().gf()/
				sp->trans(i,0).EnergyWN()) /
				(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).Emis().gf()/
				iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).EnergyWN()) *
				iso_ctrl.lgColl_excite[ipISO];;
		}
		else
			sp->trans(i,0).Coll().rate_lu_nontherm_set() = 0.;
	}

	return;
}

/*============================*/
/* evaluate collisional rates */
void iso_collide( long ipISO, long nelem )
{
	DEBUG_ENTRY( "iso_collide()" );

	/* this array stores the last temperature at which collision data were evaluated for
	 * each species of the isoelectronic sequence. */
	static double TeUsed[NISO][LIMELM]={
		{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
		{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };

	/* check that we were called with valid parameters */
	ASSERT( ipISO < NISO );
	ASSERT( nelem >= ipISO );
	ASSERT( nelem < LIMELM );
	
	t_iso_sp* sp = &iso_sp[ipISO][nelem];

	/* skip most of this routine if temperature has not changed,
	 * the check on conv.nTotalIoniz is to make sure that we redo this
	 * on the very first call in a grid calc - it is 0 on the first call */
	if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz )
	{
		ASSERT( sp->trans( iso_ctrl.nLyaLevel[ipISO], 0 ).Coll().ColUL( colliders ) >= 0. );

		if( trace.lgTrace  && (trace.lgHBug||trace.lgHeBug) )
		{
			fprintf( ioQQQ, 
				"     iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n",
				iso_ctrl.chISO[ipISO], nelem );
		}
	}
	else
	{
		TeUsed[ipISO][nelem] = phycon.te;

		if( trace.lgTrace  && (trace.lgHBug||trace.lgHeBug) )
		{
			fprintf( ioQQQ, 
				"     iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n",
				iso_ctrl.chISO[ipISO], nelem );
		}
		
		/**********************************************************
		 *                                                        *
		 * Boltzmann factors for all levels,  and                 *
		 * collisional ionization and excitation                  *
		 *                                                        *
		 **********************************************************/

		/* HION_LTE_POP	is planck^2 / (2 pi m_e k ), raised to 3/2 next */
		double factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
			(dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);

		/* term in () is stat weight of electron * ion */
		double ConvLTEPOP = pow(factor,1.5)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;

		sp->lgPopLTE_OK = true;

		// this is the maximum value of iso.PopLTE (units cm^3) that corresponds
		// to the minimum positive density values.  A smaller density will be
		// regarded as zero, and the product PopLTE*n_e*n_Z+ will also be zero.
#define MAX_POP_LTE	(MAX_DENSITY/dense.density_low_limit/dense.density_low_limit)

		/* fully define Boltzmann factors to continuum for model levels */
		for( long ipLo=0; ipLo<sp->numLevels_max; ipLo++ )
		{
			/* this Boltzmann factor is exp( +ioniz energy / Te ) */
			sp->st[ipLo].Boltzmann() = 
				dsexp(sp->st[ipLo].energy().Ryd()/phycon.te_ryd);
			sp->st[ipLo].ConBoltz() =
				dsexp(sp->fb[ipLo].xIsoLevNIonRyd/phycon.te_ryd);

			/***************************************
			 *                                     *
			 * LTE abundances for all levels       *
			 *                                     *
			 ***************************************/

			if( sp->st[ipLo].ConBoltz() > SMALLDOUBLE )
			{
				/* LTE population of given level. */
				sp->fb[ipLo].PopLTE = 
					sp->st[ipLo].g() / sp->st[ipLo].ConBoltz() * ConvLTEPOP;
				ASSERT( sp->fb[ipLo].PopLTE < BIGDOUBLE );
			}
			else
			{
				sp->fb[ipLo].PopLTE = 0.;
			}

			sp->fb[ipLo].PopLTE = MIN2( sp->fb[ipLo].PopLTE, MAX_POP_LTE );

			/* now check for any zeros - if present then matrix cannot be used */
			if( sp->fb[ipLo].PopLTE <= 0. )
			{
				sp->lgPopLTE_OK = false;
			}
		}

		iso_collisional_ionization( ipISO, nelem );

		/***********************************************************
		 *                                                         *
		 * collisional deexcitation for all lines in iso sequence  *
		 *                                                         *
		 ***********************************************************/

		if( iso_ctrl.lgColl_excite[ipISO] )
		{
			for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
			{
				for( long ipLo=0; ipLo < ipHi; ipLo++ )
				{
					for( long ipCollider = ipELECTRON; ipCollider <= ipALPHA; ipCollider++ )
					{
						double cs_temp = 0.;

						if( N_(ipHi) == N_(ipLo) && !iso_ctrl.lgColl_l_mixing[ipISO] )
							cs_temp = 0.;
						else if( N_(ipHi)-N_(ipLo) > 2 && ipCollider > ipELECTRON )
							cs_temp = 0.;
						else if( ipISO == ipH_LIKE )
							cs_temp =  HydroCSInterp( nelem , ipHi , ipLo, ipCollider );
						else if( ipISO == ipHE_LIKE )
							cs_temp = HeCSInterp( nelem , ipHi , ipLo, ipCollider );
						else
							TotalInsanity();

						cs_temp *=  iso_ctrl.lgColl_excite[ipISO];
					
						if( opac.lgCaseB_HummerStorey && N_(ipHi)==N_(ipLo)+1 && abs(L_(ipHi)-L_(ipLo))==1 )
						{
							double Aul = HydroEinstA( N_(ipLo), N_(ipHi) );
							cs_temp *= ( sp->trans(ipHi,ipLo).Emis().gf() / sp->st[ipLo].g() ) /
								( GetGF(Aul, sp->trans(ipHi,ipLo).EnergyWN(), 2.*N_(ipHi)*N_(ipHi))/(2.*N_(ipLo)*N_(ipLo)) );
						}
	
						// store electron collision strength in generic collision strength
						if( ipCollider == ipELECTRON )
							sp->trans(ipHi,ipLo).Coll().col_str() = (realnum) cs_temp;

						double reduced_mass_collider_system = dense.AtomicWeight[nelem]*ColliderMass[ipCollider]/
							(dense.AtomicWeight[nelem]+ColliderMass[ipCollider])*ATOMIC_MASS_UNIT;

						if( ipCollider == ipELECTRON )
							reduced_mass_collider_system  = ELECTRON_MASS;

						double rateCoef = cs_temp *
								pow(ELECTRON_MASS/reduced_mass_collider_system, 1.5) * COLL_CONST/
								( phycon.sqrte * (double)sp->st[ipHi].g() );
					
						if( !opac.lgCaseB_HummerStorey )
						{	
							if( ipISO == ipH_LIKE )
							{
								if( N_(ipHi) > sp->n_HighestResolved_max &&
									N_(ipLo) <= sp->n_HighestResolved_max )
								{
									rateCoef *= (8./3.)*(log(double(N_(ipHi)))+2.);
								}
							}
							else if( ipISO == ipHE_LIKE )
							{
								fixit();
								/* This is intended to be a trick to get the correct collisional excitation from
								 * collapsed levels to resolved levels.  Is it needed or does the stat weight used
								 * above handle it automatically? If it is needed, is this correct? */
								if( N_(ipHi) > sp->n_HighestResolved_max &&
									N_(ipLo) <= sp->n_HighestResolved_max )
								{
									rateCoef *= (2./3.)*(log(double(N_(ipHi)))+2.);
								}
							}
						}
	
						sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[ipCollider] =
								(realnum) rateCoef;
					}
					
					/* check for sanity */
					ASSERT( sp->trans(ipHi,ipLo).Coll().rate_coef_ul()[ipELECTRON] >= 0. );

					if( N_(ipHi) <= 5 && N_(ipLo) <= 2 )
						iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.10f, 0.30f );
					else
						iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.20f, 0.30f );
				}
			}
		}
		
		if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
		{
			fprintf( ioQQQ, "     iso_collide: %s Z=%li de-excitation rates coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
			long upper_limit = sp->numLevels_local;
			for( long ipHi=1; ipHi < upper_limit; ipHi++ )
			{
				fprintf( ioQQQ, " %li\t", ipHi );
				for( long ipLo=0; ipLo < ipHi; ipLo++ )
				{
					fprintf( ioQQQ,PrintEfmt("%10.3e", 
							sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) / dense.EdenHCorr ));
				}
				fprintf( ioQQQ, "\n" );
			}

			fprintf( ioQQQ, "     iso_collide: %s Z=%li collisional ionization coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
			for( long ipHi=0; ipHi < upper_limit; ipHi++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e",  sp->fb[ipHi].ColIoniz ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "     iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
			for( long ipHi=0; ipHi < upper_limit; ipHi++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e",  sp->st[ipHi].ConBoltz() ));
			}
			fprintf( ioQQQ, "\n" );

			fprintf( ioQQQ, "     iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
			for( long ipHi=0; ipHi < upper_limit; ipHi++ )
			{
				fprintf( ioQQQ,PrintEfmt("%10.3e",  sp->fb[ipHi].PopLTE ));
			}
			fprintf( ioQQQ, "\n" );
		}

		/* the case b hummer and storey option,
		 * this kills collisional excitation and ionization from n=1 and n=2 */
		if( opac.lgCaseB_HummerStorey )
		{
			for( long ipLo=0; ipLo<sp->numLevels_max-1; ipLo++ )
			{
				if( N_(ipLo)>=3 )
					break;
				
				sp->fb[ipLo].ColIoniz = 0.;

				for( long ipHi=ipLo+1; ipHi<sp->numLevels_max; ipHi++ )
				{
					/* don't disable 2-2 collisions */
					if( N_(ipLo)==2 && N_(ipHi)==2 )
						continue;

					sp->trans(ipHi,ipLo).Coll().col_str() = 0.;
					for( long k=0; k<ipNCOLLIDER; ++k )
						sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[k] = 0.;
				}
			}
		}
	}

	iso_suprathermal( ipISO, nelem );

	/* this must be reevaluated every time since eden can change when Te does not */
	/* save into main array - collisional ionization by thermal electrons */
	ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] = 
		sp->fb[0].ColIoniz*dense.EdenHCorr;

	/* cooling due to collisional ionization, which only includes thermal electrons */
	ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] = 
		ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]*
		rfield.anu[Heavy.ipHeavy[nelem][nelem-ipISO]-1]*EN1RYD;

	return;
}