/* 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 */
/*lines_continuum put energetics, H, and He lines into line intensity stack */
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "iso.h"
#include "geometry.h"
#include "dense.h"
#include "prt.h"
#include "coolheavy.h"
#include "phycon.h"
#include "rfield.h"
#include "predcont.h"
#include "lines_service.h"
#include "radius.h"
#include "continuum.h"
#include "lines.h"

void lines_continuum(void)
{

	double f1, 
		f2 , 
		bac , 
		flow;
	long i,nBand;

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

	/* code has all local emissivities zeroed out with cryptic comment about being
	 * situation dependent.  Why?  this is option to turn back on */
#	define KILL_CONT	FALSE

	/* some continua */
	i = StuffComment( "continua" );
	linadd( 0., (float)i , "####", 'i');

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/***********************************************************************
	 * stuff in Bac ratio - continuum above the BJ 
	 * this is trick, zeroing out saved continuum integrated so far,
	 * and adding the current version, so that the line array gives the 
	 * value in the final continuum 
	 *
	 * reflec continuum is different from others since relative to inner
	 * radius, others for for this radius 
	 *************************************************************************/

	/*TODO	2	this block of lines should have nInu, InwT, InwC like main vector of continuum points */
	/***************************************************************************
	 * "Bac " , 3646,  this is residual continuum at peak of Balmer Jump
	 * flux below - flux above                                   
	 ***************************************************************************/
	/* >>chng 00 dec 02, remove opac.tmn */
	/* >>chng 00 dec 19, remove / radius.GeoDil */
	/* extrapolated continuum above head */
	/* >>chng 01 jul 13, from ConInterOut to ConEmitOut */
	f1 = (rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1] + 
		rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/radius.r1r0sq )/
	  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;

	/* extrapolated continuum below head */
	/* >>chng 00 dec 19, remove / radius.GeoDil */
	f2 = (rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]+ 
		rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/radius.r1r0sq )/
	  rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;

	/* convert to nuFnu units */
	f1 = f1*0.250*0.250*EN1RYD*radius.r1r0sq;
	f2 = f2*0.250*0.250*EN1RYD*radius.r1r0sq;
	bac = (f1 - f2);

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;
	/* residual flux at head of Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"Bac ",'i');
	/* >>chng 03 feb 06, set to zero */
	/* emslin saves the per unit vol emissivity of a line, which is normally 
	 * what goes into linadd.  We zero this unit emissivity which was set
	 * FOR THE PREVIOUS LINE since it is so situation dependent */
	if( KILL_CONT && LineSave.ipass > 0 )
		LineSv[LineSave.nsum-1].emslin = 0.;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;
	/* total flux above head of Balmer continuum, nuFnu */
	linadd(f1/radius.dVeff,3645,"nFnu",'i');
	/* >>chng 03 feb 06, set to zero */
	/* emslin saves the per unit vol emissivity of a line, which is normally 
	 * what goes into linadd.  We zero this unit emissivity which was set
	 * FOR THE PREVIOUS LINE since it is so situation dependent */
	if( KILL_CONT && LineSave.ipass > 0 )
		LineSv[LineSave.nsum-1].emslin = 0.;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;
	/* total flux above head of Balmer continuum, nuFnu */
	linadd(f2/radius.dVeff,3647,"nFnu",'i');
	/* >>chng 03 feb 06, set to zero */
	/* emslin saves the per unit vol emissivity of a line, which is normally 
	 * what goes into linadd.  We zero this unit emissivity which was set
	 * FOR THE PREVIOUS LINE since it is so situation dependent */
	if( KILL_CONT && LineSave.ipass > 0 )
		LineSv[LineSave.nsum-1].emslin = 0.;

	/******************************************************************************
	 * "cout" , 3646,  this is outward residual continuum at peak of Balmer Jump  *
	 * equal to total in spherical geometry, half in opt thin open geo            *
	 ******************************************************************************/
	/* >>chng 00 dec 02, remove opac.tmn */
	/* >>chng 00 dec 19, remove / radius.GeoDil */
	f1 = rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/
		rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;

	/* >>chng 00 dec 19, remove / radius.GeoDil */
	f2 = rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/
		rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;

	/* net balmer jump */
	bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"cout",'i');
	/* >>chng 03 feb 06, set to zero */
	/* emslin saves the per unit vol emissivity of a line, which is normally 
	 * what goes into linadd.  We zero this unit emissivity which was set
	 * FOR THE PREVIOUS LINE since it is so situation dependent */
	if( KILL_CONT && LineSave.ipass > 0 )
		LineSv[LineSave.nsum-1].emslin = 0.;

	/*********************************************************************
	 * "cref" , 3646,  this is reflected continuum at peak of Balmer Jump*
	 * equal to zero in spherical geometry, half of total in op thin opn *
	 *********************************************************************/
	/* >>chng 00 dec 02, remove opac.tmn */
	/* >>chng 00 dec 19, remove / radius.GeoDil */
	f1 = rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/
		rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;

	f2 = rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/
		rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;

	/* net balmer jump */
	bac = (f1 - f2)*0.250*0.250*EN1RYD;

	/* memory not allocated until ipass >= 0 */
	if( LineSave.ipass > 0 )
		LineSv[LineSave.nsum].sumlin = 0.;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac)/radius.dVeff,3646,"cref",'i');
	/* >>chng 03 feb 06, set to zero */
	/* emslin saves the per unit vol emissivity of a line, which is normally 
	 * what goes into linadd.  We zero this unit emissivity which was set
	 * FOR THE PREVIOUS LINE since it is so situation dependent */
	if( KILL_CONT && LineSave.ipass > 0 )
		LineSv[LineSave.nsum-1].emslin = 0.;

	/*********************************************************************
	 * "thin" , 3646, tot optically thin continuum at peak of Balmer Jump*/
	f1 = rfield.ConEmitLocal[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1];
	f2 = rfield.ConEmitLocal[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2];
	bac = (f1 - f2)*0.250*0.250*EN1RYD;

	/* residual flux in Balmer continuum, nuFnu */
	linadd(MAX2(0.,bac),3646,"thin",'i');

	/* incident continuum nu*f_nu at H-beta, at illuminated face of cloud  */
	continuum.cn4861 /= (float)radius.dVeff;
	linadd(continuum.cn4861,4860,"Inci",'i');

	/* incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud */
	continuum.cn1216 /= (float)radius.dVeff;
	linadd(continuum.cn1216,1215,"Inci",'i');

	continuum.cn1216 *= (float)radius.dVeff;
	continuum.cn4861 *= (float)radius.dVeff;

	if( LineSave.ipass > 0 )
	{
		continuum.cn4861 = 0.;
		continuum.cn1216 = 0.;
	}

	/* integrated Balmer continuum emission */
	flow = (iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2p][ipRecRad] + iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2s][ipRecRad])*
	  iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2p][ipRecEsc]*dense.eden*dense.xIonDense[ipHYDROGEN][1]*
	  5.45e-12;
	linadd(flow,0,"Ba C",'i');

	/* >>chng 96 june 13, had been hrec(3,1,ipRecEsc)
	 * Paschen continuum emission */
	flow = iso.RadRecomb[ipH_LIKE][ipHYDROGEN][3][ipRecRad]*iso.RadRecomb[ipH_LIKE][ipHYDROGEN][3][ipRecEsc]*
	  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
	linadd(flow,0,"PA C",'i');

	/* >>chng 05 oct 16, add this output option, these are a
	 * series of continuum bands defined in the file bands_continuum.dat
	 * this makes it possible to enter any integrated total emission
	 * into the emission-line stack */
	for( nBand=0; nBand<continuum.nContBand; ++nBand )
	{
		double total;
		/* find total emission over band */
		total = 0.;
		for( i=continuum.ipContBandLow[nBand]; i<continuum.ipContBandHi[nBand]; ++i )
		{
			total += rfield.anu[i-1] *
				(
				/* the attenuated incident continuum */
				rfield.flux[i-1] +

				/* the outward emitted continuum radiation field */
				rfield.ConEmitOut[i-1] + 

				/* outward emitted lines */
				rfield.outlin[i-1]*rfield.anu[i-1]*geometry.covgeo
				);
		}

		/* convert to physical units */
		total *= EN1RYD*radius.r1r0sq/radius.dVeff;
		/* memory not allocated until ipass >= 0 */
		if( LineSave.ipass > 0 )
			LineSv[LineSave.nsum].sumlin = 0.;

		linadd( total , 
			continuum.ContBandWavelength[nBand] , 
			continuum.chContBandLabels[nBand],
			'i' );
	}

	/* net free-free cooling, ALL species, free-free heating subtracted,
	 * so nearly cancels with cooling in lte */
	linadd(MAX2(0.,CoolHeavy.brems_cool_net),0,"HFFc",'c');

	/* net free-free heating, nearly cancels with cooling in lte */
	linadd(MAX2(0.,-CoolHeavy.brems_cool_net),0,"HFFh",'h');

	/* H brems (free-free) cooling */
	linadd(CoolHeavy.brems_cool_h,0,"H FF",'i');

	/* total free-free heating */
	linadd(CoolHeavy.brems_heat_total,0,"FF H",'i');

	/* He brems emission */
	linadd(CoolHeavy.brems_cool_he,0,"HeFF",'i');

	/* He recombination cooling */
	linadd(CoolHeavy.herec,0,"HeFB",'c');

	/* heavy element recombination cooling */
	linadd(CoolHeavy.heavfb,0,"MeFB",'c');

	/* meavy elements (metals) brems cooling, heat not subtracted */
	linadd(CoolHeavy.brems_cool_metals,0,"MeFF",'i');

	/* total brems emission - total cooling but not minus heating */
	linadd(CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he+CoolHeavy.brems_cool_metals,0,"ToFF",'i');

	/* part of H brems, in x-ray beyond 0.5KeV */
	linadd((CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he)*sexp(5.8e6/phycon.te),0,"FF X"
	  ,'i');

	/* electron - electron brems */
	linadd(CoolHeavy.eebrm,0,"eeff",'c');

	/* predict emitted continuum at series of continuum points */
	/* variables are located in predcont.h, 
	 * NPREDCONT number of continuum points,
	 * next two arrays are defined in zerologic.c
	 * EnrPredCont - energy (Ryd) where we want to predict the continuum,
	 *               these are set in zerologic.c,
	 * ipPreCont -   pointers to energy within continuum array, 
	 *               these are set in ContCreatePointers
	 *
	 * ths array of continuum energies where this will be done is
	 * stored in the array EnrPredCont which lives in predcont.h
	 * and is initialized in zerologic.c
	 *
	 * the entry nFnu will only be printed if the command
	 * print diffuse continuum
	 * is entered - 
	 *
	 * this code should be kept parallel with that in dopunch, where
	 * punch continuum is produced, since two must agree */
	for( i=0; i < NPREDCONT; i++ )
	{
		double SourceTransmitted , Cont_nInu ;
		double SourceReflected, 
		  DiffuseOutward,
		  DiffuseInward;
		double renorm;

		/* put wavelength in Angstroms into dummy structure, so that we can use iWavLen
		 * to get a proper wavelength with units, continuum energies are stored in EnrPredCont */
		TauDummy.WLAng = (float)(RYDLAM/EnrPredCont[i]);
		/*lamba = iWavLen(&TauDummy , &chUnits , &chShift );*/

		/* >>chng 00 dec 02, there were three occurances of /opac.tmn which had the
		 * effect of raising the summed continuum by the local opacity correction factor.
		 * in the case of the Lyman continuum this raised the reported value by orders
		 * of magnitude.  There have been commented out in the following for now. */
		/* reflected toal continuum (dif+incident emitted inward direction) */

		/* >>chng 00 dec 08, implement the "set nFnu [SOURCE_REFLECTED] ... command, PvH */
		/* >>chng 00 dec 19, remove / radius.GeoDil */
		renorm = rfield.anu2[ipPredCont[i]]*EN1RYD/rfield.widflx[ipPredCont[i]]/*/radius.GeoDil*/;
		  
		/* this is the reflected diffuse continuum */
		if( prt.lgDiffuseInward )
		{
			DiffuseInward = rfield.ConEmitReflec[ipPredCont[i]]*renorm;
			/* /opac.tmn[ipPredCont[i]]*/
		}
		else
		{
			DiffuseInward = 0.;
		}

		/* the outward diffuse continuum */
		if( prt.lgDiffuseOutward )
		{
			DiffuseOutward = rfield.ConEmitOut[ipPredCont[i]]*renorm*radius.r1r0sq ;
		}
		else
		{
			DiffuseOutward = 0. ;
		}

		/* reflected part of INCIDENT continuum (only incident, not diffuse, which was above) */
		if( prt.lgSourceReflected )
		{
			SourceReflected =  rfield.ConRefIncid[ipPredCont[i]]*renorm;
			/* /opac.tmn[ipPredCont[i]]*/
		}
		else
		{
			SourceReflected =  0.;
		}

		/* the attenuated incident continuum */
		if( prt.lgSourceTransmitted )
		{
			SourceTransmitted = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq ;
		}
		else
		{
			SourceTransmitted = 0. ;
		}

		/* memory has not been allocated until ipass >= 0, so must not access this element,
		 * this element will be used to save the following quantity */
		if( LineSave.ipass > 0 )
			LineSv[LineSave.nsum].sumlin = 0.;

		/* >>chng 00 dec 02, did not include SourceReflected */
		/* total continuum produced by cloud at selected energy points */
		linadd((DiffuseInward+SourceReflected+DiffuseOutward+SourceTransmitted)/radius.dVeff,
			TauDummy.WLAng,"nFnu",'i' );

		/* emslin saves the per unit vol emissivity of a line, which is normally 
		 * what goes into linadd.  We zero this unit emissivity which was set
		 * FOR THE PREVIOUS LINE since it is so situation dependent */
		if( KILL_CONT && LineSave.ipass > 0 )
			LineSv[LineSave.nsum-1].emslin = 0.;

		/* this is the normal set to zero to trick the NEXT line into going in properly */
		if( LineSave.ipass > 0 )
			LineSv[LineSave.nsum].sumlin = 0.;

		/* the nsum-1 -- emslin and nsum -- sumlin is not a bug, look above - they do
		 * different things to different saves */

		Cont_nInu = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq +
			rfield.ConRefIncid[ipPredCont[i]]*renorm;
		/* transmitted and reflected indicent continuum, at selected energy points */
		linadd( Cont_nInu/radius.dVeff,
			TauDummy.WLAng,"nInu",'i' );

		/* emslin saves the per unit vo emissivity of a line, which isnormally 
		 * what goes into linadd.  We zero this unit emissivity since it is so situation dependent */
		if( KILL_CONT && LineSave.ipass > 0 )
			LineSv[LineSave.nsum-1].emslin = 0.;

		/* memory has not been allocated until ipass >= 0 */
		if( LineSave.ipass > 0 )
			LineSv[LineSave.nsum].sumlin = 0.;

		/* >>chng 00 dec 02, did not include SourceReflected */
		/* total reflected continuum, total inward emission plus
		 * reflected (XXdiffuseXX) total continuum */
		linadd( (DiffuseInward+SourceReflected)/radius.dVeff,TauDummy.WLAng,"InwT",'i');

		if( KILL_CONT && LineSave.ipass > 0 )
			LineSv[LineSave.nsum-1].emslin = 0.;

		/* memory has not been allocated until ipass >= 0 */
		if( LineSave.ipass > 0 )
			LineSv[LineSave.nsum].sumlin = 0.;

		/* reflected incident continuum (only incident) */
		linadd(SourceReflected/radius.dVeff,TauDummy.WLAng,"InwC",'i');

		if( KILL_CONT && LineSave.ipass > 0 )
			LineSv[LineSave.nsum-1].emslin = 0.;
	}

#	ifdef DEBUG_FUN
	fputs( " <->LineSet1()\n", debug_fp );
#	endif
	return;
}
#undef	KILL_CONT

