/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*cdSPEC returns the spectrum needed for Keith Arnaud's XSPEC */
/*Spec_cont called by cdSPEC to generate actual spectrum */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "geometry.h"
#include "radius.h"
#include "rfield.h"
#include "opacity.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 
 */

/* array index within energy array in Cloudy grids -  this has file scope */
static long int iplo , iphi;

/* energies of lower and upper bound of XSPEC cell we are considering */
static double Elo , Ehi;

/* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
static double Spec_cont( float cont[] );

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 energy of the lower edge of each cell 
	 * (in Ryd to be consistent with all the rest of Cloudy) */
	double EnergyLow[] , 

	/* 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 */
	float *cont , 
		refac;
	long int ncell , j;

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

#	ifdef DEBUG_FUN
	fputs( "<+>cdSPEC()\n", debug_fp );
#	endif

	if( nOption == 1 )
	{
		/* this is the incident continuum, col 2 of punch continuum command */
		cont = rfield.FluxSave;
		lgFREE = FALSE;
	}
	else if( nOption == 2 )
	{
		/* the attenuated transmitted continuum, no diffuse emission,
		 * col 3 of punch continuum command */
		cont = rfield.flux;
		lgFREE = FALSE;
	}
	else if( nOption == 3 )
	{
		/* reflected incident continuum, col 6 of punch continuum command */
		lgFREE = FALSE;
		cont = rfield.ConRefIncid ;
	}
	else if( nOption == 4 )
	{
		/* diffuse continuous emission outward direction  */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 4 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		refac = (float)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux ; ++j)
		{
			cont[j] = rfield.ConEmitOut[j]*refac ;
		}
	}
	else if( nOption == 5 )
	{
		/* reflected diffuse continuous emission */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 5 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		refac = (float)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux ; ++j)
		{
			cont[j] = rfield.ConEmitReflec[j]*refac  ;
		}
	}
	else if( nOption == 6 )
	{
		/* all outward lines,   */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 6 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		/* correct back to inner radius */
		refac = (float)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[j] *rfield.widflx[j]/rfield.anu[j]*refac ;
		}
	}
	else if( nOption == 7 )
	{
		/* all reflected lines */
		if( geometry.lgSphere )
		{
			refac = 0.;
		}
		else
		{
			refac = 1.;
		}

		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 7 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* 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[j] *rfield.widflx[j]/rfield.anu[j]*refac ;
		}
	}
	else
	{
		fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
		puts( "[Stop in cdSPEC]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this is array index used in Spec_cont */
	iplo = 0;
	iphi = 0;
	/* now generate the continua */
	for( ncell = 0; ncell < nEnergy-1; ++ncell )
	{
		/* Spec_cont will find the spectrum between these two energies */
		Elo = EnergyLow[ncell];
		Ehi = EnergyLow[ncell+1];
		ReturnedSpectrum[ncell] = Spec_cont( cont );
	}

	/* need to free the vector once done if this flag was set */
	if( lgFREE )
	{
		free(cont) ;
	}
#	ifdef DEBUG_FUN
	fputs( " <->cdSPEC()\n", debug_fp );
#	endif
	return;
}

/* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
static double Spec_cont( float cont[] )
{
	double sum;
	long int i;

	/* iplo and iphi are set to zero before this is called to generate a full spectrum,
	 * since energies are in increasing order, after the first call it should only
	 * be necessary to increase ip above where it now is */
	/* return zero if continuum does not go this high */
	if( iplo >= rfield.nflux )
		return 0.;

	/* now skip out to where continuum starts, often we will already be there
	 * if this is a successive call */
	/* iplo should be index where Elo is greater than anu-widflx/1. */
	while( (iplo < rfield.nflux)  &&  
		!(Elo <= (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) &&
		 (Elo >= (rfield.AnuOrg[iplo]-rfield.widflx[iplo]/2.)  )) )
		++iplo;

	iphi = iplo;
	/* iphi should be index where Ehi is greater than anu-widflx/1. */
	while( (iphi < rfield.nflux)  &&  
		!(Ehi <= (rfield.AnuOrg[iphi+1]-rfield.widflx[iphi+1]/2.) &&
		 (Ehi >= (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.)  )) )
		++iphi;

	sum = 0.;
	if( iphi-iplo>1 )
	{
		/* this branch when more than two cells are involved - 
		 * begin by summing intermediate cells */
		for( i=iplo+1; i<iphi; ++i )
		{
			sum += cont[i] * rfield.anu2[i];
		}
	}
	ASSERT( sum>=0.);

	/* at this point sum will often be zero, and we need to add on flux from
	 * the iplo and iphi cells - there are two possibilities - iplo and iphi can
	 * be the same cell, or different cells */
	if( iplo == iphi )
	{
		/* we want only a piece of the ip cell */
		sum = cont[iplo]/rfield.widflx[iplo]*(Ehi-Elo) * rfield.anu2[iplo];
		ASSERT( sum>=0.);
	}
	else
	{
		/* we want to add on a fraction of the low and high energy cells */
		double piece;

		/* the low energy cell */
		piece = cont[iplo]/rfield.widflx[iplo]*( (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) - Elo ) * 
			rfield.anu2[iplo];
		ASSERT( piece >= 0.);
		sum += piece;

		/* now the high energy cell */
		piece = cont[iphi]/rfield.widflx[iphi]*(Ehi - (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) ) * 
			rfield.anu2[iphi];
		ASSERT( piece >= 0.);
		sum += piece;
		ASSERT( sum>=0.);
	}

	/* now divide by total energy width */
	sum *= EN1RYD / (Ehi - Elo );

	ASSERT( sum >=0. );
	return sum;

#	if 0
	/* at this point Elo is less than the upper edge energy of the ip cell */
	/* >>chng 01 jul 23, change logic, use this branch if desired width smaller than
	 * intrinsic cloudy width */
	/*if( rfield.anu[ip]+rfield.widflx[ip]/2. > Ehi )*/
	if( (Ehi-Elo) < rfield.widflx[ip] )
	{
		/* this is where we want the flux evaluated */
		double Emiddle = (Elo + Ehi)/2.;
		double f1 , f2 ,finterp;
		ASSERT( (Elo >= (rfield.anu[ip]+rfield.widflx[ip]/2.
		/* this is case where requested cell is within a single Cloudy cell -
		 * do linear interpolation */
		/* >>chng 01 jul 23, change to linear interpolation */
		/*return( cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD );*/
		f1 = cont[ip-1] /rfield.widflx[ip-1] *rfield.anu2[ip-1]*EN1RYD;
		f2 = cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD;

		finterp = f1 + (f2-f1) /rfield.widflx[ip] *
			fabs(Emiddle-(rfield.anu[ip]-rfield.widflx[ip]/2.));
		return( finterp );
	}
	else
	{
		/* this is case where we will add on more than one Cloudy cell */
		sum = cont[ip] * rfield.anu2[ip]* EN1RYD ;
		/* energy of upper edge of Cloudy cell */
		EcHi = ( rfield.anu[ip] + rfield.widflx[ip]/2. + rfield.anu[ip+1] - rfield.widflx[ip+1]/2.)/2.;
		/* mutiply this by energy width of this first cell that counts to final integral */
		sum *= (EcHi - Elo )/rfield.widflx[ip];
	}

	/* increment index since we will now work on next cell */
	++ip;

	/* now add up total cells, this loop for cloudy cells totally within desired bins,
	 * we will divide by energy width later, so factor of widflx is missing here */
	while( Ehi > rfield.anu[ip]+rfield.widflx[ip]/2. )
	{
		sum += cont[ip] * rfield.anu2[ip]*EN1RYD ;
		++ip;
	}

	/* energy of upper edge of Cloudy cell */
	EcLo = ( rfield.anu[ip] - rfield.widflx[ip]/2. + rfield.anu[ip-1] + rfield.widflx[ip-1]/2.)/2.;

	/* now finish up with last cell, whose upper bound is higher than desired bin */
	sum += cont[ip] * rfield.anu2[ip] * EN1RYD *
		(Ehi - EcLo )/rfield.widflx[ip];

	/* now divide by total energy width */
	sum /= (Ehi - Elo );

	return(sum);
#	endif
}

/* 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 continuous emission outward direction
	 * 5	diffuse continuous emission, reflected
	 *
	 * 6	collisional+recombination lines, outward
	 * 7	collisional+recombination lines, reflected
	 *******************************************
	 * NB - for now 8 is this:
	 * 8    exp(-tau) to the illuminated face 
	 *******************************************
	 * 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 */
	float ReturnedSpectrum[] )

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

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

#	ifdef DEBUG_FUN
	fputs( "<+>cdSPEC2()\n", debug_fp );
#	endif

	if( nOption == 1 )
	{
		/* this is the incident continuum, col 2 of punch continuum command */
		cont = rfield.FluxSave;
		lgFREE = FALSE;
	}
	else if( nOption == 2 )
	{
		/* the attenuated transmitted continuum, no diffuse emission,
		 * col 3 of punch continuum command */
		cont = rfield.flux;
		lgFREE = FALSE;
	}
	else if( nOption == 3 )
	{
		/* reflected incident continuum, col 6 of punch continuum command */
		lgFREE = FALSE;
		cont = rfield.ConRefIncid ;
	}
	else if( nOption == 4 )
	{
		/* diffuse continuous emission outward direction  */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 4 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		refac = (float)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux ; ++j)
		{
			cont[j] = rfield.ConEmitOut[j]*refac ;
		}
	}
	else if( nOption == 5 )
	{
		/* reflected diffuse continuous emission */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 5 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		refac = (float)radius.r1r0sq*geometry.covgeo;
		for( j=0; j<rfield.nflux ; ++j)
		{
			cont[j] = rfield.ConEmitReflec[j]*refac  ;
		}
	}
	else if( nOption == 6 )
	{
		/* all outward lines,   */
		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 6 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* need to free the vector once done */
		lgFREE = TRUE;
		/* correct back to inner radius */
		refac = (float)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 
			 *>>chng 05 sep 29, RP add ConEmitOut to following, is this correct?
			 * is this just lines of also continuum? */
			cont[j] = (rfield.ConEmitOut[j]+rfield.outlin[j]*rfield.widflx[j]/rfield.anu[j])*refac ;
		}
	}
	else if( nOption == 7 )
	{
		/* all reflected lines */
		if( geometry.lgSphere )
		{
			refac = 0.;
		}
		else
		{
			refac = 1.;
		}

		if( (cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper ) ) == NULL)
		{
			fprintf(ioQQQ," cdSPEC could not MALLOC nOption 7 \n");
			cdEXIT(EXIT_FAILURE);
		}
		/* 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[j] *rfield.widflx[j]/rfield.anu[j]*refac ;
		}
	}
	else if( nOption == 8 )
	{
		/* this is exp(-tau) */
		for( j=0; j<nEnergy ; ++j)
		{
			ReturnedSpectrum[j] = opac.ExpmTau[j];
		}
		return;
	}
	else
	{
		fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
		puts( "[Stop in cdSPEC]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this is array index used in Spec_cont */
	iplo = 0;
	iphi = 0;
	/* now generate the continua */
	for( ncell = 0; ncell < nEnergy; ++ncell )
	{
        if( ncell >= rfield.nflux )
		{
			ReturnedSpectrum[ncell] = 0.f;
		}
		else
		{
			ReturnedSpectrum[ncell] = cont[ncell] * rfield.anu2[ncell] * ((float)EN1RYD) / rfield.widflx[ncell];
		}
		ASSERT( ReturnedSpectrum[ncell] >=0.f );
	}

	/* need to free the vector once done if this flag was set */
	if( lgFREE )
	{
		free(cont) ;
	}
#	ifdef DEBUG_FUN
	fputs( " <->cdSPEC2()\n", debug_fp );
#	endif
	return;
}
