/* 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 "transition.h"
#include "version.h"
#include "dense.h"
#include "elementnames.h"
#include "lines.h"
#include "lines_service.h"
#include "opacity.h"
#include "phycon.h"
#include "radius.h"
#include "rfield.h"
#include "rt.h"
#include "taulines.h"
#include "conv.h"

/*outline - adds line photons to reflin and outlin */
/*PutLine enter local line intensity into the intensity stack for eventual printout */
/*PutExtra enter and 'extra' intensity source for some line */
/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
/*TexcLine derive excitation temperature of line from contents of line array */
/*transition::Zero zeros out transition */
/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
/*OccupationNumberLine - derive the photon occupation number at line center for any line */
/*MakeCS compute collision strength by g-bar approximations */
/*gbar1 compute g-bar collision strength using Mewe approximations */
/*gbar0 compute g-bar gaunt factor for neutrals */
/*emit_frac returns fraction of populations the produce emission */
/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
/*chLineLbl use information in line transfer arrays to generate a line label */
/*PutCS enter a collision strength into an individual line vector */

/*outline - adds line photons to reflin and outlin */
void TransitionProxy::outline_resonance( ) const
{
	bool lgDoChecks = true;
	outline(Emis().ColOvTot(), lgDoChecks);
}

/*outline - adds line photons to reflin and outlin */
void TransitionProxy::outline( double nonScatteredFraction, 
										 bool lgDoChecks  ) const
{
	long int ip = ipCont()-1;

	DEBUG_ENTRY( "TransitionProxy::outline()" );

	if ( 0 && lgDoChecks)
	{
		// nothing to do if very small photon flux or below plasma frequency
		ASSERT( Emis().phots() >= 0. );
		if( Emis().phots() < SMALLFLOAT || EnergyErg() / EN1RYD <= rfield.plsfrq )
			return;

		ASSERT( Emis().FracInwd() >= 0. );
		ASSERT( radius.BeamInIn >= 0. );
		ASSERT( radius.BeamInOut >= 0. );
		double Ptot = Emis().Pesc() + Emis().Pelec_esc();
		double PhotEmit = Emis().Aul()*Ptot*(*Hi()).Pop();
		// do not assert accuracy if close to bounds of fp precision
		double error = MAX2( SMALLFLOAT*1e3 , 3e-1*PhotEmit );
		// do not assert if do not have valid solution
		bool lgGoodSolution = conv.lgConvEden && conv.lgConvIoniz() &&
			conv.lgConvPops && conv.lgConvPres && conv.lgConvTemp;
		// see ticket #135 rt inconsistent results
		// this assert trips on a regular basis.  The fix is to rewrite
		// the RT and level population solvers so that the level population
		// and OTS rates are done simultaneously rather than sequentially
		// For now do not throw the assert on a release version - we know
		// about this problem
		ASSERT( t_version::Inst().lgRelease || !lgGoodSolution || 
			Ptot < 1e-3 ||  fp_equal_tol( Emis().phots(), PhotEmit , error ));
	}

	bool lgTransStackLine = true;
	outline_base(Emis().dampXvel(), Emis().damp(), lgTransStackLine, ip, Emis().phots(), Emis().FracInwd(),
					 nonScatteredFraction);
}

/*emit_frac returns fraction of populations the produce emission */
double emit_frac(const TransitionProxy& t)
{
	DEBUG_ENTRY( "emit_frac()" );

	ASSERT( t.ipCont() > 0 );

	/* collisional deexcitation and destruction by background opacities
	 * are loss of photons without net emission */
	double deexcit_loss = t.Coll().col_str() * dense.cdsqte + t.Emis().Aul()*t.Emis().Pdest();
	/* this is what is observed */
	double rad_deexcit = t.Emis().Aul()*(t.Emis().Pelec_esc() + t.Emis().Pesc());
	return rad_deexcit/(deexcit_loss + rad_deexcit);
}

/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
void DumpLine(const TransitionProxy& t)
{
	char chLbl[110];

	DEBUG_ENTRY( "DumpLine()" );

	ASSERT( t.ipCont() > 0 );

	/* routine to print contents of line arrays */
	strcpy( chLbl, "DEBUG ");
	strcat( chLbl, chLineLbl(t) );

	fprintf( ioQQQ, 
		"%10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n", 
	  chLbl, 
	  phycon.te, 
	  dense.eden, 
	  t.Coll().col_str(), 
	  t.Emis().Aul(), 
	  TexcLine(t), 
	  t.Coll().cool(), 
	  t.Coll().heat() ,
	  opac.opacity_abs[t.ipCont()-1],
	  opac.albedo[t.ipCont()-1]);

	fprintf( ioQQQ, 
		"Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n", 
	  t.Emis().TauIn(), 
	  t.Emis().TauTot(), 
	  t.Emis().Pesc(), 
	  t.Emis().Pelec_esc(), 
	  t.Emis().Pdest(), 
	  t.Emis().pump(), 
	  t.Emis().ots(), 
	  (*t.Lo()).Pop(), 
	  (*t.Hi()).Pop() ,
	  t.Emis().PopOpc() );
	return;
}


/*OccupationNumberLine - derive the photon occupation number at line center for any line */
double OccupationNumberLine(const TransitionProxy& t)
{
	double OccupationNumberLine_v;

	DEBUG_ENTRY( "OccupationNumberLine()" );

	ASSERT( t.ipCont() > 0 );

	/* routine to evaluate line photon occupation number - 
	 * return negative number if line is maser */
	if( fabs(t.Emis().PopOpc()) > SMALLFLOAT )
	{
		/* the lower population with correction for stimulated emission */
		OccupationNumberLine_v = ( (*t.Hi()).Pop() / (*t.Hi()).g() ) /
			( t.Emis().PopOpc() / (*t.Lo()).g() )  *
			(1. - t.Emis().Pesc());
	}
	else
	{
		OccupationNumberLine_v = 0.;
	}
	/* return value is not guaranteed to be positive - negative if
	 * line mases */
	return( OccupationNumberLine_v );
}

/*TexcLine derive excitation temperature of line from contents of line array */
double TexcLine(const TransitionProxy& t)
{
	double TexcLine_v;

	DEBUG_ENTRY( "TexcLine()" );

	/* routine to evaluate line excitation temp using contents of line array
	 * */
	if( (*t.Hi()).Pop() * (*t.Lo()).Pop() > 0. )
	{
		TexcLine_v = ( (*t.Hi()).Pop() / (*t.Hi()).g() )/( (*t.Lo()).Pop() / (*t.Lo()).g() );
		TexcLine_v = log(TexcLine_v);
		/* protect against infinite temp limit */
		if( fabs(TexcLine_v) > SMALLFLOAT )
		{
			TexcLine_v = - t.EnergyK() / TexcLine_v;
		}
	}
	else
	{
		TexcLine_v = 0.;
	}
	return( TexcLine_v );
}

/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
void chIonLbl(char *chIonLbl_v, const TransitionProxy& t)
{
	DEBUG_ENTRY( "chIonLbl()" );

	/* function to use information within the line array
	 * to generate an ion label, giving element and
	 * ionization stage
	 * */
	if( (*t.Hi()).nelem() < 0 )
	{
		/* this line is to be ignored */
		if( (*t.Hi()).chLabel()[0]=='\0' )
			strcpy( chIonLbl_v, "Dumy" );
		else 
			strcpy( chIonLbl_v, (*t.Hi()).chLabel() );
	}
	else
	{
		chIonLbl( chIonLbl_v, (*t.Hi()).nelem(), (*t.Hi()).IonStg() );
	}
	/* chIonLbl is four char null terminated string */
	return/*( chIonLbl_v )*/;
}

void chIonLbl(char *chIonLbl_v, const long& nelem, const long& IonStg)
{
	DEBUG_ENTRY( "chIonLbl()" );

	// NB - nelem passed here is expected to be on the physical scale, not C (so hydrogen = 1).
	ASSERT( nelem >= 1 );
	ASSERT( nelem <= LIMELM );
	/* ElmntSym.chElementSym is null terminated, 2 ch + null, string giving very
	 * short form of element name */
	strcpy( chIonLbl_v , elementnames.chElementSym[nelem-1] );
	/* chIonStage is four char null terminated string, starting with "_1__" */
	strcat( chIonLbl_v, elementnames.chIonStage[IonStg-1] );
	return;
}

/*chLineLbl use information in line transfer arrays to generate a line label */
/* this label is null terminated */
/* ContCreatePointers has test this with full range of wavelengths */
char* chLineLbl(const TransitionProxy& t)
{
	static char chLineLbl_v[11];
	static char chSpecies[5];

	DEBUG_ENTRY( "chLineLbl()" );

	if( (*t.Hi()).nelem() < 1 && (*t.Hi()).IonStg() < 1 )
	{
		sprintf( chSpecies, "%4.4s", (*t.Hi()).chLabel() );
	}
	else
	{
		ASSERT( (*t.Hi()).nelem() >= 1 );
		ASSERT( (*t.Hi()).IonStg() >= 1 && (*t.Hi()).IonStg() <= (*t.Hi()).nelem() + 1 );
		sprintf( chSpecies, "%2.2s%2.2s",
			elementnames.chElementSym[(*t.Hi()).nelem() -1], 
			elementnames.chIonStage[(*t.Hi()).IonStg()-1] ); 
	}

	/* function to use information within the line array
	 * to generate a line label, giving element and
	 * ionization stage
	 * rhs are set in large block data */

	/* NB this function is profoundly slow due to sprintf statement
	 * also - it cannot be evaluated within a write statement itself*/

	if( t.WLAng() > (realnum)INT_MAX )
	{
		sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies, 
		   (int)(t.WLAng()/1e8), 'c' );
	}
	else if( t.WLAng() > 9999999. )
	{
		/* wavelength is very large, convert to centimeters */
		sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies, 
			t.WLAng()/1e8, 'c' );
	}
	else if( t.WLAng() > 999999. )
	{
		/* wavelength is very large, convert to microns */
		sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies, 
			(int)(t.WLAng()/1e4), 'm' );
	}
	else if( t.WLAng() > 99999. )
	{
		/* wavelength is very large, convert to microns */
		sprintf( chLineLbl_v, "%4.4s%5.1f%c", chSpecies, 
			t.WLAng()/1e4, 'm' );
	}
	else if( t.WLAng() > 9999. )
	{
		sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies, 
		   t.WLAng()/1e4, 'm' );
	}
	else if( t.WLAng() >= 100. )
	{
		sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies, 
		   (int)t.WLAng(), 'A' );
	}
	/* the following two formats should be changed for the
	 * wavelength to get more precision */
	else if( t.WLAng() >= 10. )
	{
		sprintf( chLineLbl_v, "%4.4s%5.1f%c", chSpecies, 
		   t.WLAng(), 'A' );
	}
	else
	{
		sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies, 
		   t.WLAng(), 'A' );
	}
	/* make sure that string ends with null character - this should
	 * be redundant */
	ASSERT( chLineLbl_v[10]=='\0' );
	return( chLineLbl_v );
}

/*PutCS enter a collision strength into an individual line vector */
void PutCS(double cs, 
  const TransitionProxy& t)
{
	DEBUG_ENTRY( "PutCS()" );

	/* collision strength must be non-negative */
	ASSERT( cs > 0. );

	t.Coll().col_str() = (realnum)cs;

	return;
}

void GenerateTransitionConfiguration( const TransitionProxy &t, char *chComment )
{
	strcpy( chComment, (*t.Lo()).chConfig() );
	strcat( chComment, " - " );
	strcat( chComment, (*t.Hi()).chConfig() );
	return;
}

static realnum ExtraInten;

/*PutLine enter local line intensity into the intensity stack for eventual printout */
STATIC void PutLine_base(const TransitionProxy& t, const char *chComment, const char *chLabelTemp, bool lgLabel)
{
	DEBUG_ENTRY( "PutLine_base()" );

	char chLabel[5];
	double xIntensity,
		other,
		xIntensity_in;
		
	/* routine to use line array data to generate input
	 * for emission line array */
	ASSERT( t.ipCont() > 0. );
		
	if (lgLabel == true)
	{
		strncpy( chLabel, chLabelTemp, 4 );
		chLabel[4] = '\0';
	}

	/* if ipass=0 then we must generate label info since first pass
	 * gt.0 then only need line intensity data */
	if( LineSave.ipass == 0 )
	{
		if (lgLabel == false)
		{
			/* these variables not used by linadd if ipass ne 0 */
			/* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" */
			chIonLbl(chLabel, t);
		}
		xIntensity = 0.;
	}
	else
	{
		/* both the counting and integrating modes comes here */
		/* not actually used so set to safe value */
		chLabel[0] = '\0';
		/* total line intensity or luminosity 
		 * these may not be defined in initial calls so define here */
		xIntensity = t.Emis().xIntensity() + ExtraInten;
	}

	/* initial counting case, where ipass == -1, just ignored above, call linadd below */
	
	/* ExtraInten is option that allows extra intensity (i.e., recomb)
	 * to be added to this line  with Call PutExtra( exta )
	 * in main lines this extra
	 * contribution must be identified explicitly */
	ExtraInten = 0.;
	/*linadd(xIntensity,wl,chLabel,'i');*/
	/*lindst add line with destruction and outward */
	rt.fracin = t.Emis().FracInwd();
	lindst(xIntensity, 
			 t.WLAng(), 
			 chLabel, 
			 t.ipCont(), 
			 /* this is information only - has been counted in cooling already */
			 'i', 
			 /* do not add to outward beam, also done separately */
			 false,
			 chComment);
	rt.fracin = 0.5;
		
	/* inward part of line - do not move this away from previous lines
	 * since xIntensity is used here */
	xIntensity_in = xIntensity*t.Emis().FracInwd();
	ASSERT( xIntensity_in>=0. );
	linadd(xIntensity_in,t.WLAng(),"Inwd",'i',chComment);
	
	/* cooling part of line */
	other = t.Coll().cool();
	linadd(other,t.WLAng(),"Coll",'i',chComment);
	
	/* fluorescent excited part of line */
	double radiative_branching;
	enum { lgNEW = true };
	if (lgNEW)
	{
		// Improved two-level version of radiative branching ratio
		const double AulEscp = t.Emis().Aul()*(t.Emis().Pesc()+t.Emis().Pelec_esc());
		// Would be better to include all outward transition processes from the
		// line, to cater for the general non-two-level case
		const double sinkrate = AulEscp + t.Emis().Aul()*t.Emis().Pdest() + t.Coll().ColUL( colliders );
		if (sinkrate > 0.0) 
		{
			radiative_branching = AulEscp/sinkrate;
		}
		else
		{
			radiative_branching = 0.;
		}
	}
	else
	{
		// This is the excitation ratio not the de-excitation ratio according
		// to its specification
		radiative_branching = (1.-t.Emis().ColOvTot());
	}

	other = (*t.Lo()).Pop() * t.Emis().pump() * radiative_branching * t.EnergyErg();
	linadd(other,t.WLAng(),"Pump",'i',chComment);
		

	/* heating part of line */
	other = t.Coll().heat();
	linadd(other,t.WLAng(),"Heat",'i',chComment);
}

/*PutLine enter local line intensity into the intensity stack for eventual printout */
void PutLine(const TransitionProxy& t, const char *chComment, const char *chLabelTemp)
{	
	const bool lgLabel = true;
	DEBUG_ENTRY( "PutLine()" );
	PutLine_base(t, chComment, chLabelTemp, lgLabel);
}

/*PutLine enter local line intensity into the intensity stack for eventual printout */
void PutLine(const TransitionProxy& t,
	     const char *chComment)
{
	const bool lgLabel = false;
	const char *chLabelTemp = NULL;
	DEBUG_ENTRY( "PutLine()" );
	PutLine_base(t, chComment, chLabelTemp, lgLabel);
}

/* ==================================================================== */
/*PutExtra enter and 'extra' intensity source for some line */
void PutExtra(double Extra)
{

	DEBUG_ENTRY( "PutExtra()" );

	ExtraInten = (realnum)Extra;
	return;
}

void TransitionProxy::Junk( void ) const
{

	DEBUG_ENTRY( "TransitionProxy::Junk()" );

		/* wavelength, usually in A, used for printout */
	WLAng() = -FLT_MAX;

	/* transition energy in wavenumbers */
	EnergyWN() = -FLT_MAX;

	/* array offset for radiative transition within continuum array 
	 * is negative if transition is non-radiative. */
	ipCont() = -10000;

	CollisionJunk( Coll() );

	/* set these equal to NULL first. That will cause the code to crash if
	 * the variables are ever used before being deliberately set. */
	ipEmis() = -1;
	
	setLo(-1);
	setHi(-1);
	return;
}

/*TransitionZero zeros out transition array at start of calculation, sets
 * optical depths to initial values */
void TransitionProxy::Zero( void ) const
{

	DEBUG_ENTRY( "TransitionProxy::Zero()" );

	CollisionZero( Coll() );

	::Zero( *Lo() );
	::Zero( *Hi() );
	EmLineZero( Emis() );
	TauZero( Emis() );

	return;
}

/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
void LineConvRate2CS( const TransitionProxy& t , realnum rate )
{

	DEBUG_ENTRY( "LineConvRate2CS()" );

	/* return is collision strength, convert from collision rate from 
	 * upper to lower, this assumes pure electron collisions, but that will
	 * also be assumed by anything that uses cs, for self-consistency */
	t.Coll().col_str() = rate * (*t.Hi()).g() / (realnum)dense.cdsqte;

	/* change assert to non-negative - there can be cases (Iin H2) where cs has
	 * underflowed to 0 on some platforms */
	ASSERT( t.Coll().col_str() >= 0. );
	return;
}

/*gbar0 compute g-bar gaunt factor for neutrals */
STATIC void gbar0(double ex, 
	  realnum *g)
{
	double a, 
	  b, 
	  c, 
	  d, 
	  y;

	DEBUG_ENTRY( "gbar0()" );

	/* written by Dima Verner
	 *
	 * Calculation of the effective Gaunt-factor by use of 
	 * >>refer	gbar	cs	Van Regemorter, H., 1962, ApJ 136, 906
	 * fits for neutrals
	 *  Input parameters: 
	 * ex - energy ryd - now K
	 * t  - temperature in K
	 *  Output parameter:
	 * g  - effective Gaunt factor
	 * */

	/* y = ex*157813.7/te */
	y = ex/phycon.te;
	if( y < 0.01 )
	{
		*g = (realnum)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
	}
	else
	{
		if( y > 10.0 )
		{
			*g = (realnum)(0.066/sqrt(y));
		}
		else
		{
			a = 1.5819068e-02;
			b = 1.3018207e00;
			c = 2.6896230e-03;
			d = 2.5486007e00;
			d = log(y/c)/d;
			*g = (realnum)(a + b*exp(-0.5*POW2(d)));
		}
	}
	return;
}

/*gbar1 compute g-bar collision strength using Mewe approximations */
STATIC void gbar1(double ex, 
	  realnum *g)
{
	double y;

	DEBUG_ENTRY( "gbar1()" );

	/*	*written by Dima Verner
	 *
	 * Calculation of the effective Gaunt-factor by use of 
	 * >>refer	gbar	cs	Mewe,R., 1972, A&A 20, 215
	 * fits for permitted transitions in ions MgII, CaII, FeII (delta n = 0)
	 * Input parameters: 
	 * ex - excitation energy in Ryd - now K
	 * te  - temperature in K
	 * Output parameter:
	 * g  - effective Gaunt factor
	 */

	/* y = ex*157813.7/te */
	y = ex/phycon.te;
	*g = (realnum)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
	return;
}

/*MakeCS compute collision strength by g-bar approximations */
void MakeCS(const TransitionProxy& t)
{
	long int ion;
	double Abun, 
	  cs;
	realnum
	  gbar;

	DEBUG_ENTRY( "MakeCS()" );

	/* routine to get cs from various approximations */

	/* check if abundance greater than 0 */
	ion = (*t.Hi()).IonStg();

	//This is the oscillator strength limit where larger values are assumed to be for allowed transitions.
	static double gfLimit = 1e-8;

	Abun = dense.xIonDense[ (*t.Hi()).nelem() -1 ][ ion-1 ];
	if( Abun <= 0. )
	{
		gbar = 1.;
	}
	else if( t.Emis().gf() >= gfLimit )
	{
		/* check if neutral or ion */
		if( ion == 1 )
		{
			/* neutral - compute gbar for eventual cs */
			gbar0(t.EnergyK(),&gbar);
		}
		else
		{
			/* ion - compute gbar for eventual cs */
			gbar1(t.EnergyK(),&gbar);
		}
	}
	else
	{
		//Mewe72 provides a gbar estimate for forbidden transitions.
		gbar = 0.15;
	}

	/* above was g-bar, convert to cs */
	cs = gbar*(14.5104/WAVNRYD)*t.Emis().gf()/t.EnergyWN();

	/* stuff the cs in place */
	t.Coll().col_str() = (realnum)cs;
	return;
}

void TransitionProxy::AddLine2Stack() const
{
	DEBUG_ENTRY( "AddLine2Stack()" );

	ASSERT( lgLinesAdded == false );

	size_t newsize = m_list->Emis.size()+1; 
	m_list->Emis.resize(newsize);
	ipEmis() = newsize-1;
	this->resetEmis();
}

void TransitionProxy::AddLoState() const
{
	DEBUG_ENTRY( "AddLoState()" );

	ASSERT( !lgStatesAdded );

	m_list->states->resize(m_list->states->size()+1);

	setLo(m_list->states->size()-1);
}

void TransitionProxy::AddHiState() const
{
	DEBUG_ENTRY( "AddHiState()" );

	ASSERT( !lgStatesAdded );

	m_list->states->resize(m_list->states->size()+1);

	setHi(m_list->states->size()-1);
}