/* 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 "ipoint.h"
#include "rfield.h"
#include "two_photon.h"

void TwoPhotonSetup( vector<two_photon> &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem )
{
	DEBUG_ENTRY( "TwoPhotonSetup()" );

	tnu_vec.resize( tnu_vec.size() + 1 );
	two_photon &tnu = tnu_vec.back();

	tnu.ipHi = ipHi;
	tnu.ipLo = ipLo;
	tnu.AulTotal = Aul;
	tnu.Pop = &(*tr.Hi()).Pop();
	tnu.E2nu = tr.EnergyRyd();
	tnu.induc_dn_max = 0.;
	
	/* pointer to the Lya energy */
	tnu.ipTwoPhoE = ipoint(tnu.E2nu);
	while( rfield.anu[tnu.ipTwoPhoE] > tnu.E2nu )
	{
		--tnu.ipTwoPhoE;
	}
	tnu.ipSym2nu.resize( tnu.ipTwoPhoE );
	tnu.As2nu.resize( tnu.ipTwoPhoE );
	tnu.local_emis.resize( tnu.ipTwoPhoE );

	/* >>chng 02 aug 14, change upper limit to full Lya energy */
	for( long i=0; i < tnu.ipTwoPhoE; i++ )
	{
		/* energy is symmetric energy, the other side of half E2nu */
		double energy = tnu.E2nu - rfield.anu[i];
		/* this is needed since mirror image of cell next to two-nu energy
		* may be slightly negative */
		energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
		/* find index for this symmetric energy */
		tnu.ipSym2nu[i] = ipoint(energy);
		while( rfield.anu[tnu.ipSym2nu[i]] > tnu.E2nu ||
			tnu.ipSym2nu[i] >= tnu.ipTwoPhoE)
		{
			--tnu.ipSym2nu[i];
		}
		ASSERT( tnu.ipSym2nu[i] >= 0 );
	}

	double SumShapeFunction = 0., Renorm= 0.;

	/* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
	 * to include that in the following sum */
	for( long i=0; i < tnu.ipTwoPhoE; i++ )
	{
		double ShapeFunction;

		ASSERT( rfield.anu[i]<=tnu.E2nu );
		ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/tnu.E2nu, ipISO, nelem )*rfield.widflx[i]/tnu.E2nu;
		SumShapeFunction += ShapeFunction;

		/* >>refer	HI	2nu	Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
		/* As2nu will add up to the A, so its units are s-1	*/ 
		tnu.As2nu[i] = (realnum)( tnu.AulTotal * ShapeFunction );
	}

	/* The spline function in twophoton.c causes a bit of an error in the integral of the
	 * shape function.  So we renormalize the integral to 1.	*/
	Renorm = 1./SumShapeFunction;

	for( long i=0; i < tnu.ipTwoPhoE; i++ )
	{
		tnu.As2nu[i] *= (realnum)Renorm;
	}

	/* The result should be VERY close to 1.	*/
	ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );

	return;
}

void CalcTwoPhotonRates( two_photon& tnu, bool lgDoInduced )
{
	DEBUG_ENTRY( "CalcTwoPhotonRates()" );
			
	/* this could fail when pops very low and Pop2Ion is zero */
	ASSERT( tnu.ipTwoPhoE>0 && tnu.E2nu>0. );

	double sum = 0.;
	tnu.induc_up = 0.;
	tnu.induc_dn = 0.;
	/* two photon emission, ipTwoPhoE is 
	* continuum array index for Lya energy */
	for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
	{
		// We do not assert rfield.anu[nu]<=tnu.E2nu because the maximum 
		// index could be set to point to the next highest bin.
		ASSERT( rfield.anu[nu] < 1.01*tnu.E2nu || rfield.anu[nu-1]<tnu.E2nu );

		// As2nu[nu] is transition probability A per bin
		// So sum is the total transition probability
		sum += tnu.As2nu[nu];
					
		// only include this if induced processes turned on,
		// otherwise inconsistent with rate solver treatment.
		if( lgDoInduced )
		{
			// factor 0.5 accounts for fact that photons must be anti-aligned.
			double rate_up = 0.5 * tnu.As2nu[nu] *
				rfield.SummedOcc[nu] * rfield.SummedOcc[tnu.ipSym2nu[nu]-1];
			tnu.induc_up += rate_up;
			tnu.induc_dn += rate_up + tnu.As2nu[nu] *
				(rfield.SummedOcc[nu] + rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
		}
	}

	/* a sanity check on the code, see Spitzer and Greenstein (1951), eqn 4.	*/
	/* >>refer	HI	2nu	Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
	ASSERT( fabs( 1.f - (realnum)sum/tnu.AulTotal ) < 0.01f );

	return;
}

void CalcTwoPhotonEmission( two_photon& tnu, bool lgDoInduced )
{
	DEBUG_ENTRY( "CalcTwoPhotonEmission()" );
			
	/* this could fail when pops very low and Pop2Ion is zero */
	ASSERT( tnu.ipTwoPhoE>0 );

	/* two photon emission, ipTwoPhoE is 
	 * continuum array index for Lya energy */
	for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
	{
		// Pop has dimension cm^-3.  The factor of 2 is for two photons per 
		// transition. Thus two_photon_emiss has dimension photons cm-3 s-1.
		tnu.local_emis[nu] = 2.f * (realnum)(*tnu.Pop) * tnu.As2nu[nu];
	}

	// only include this if induced processes turned on,
	// otherwise inconsistent with rate solver treatment.
	if( lgDoInduced )
	{
		for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
		{
			// this is the total rate (in this energy bin)
			tnu.local_emis[nu] *= (1.f + rfield.SummedOcc[nu]) *
				(1.f+rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
		}
	}

	return;
}

/* option to print hydrogen and helium two-photon emission coefficients?	*/
void PrtTwoPhotonEmissCoef( const two_photon& tnu, const double& densityProduct )
{
	DEBUG_ENTRY( "PrtTwoPhotonEmissCoef()" );
				
	fprintf( ioQQQ, "\ny\tGammaNot(2q)\n");

	for( long yTimes20=1; yTimes20<=10; yTimes20++ )
	{
		double y = yTimes20/20.;

		fprintf( ioQQQ, "%.3e\t", (double)y );	

		long i = ipoint(y*tnu.E2nu);
		fprintf( ioQQQ, "%.3e\n", 
			8./3.*HPLANCK*(*tnu.Pop)/densityProduct*y*tnu.As2nu[i]*tnu.E2nu/rfield.widflx[i] );
	}

	return;
}