/* 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 */
/*RT_continuum attenuation of diffuse and beamed continua */
/*pnegopc save negative opacities on io unit, iff 'set negopc' command was given */
#include "cddefines.h"
#include "rfield.h"
#include "opacity.h"
#include "dense.h"
#include "colden.h"
#include "opacity.h"
#include "geometry.h"
#include "trace.h"
#include "radius.h"
#include "iso.h"
#include "hextra.h"

/*cmshft - so Compton scattering shift of spectrum 
 * this code is a placeholder */
STATIC void cmshft(void)
{
	DEBUG_ENTRY( "cmshft()" );

	/* first check whether Compton scattering is in as heat/cool */
	if( !rfield.lgComptonOn )
	{ 
		return;
	}

	if( rfield.lgComptonOn )
	{ 
		return;
	}

	/* do reshuffle */
	for( long i=0; i < rfield.nflux; i++ )
	{
		continue;
	}
	return;
}


#if !defined(NDEBUG)
/*pnegopc save negative opacities on io unit, iff 'set negopc' command was given */
STATIC void pnegopc(void)
{
	FILE *ioFile;

	DEBUG_ENTRY( "pnegopc()" );

	if( opac.lgNegOpacIO )
	{
		/* option to save negative opacities */
		ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY );
		for( long i=0; i < rfield.nflux; i++ )
		{
			fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i], 
				opac.opacity_abs[i] );
		}
		fclose( ioFile);
	}
	return;
}
#endif

/*RT_continuum attenuation of diffuse and beamed continua */
void RT_continuum(void)
{

	DEBUG_ENTRY( "RT_continuum()" );

	if( trace.lgTrace && trace.lgConBug )
	{
		fprintf( ioQQQ, " Energy, flux, OTS:\n" );
		for( long i=0; i < rfield.nflux; i++ )
		{
			fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 
				rfield.flux[0][i] + rfield.outlin[0][i] + rfield.ConInterOut[i], 
				rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
		}
		fprintf( ioQQQ, "\n" );
	}

	/* begin sanity check in debug mode */
#	if !defined(NDEBUG)
	bool lgFlxNeg = false;
	for( long i=0; i < rfield.nflux; i++ )
	{
		if( rfield.flux[0][i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 
				rfield.flux[0][i], rfield.anu[i], i );
			lgFlxNeg = true;
		}
		if( rfield.otscon[i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 
				rfield.otscon[i], rfield.anu[i], i );
			lgFlxNeg = true;
		}
		if( opac.tmn[i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
			fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
				opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
			lgFlxNeg = true;
		}
		if( rfield.otslin[i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
				rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i]  );
			lgFlxNeg = true;
		}
		if( rfield.outlin[0][i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
				rfield.outlin[0][i], rfield.anu[i], i, rfield.chLineLabel[i]  );
			lgFlxNeg = true;
		}
		if( rfield.ConInterOut[i] < 0. )
		{
			fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
				rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i]  );
			lgFlxNeg = true;
		}
		if( opac.opacity_abs[i] < 0. )
		{
			opac.lgOpacNeg = true;
			/* this sub will save negative opacities on io unit,
			 * iff 'set negopc' command was given */
			pnegopc();
		}
	}
	if( lgFlxNeg )
	{
		fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n", 
			nzone );
		ShowMe();
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check*/
#	endif

	/* rfield.lgOpacityFine flag set false with no fine opacities command 
	 * tests show that always evaluating this changes fast run of
	 * pn_paris from 26.7 sec to 35.1 sec 
	 * but must always update fine opacities since used for transmission */
	if( rfield.lgOpacityFine )
	{
		/* increment the fine opacity array */
		for( long i=0; i<rfield.nfine; ++i )
		{
			realnum tauzone = rfield.fine_opac_zone[i]*(realnum)radius.drad_x_fillfac;
			rfield.fine_opt_depth[i] += tauzone;
		}
		rfield.trans_coef_total_stale = true;
	}

	// ratio of inner to outer radii, at this point
	// radius is the outer radius of this zone 
	double DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
		radius.Radius);

	rfield.EnergyIncidCont = 0.;
	rfield.EnergyDiffCont = 0.;

	// this loop should not be to <= nflux since highest cell is for
	// continuum unit integration
	// scattering opacities are included in energy exchange here in the
	// sphere case, since photons diffuse out of the closed sphere.
	// scattering opacities are not included as extinction sources in the
	// sphere case
	for( long i=0; i < rfield.nflux; i++ )
	{
		double dTau_abs = opac.opacity_abs[i]*radius.drad_x_fillfac;
		double dTau_sct = opac.opacity_sct[i]*radius.drad_x_fillfac;

		// sum total continuous optical depths
		opac.TauAbsGeo[0][i] += (realnum)(dTau_abs);
		opac.TauScatGeo[0][i] += (realnum)(dTau_sct);

		// following only optical depth to illuminated face
		opac.TauAbsFace[i] += (realnum)(dTau_abs);
		opac.TauScatFace[i] += (realnum)(dTau_sct);

		// these are total in inward direction, large if spherical
		opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];

		// attenuation of flux by optical depths IN THIS ZONE 
		// DirectionalCosin is 1/COS(theta), is usually 1, reset with illuminate command,
		// option for illumination of slab at an angle 
		opac.ExpZone[i] = sexp(dTau_abs*geometry.DirectionalCosin);

		// e(-tau) in inward direction, up to illuminated face
		opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];

		// e2(tau) in inward direction, up to illuminated face 
		opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
		ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );

		// on second and later iterations define outward E2
		if( iteration > 1 )
		{
			// e2 from current position to outer edge of shell 
			realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
			opac.E2TauAbsOut[i] = (realnum)e2( tau );
			ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
		}

		// DilutionHere is square of ratio of inner to outer radius
		double AttenuationDilutionFactor = opac.ExpZone[i]*DilutionHere;
		ASSERT( AttenuationDilutionFactor <= 1.0 );

		// continuum has three parts 
		rfield.flux_beam_const[i] *= (realnum)AttenuationDilutionFactor;
		rfield.flux_beam_time[i] *= (realnum)AttenuationDilutionFactor;
		rfield.flux_isotropic[i] *= (realnum)AttenuationDilutionFactor;
		rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
			rfield.flux_isotropic[i];

		// update SummedCon here since flux changed
		rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];

		// outward lines and diffuse continua 
		rfield.outlin[0][i] *= (realnum)AttenuationDilutionFactor;
		rfield.outlin_noplot[i] *= (realnum)AttenuationDilutionFactor;

		// interactive outward diffuse continuum 
		TestCode();// move this over from rt_diffuse - this preserves
				   // the original order and is incorrect 
		rfield.ConInterOut[i] += rfield.DiffuseEscape[i]*(realnum)radius.dVolOutwrd;
		rfield.ConInterOut[i] *= (realnum)AttenuationDilutionFactor;

		// this is not the interacting continuum 
		rfield.ConEmitOut[0][i] *= (realnum)AttenuationDilutionFactor;
		rfield.ConEmitOut[0][i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i]/**/;

		// set occupation numbers, first attenuated incident continuum 
		rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];

		// local diffuse continua 
		rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];

		// outward diffuse continuum 
		rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[0][i]*rfield.convoc[i];

		// integrated energy flux, ergs s^-1 cm^-2 
		rfield.EnergyIncidCont += rfield.flux[0][i]*rfield.anu[i];
		rfield.EnergyDiffCont += (rfield.outlin[0][i] + rfield.outlin_noplot[i] +
			rfield.ConInterOut[i])* rfield.anu[i];
	}

	// convert Ryd to erg 
	rfield.EnergyIncidCont *= (realnum)EN1RYD;
	rfield.EnergyDiffCont *= (realnum)EN1RYD;

	// sanity check, compare total Lyman continuum optical depth 
	// with amount of extinction there
	// this is amount continuum attenuated to illuminated face, 
	// but only do test if flux positive, not counting scattering opacity,
	// and correction for spherical dilution not important 
	if( rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]>SMALLFLOAT &&
		(rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
		SDIV(rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) ) > SMALLFLOAT && 
		!opac.lgScatON &&
		radius.depth/radius.Radius < 1e-4 )
	{
		// ratio of current to incident continuum, converted to optical depth
		/* >>chng 99 apr 23, this crashed on alpha due to underflow to zero in argy
		* to log.  defended two ways - above checks that ratio of fluxes is large enough, 
		* and here convert to double.
		* error found by Peter van Hoof */
		double tau_effec = -log((double)rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
			(double)opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
			(double)rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]);

		// this is computed absorption optical depth to illuminated face
		double tau_true = opac.TauAbsFace[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]*geometry.DirectionalCosin;

		// first test is relative error, second is to absolute error and comes
		// in for very small tau, where differences are in the round-off 
		if( fabs( tau_effec - tau_true ) / MAX2(tau_effec , tau_true) > 0.01 && 
			// for very small inner optical depths, the tmn correction is major,
			// and this test is not precise
			fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) )
		{
			// in print below must add extra HI col den since this will later be
			// incremented in RT_tau_inc 
			fprintf( ioQQQ,
				" PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
				nzone, 
				tau_effec , 
				tau_true ,
				6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
			TotalInsanity();
		}
	}

	// do scattering opacity, not included when sphere is set
	if( opac.lgScatON )
	{
		for( long i=0; i < rfield.nflux; i++ )
		{
			// Lightman and White equation 11 in small epsilon limit,
			// >>refer	continuum	RT	Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */
			double AttenuationScatteringFactor = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
			ASSERT( AttenuationScatteringFactor <= 1.0 );
			rfield.flux_beam_const[i] *= (realnum)AttenuationScatteringFactor;
			rfield.flux_beam_time[i] *= (realnum)AttenuationScatteringFactor;
			rfield.flux_isotropic[i] *= (realnum)AttenuationScatteringFactor;
			rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
				rfield.flux_isotropic[i];

			rfield.ConInterOut[i] *= (realnum)AttenuationScatteringFactor;
			rfield.ConEmitOut[0][i] *= (realnum)AttenuationScatteringFactor;
			rfield.outlin[0][i] *= (realnum)AttenuationScatteringFactor;
			rfield.outlin_noplot[i] *= (realnum)AttenuationScatteringFactor;
		}
	}

	// this dilution is needed to conserve volume in spherical models.  tests such
	// as parispn.in will fault if this is removed
	realnum Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );

	// this is a unit of energy that will be passed through the code as a test
	// that all integrations are carried out.  A similar test is set in lineset1
	// and verified in PrtFinal.  The opacity at this cell is zero so only
	// geometrical dilution will affect the integral
	// Radius is currently outer edge of zone, so radius-drad/2 is radius
	// of center of zone 
	rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution;
	rfield.DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;
	// must do unit integration somewhere 
	rfield.ConInterOut[rfield.nflux] += 
		rfield.DiffuseEscape[rfield.nflux]*(realnum)radius.dVolOutwrd;

	// opacity should be zero at this energy so J not changed elsewhere
	ASSERT( opac.opacity_abs[rfield.nflux] == 0. );

	// placeholder code for Compton scattering 
	cmshft();

	// attenuate neutrons if they are present 
	if( hextra.lgNeutrnHeatOn )
	{
		// correct for optical depth effects 
		hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*
			dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
		// correct for spherical effects 
		hextra.totneu *= (realnum)DilutionHere;
	}

	// following radiation factors are extinguished by 1/r**2ilution, electron
	// scattering by free and bound electrons

	// do all emergent spectrum from illuminated face if model is NOT spherical
	if( !geometry.lgSphere )
	{
		double Reflec_Diffuse_Cont;

		// emission starting at the the plasma frequency
		for( long i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
		{
			if( opac.TauAbsGeo[0][i] < 30. )
			{
				// ConEmitLocal is diffuse emission per unit vol, fill factor
				// the 1/2 comes from isotropic emission 
				Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
					radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;

				// ConEmitReflec - reflected diffuse continuum
				rfield.ConEmitReflec[0][i] += (realnum)(Reflec_Diffuse_Cont);

				// the reflected part of the incident continuum
				rfield.ConRefIncid[0][i] += (realnum)(rfield.flux[0][i]*opac.opacity_sct[i]*
					radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);

				// reflected line emission 
				rfield.reflin[0][i] += (realnum)(rfield.outlin[0][i]*opac.opacity_sct[i]*
					radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
			}
		}
	}
	return;
}