/* 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 */
/*cdSPEC returns the spectrum needed for Keith Arnaud's XSPEC */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "geometry.h"
#include "radius.h"
#include "rfield.h"
#include "opacity.h"
#include "grid.h"

/* 
 * this routine returns the spectrum needed for Keith Arnaud's XSPEC
 * X-Ray analysis code.  It should be called after cdDrive has successfully
 * computed a model.  the calling routine must ensure that the vectors
 * have enough space to store the resulting spectrum, 
 * given the bounds and energy resolution 
 */

void cdSPEC( 
	/* option - the type of spectrum to be returned
	 * 1	the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
	 *
	 * 2	the attenuated incident continuum, same units
	 * 3	the reflected continuum, same units
	 *
	 * 4	diffuse continuous emission outward direction
	 * 5	diffuse continuous emission, reflected
	 *
	 * 6	collisional+recombination lines, outward
	 * 7	collisional+recombination lines, reflected
	 * 
	 * 8	pumped lines, outward <= not implemented
	 * 9	pumped lines, reflected <= not implemented
	 *
	 *		all lines and continuum emitted by the cloud assume full coverage of 
	 *		continuum source */
	int nOption ,

	/* the number of cells + 1*/
	long int nEnergy ,

	/* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
	double ReturnedSpectrum[] )

{
	/* this pointer will bet set to one of the cloudy continuum arrays */
	realnum *cont , 
		refac;
	long int ncell , j;

	/* flag saying whether we will need to free cont at the end */
	bool lgFREE;

	DEBUG_ENTRY( "cdSPEC()" );

	ASSERT( nEnergy <= rfield.nflux );

	if( nOption == 1 )
	{
		/* this is the incident continuum, col 2 of save continuum command */
		cont = rfield.flux_total_incident[0];
		lgFREE = false;
	}
	else if( nOption == 2 )
	{
		/* the attenuated transmitted continuum, no diffuse emission,
		 * col 3 of save continuum command */
		cont = rfield.flux[0];
		lgFREE = false;
	}
	else if( nOption == 3 )
	{
		/* reflected incident continuum, col 6 of save continuum command */
		lgFREE = false;
		cont = rfield.ConRefIncid[0];
	}
	else if( nOption == 4 )
	{
		/* diffuse continuous emission outward direction  */
		cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );

		/* need to free the vector once done */
		lgFREE = true;
		refac = (realnum)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux; ++j)
		{
			cont[j] = rfield.ConEmitOut[0][j]*refac;
		}
	}
	else if( nOption == 5 )
	{
		/* reflected diffuse continuous emission */
		cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
		/* need to free the vector once done */
		lgFREE = true;
		refac = (realnum)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux; ++j)
		{
			cont[j] = rfield.ConEmitReflec[0][j]*refac;
		}
	}
	else if( nOption == 6 )
	{
		/* all outward lines,   */
		cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
		/* need to free the vector once done */
		lgFREE = true;
		/* correct back to inner radius */
		refac = (realnum)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux; ++j)
		{
			/* units of lines here are to cancel out with tricks applied to continuum cells
			 * when finally extracted below */
			cont[j] = rfield.outlin[0][j] *rfield.widflx[j]/rfield.anu[j]*refac;
		}
	}
	else if( nOption == 7 )
	{
		/* all reflected lines */
		if( geometry.lgSphere )
		{
			refac = 0.;
		}
		else
		{
			refac = 1.;
		}

		cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
		/* need to free the vector once done */
		lgFREE = true;
		for( j=0; j<rfield.nflux; ++j)
		{
			/* units of lines here are to cancel out with tricks applied to continuum cells
			 * when finally extracted below */
			cont[j] = rfield.reflin[0][j] *rfield.widflx[j]/rfield.anu[j]*refac;
		}
	}
	else
	{
		fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
		cdEXIT(EXIT_FAILURE);
	}

	/* now generate the continua */
	for( ncell = 0; ncell < nEnergy-1; ++ncell )
	{
		ReturnedSpectrum[ncell] = cont[ncell] * EN1RYD * rfield.anu2[ncell] / rfield.widflx[ncell];
	}

	/* need to free the vector once done if this flag was set */
	if( lgFREE )
	{
		free(cont);
	}
	return;
}


/* returns in units photons/cm^2/s/bin */
void cdSPEC2( 
	/* option - the type of spectrum to be returned
	 * 1	the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
	 *
	 * 2	the attenuated incident continuum, same units
	 * 3	the reflected continuum, same units
	 *
	 * 4	diffuse emission, lines + continuum, outward
	 * 5	diffuse emission, lines + continuum, reflected
	 *
	 * 6	diffuse continuous emission outward direction
	 * 7	diffuse continuous emission, reflected
	 *
	 * 8	total transmitted, lines and continuum
	 * 9	total reflected, lines and continuum
	 *
	 *10    exp(-tau) to the illuminated face
	 *
	 *		all lines and continuum emitted by the cloud assume full coverage of 
	 *		continuum source */
	int nOption ,

	/* the number of cells */
	long int nEnergy,

	/* the index of the lowest and highest energy bounds to use. */
	long ipLoEnergy,
	long ipHiEnergy,

	/* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
	realnum ReturnedSpectrum[] )

{
	realnum refac;

	DEBUG_ENTRY( "cdSPEC2()" );

	ASSERT( ipLoEnergy >= 0 );
	ASSERT( ipLoEnergy < ipHiEnergy );
	ASSERT( ipHiEnergy < rfield.nupper );
	ASSERT( nEnergy == (ipHiEnergy-ipLoEnergy+1) );
	ASSERT( nEnergy >= 2 );

	ASSERT( nOption <= NUM_OUTPUT_TYPES );

	const realnum *trans_coef_total = rfield.getCoarseTransCoef();

	for( long i = 0; i < nEnergy; i++ )
	{
		long j = ipLoEnergy + i;

		if( j >= rfield.nflux )
		{
			ReturnedSpectrum[i] = SMALLFLOAT;
			continue;
		}
		
		if( nOption == 0 )
		{
			/* the attenuated incident continuum */
			realnum flxatt = rfield.flux[0][j]*
				(realnum)radius.r1r0sq * trans_coef_total[j];

			/* the outward emitted continuum */
			realnum conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j])*
				(realnum)radius.r1r0sq*geometry.covgeo;

			/* the reflected continuum */
			realnum flxref = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
				rfield.reflin[0][j];

			ReturnedSpectrum[i] = flxatt + conem + flxref;
		}
		else if( nOption == 1 )
		{
			/* this is the incident continuum, col 2 of save continuum command */
			ReturnedSpectrum[i] = rfield.flux_total_incident[0][j];
		}
		else if( nOption == 2 )
		{
			/* the attenuated transmitted continuum, no diffuse emission,
			 * col 3 of save continuum command */
			ReturnedSpectrum[i] = rfield.flux[0][j]*
				(realnum)radius.r1r0sq * trans_coef_total[j];
		}
		else if( nOption == 3 )
		{
			/* reflected incident continuum, col 6 of save continuum command */
			ReturnedSpectrum[i] = rfield.ConRefIncid[0][j];
		}
		else if( nOption == 4 )
		{
			/* all outward diffuse emission */
			/* correct back to inner radius */
			refac = (realnum)radius.r1r0sq*geometry.covgeo;
			ReturnedSpectrum[i] = (rfield.outlin[0][j]+rfield.ConEmitOut[0][j])*refac;
		}
		else if( nOption == 5 )
		{
			/* all reflected diffuse emission */
			if( geometry.lgSphere )
			{
				refac = 0.;
			}
			else
			{
				refac = 1.;
			}

			ReturnedSpectrum[i] = (rfield.reflin[0][j]+rfield.ConEmitReflec[0][j])*refac;
		}
		else if( nOption == 6 )
		{
			/* all outward line emission */
			/* correct back to inner radius */
			refac = (realnum)radius.r1r0sq*geometry.covgeo;
			ReturnedSpectrum[i] = rfield.outlin[0][j]*refac;
		}
		else if( nOption == 7 )
		{
			/* all reflected line emission */
			if( geometry.lgSphere )
			{
				refac = 0.;
			}
			else
			{
				refac = 1.;
			}

			ReturnedSpectrum[i] = rfield.reflin[0][j]*refac;
		}
		else if( nOption == 8 )
		{
			/* total transmitted continuum */
			/* correct back to inner radius */
			refac = (realnum)radius.r1r0sq*geometry.covgeo;
			ReturnedSpectrum[i] = (rfield.ConEmitOut[0][j]+ rfield.outlin[0][j])*refac
				+ rfield.flux[0][j]*(realnum)radius.r1r0sq*trans_coef_total[j];
		}
		else if( nOption == 9 )
		{
			/* total reflected continuum */
			ReturnedSpectrum[i] = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
				rfield.reflin[0][j];
		}
		else if( nOption == 10 )
		{
			/* this is exp(-tau) */
			/* This is the TOTAL attenuation in both the continuum and lines.  
			 * Jon Miller discovered that the line attenuation was missing in 07.02 */
			ReturnedSpectrum[i] = opac.ExpmTau[j]*trans_coef_total[j];
		}
		else
		{
			fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
			cdEXIT(EXIT_FAILURE);
		}

		ASSERT( ReturnedSpectrum[i] >=0.f );
	}

	return;
}