/* 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 */
/*RT_diffuse evaluate local diffuse emission for this zone,
 * fill in ConEmitLocal array with diffuse emission,
 * OTS rates for this zone were set in RT_OTS,
 * called by Cloudy, this routine adds energy to the outward beam */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "grains.h"
#include "grainvar.h"
#include "iso.h"
#include "dense.h"
#include "opacity.h"
#include "trace.h"
#include "rfield.h"
#include "phycon.h"
#include "hmi.h"
#include "radius.h"
#include "atmdat.h"
#include "heavy.h"
#include "atomfeii.h"
#include "lines_service.h"
#include "h2.h"
/* only used during debug print */
#include "ipoint.h"
#include "rt.h"

#define	TwoS	(1+ipISO)

void RT_diffuse(void)
{
	/* If this flag is set to TRUE, a table of 2 photon emission coefficients
	 * comparable to Brown and Mathews (1971) will be printed.	*/
	static long lgPrt2PhtEmisCoef = FALSE;

	/* DiffuseEscape is the spectrum of diffuse emission that escapes this zone,
	 * at end of this routine part is thrown into the outward beam -
	 * units are photons s-1 cm-3 */
	static float *DiffuseEscape/*[NC_ELL]*/ ;
	static int lgMustInit=TRUE;

	long int i=-100000, 
	  ip=-100000, 
	  ipISO=-100000,
	  ipHi=-100000, 
	  ipLo=-100000, 
	  ipla=-100000,
	  nelem=-100000,
	  ion=-100000,
	  limit=-100000, 
	  n=-100000,
	  nu=-10000;

	double EnergyLimit;

	double EdenAbund, 
	  difflya, 
	  xInWrd,
	  arg, 
	  fac, 
	  factor, 
	  gamma, 
	  gion, 
	  gn, 
	  photon, 
	  sum,
	  Sum1level,
	  SumCaseB;

	float Dilution , fach;

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

	if( lgAbort )
	{
		/* quit if we are aborting */
#		ifdef DEBUG_FUN
		fputs( " <->RT_diffuse()\n", debug_fp );
#		endif
		return;
	}

	/* one time creation of space for ThowOut array */
	if( lgMustInit )
	{
		lgMustInit = FALSE;
		if( (DiffuseEscape = (float*)MALLOC((unsigned)rfield.nupper*sizeof(float))  ) == NULL )
			BadMalloc();
	}

	/* many arrays were malloced to nupper, and we will add unit flux to [nflux] -
	 8 this must be true to work */
	ASSERT(rfield.nflux < rfield.nupper );

	/* this routine evaluates the local diffuse fields
	 * it fills in all of the following vectors */
	memset(DiffuseEscape                    , 0 , (unsigned)rfield.nupper*sizeof(float) );
	memset(rfield.ConEmitLocal         , 0 , (unsigned)rfield.nupper*sizeof(float) );
	memset(rfield.ConOTS_local_photons , 0 , (unsigned)rfield.nupper*sizeof(float) );
	memset(rfield.TotDiff2Pht          , 0 , (unsigned)rfield.nupper*sizeof(float) );

	/* loop over iso-sequences of all elements 
	 * to add all recombination continua and lines*/
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		/* >>chng 01 sep 23, rewrote for iso sequences */
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			/* this will be the sum of recombinations to all excited levels */
			SumCaseB = 0.;

			/* the product of the densities of the parent ion and electrons */
			EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];

			/* recombination continua for all iso seq - 
			 * if this stage of ionization exists */
			if( dense.IonHigh[nelem] >= nelem+1-ipISO  )
			{
				/* loop over all levels to include recombination diffuse continua,
				 * pick highest energy continuum point that opacities extend to 
				 * for ground continuum, go up to highest defined Boltzmann factor,
				 * at bottom of loop will be reset to ground state photo edge */
				ipHi = rfield.nflux;
				for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
				{
					Sum1level = 0.;

					/* >>chng 02 nov 20 - pull the plug on old he treatment */
					/* the number is (2 pi me k/h^2) ^ -3/2 * 8 pi/c^2 / ge - it includes
					 * the stat weight of the free electron in the demominator */
					/* TODO	2	This doesn't really seem to be the expression above!!!	*/
					/*gamma = 2.0618684e11*iso.stat[ipISO][nelem][n]/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte;*/
					gamma = 0.5*MILNE_CONST*iso.stat[ipISO][nelem][n]/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte;

					/* loop over all recombination continua 
					 * escaping part of recombinations are added to rfield.ConEmitLocal 
					 * added to ConInterOut at end of routine */
					for( nu=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; nu < ipHi; nu++ )
					{
						/* dwid used to adjust where within WIDFLX exp is evaluated -
						 * weighted to lower energy due to exp(-energy/T) */
						double dwid = 0.2;

						/* this is the term in the negative exponential Boltzmann factor
						 * for continuum emission */
						arg = (rfield.anu[nu]-iso.xIsoLevNIonRyd[ipISO][nelem][n]+
							rfield.widflx[nu]*dwid)/phycon.te_ryd;
						arg = MAX2(0.,arg);
						/* this is the limit sexp normally uses Boltzann */
						/*TODO	2	Steve Kraemer says we're lacking radio recombination continua.
						 * Is this why?  --RLP */
						if( arg > SEXP_LIMIT ) 
							break;

						/* photons is photons cm3 s-1 per cell */
						photon = gamma*exp(-arg)*rfield.widflx[nu]*
							opac.OpacStack[nu-iso.ipIsoLevNIonCon[ipISO][nelem][n] + 
							iso.ipOpac[ipISO][nelem][n]]*rfield.anu2[nu];
						/*ASSERT( photon >= 0. );*/
						Sum1level += photon;
						/* total local diffuse emission */
						rfield.ConEmitLocal[nu] += (float)(photon*EdenAbund);
						/* iso.RadRecomb[ipISO][nelem][n][ipRecEsc] is escape probability
						 * DiffuseEscape is local emission that escapes this zone */
						DiffuseEscape[nu] += 
							(float)(photon*EdenAbund*iso.RadRecomb[ipISO][nelem][n][ipRecEsc]);
					}
					/* this will be used below to confirm case B sum */
					if( n > 0 )
					{
						/* SumCaseB will be sum to all excited */
						SumCaseB += Sum1level;
					}

					/* on entry to this loop ipHi was set to the upper limit of the code,
					 * and ground state recom continua for all energies was added.  For
					 * excited states (which are the ones tha will be done after
					 * first pass through) only go to ground state threshold,
					 * since that will be so much stronger than excited state recom*/
					ipHi = iso.ipIsoLevNIonCon[ipISO][nelem][0]-1;
				}
				/* this is check on self-consistency */
				iso.CaseBCheck[ipISO][nelem] = MAX2(iso.CaseBCheck[ipISO][nelem],
				  (float)(SumCaseB/iso.RadRec_caseB[ipISO][nelem]));

				/* this add line emission from the model atoms,
				 * do not include very highest level since disturbed by topoff */
				for( ipLo=0; ipLo < (iso.numLevels_max[ipISO][nelem] - 2); ipLo++ )
				{
					for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem] - 1; ipHi++ )
					{
						/* must not include fake transtions (2s-1s, two photon, or truely
						 * forbidden transitions */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont < 1 )
							continue;

						/* pointer to line energy in continuum array */
						ip = EmisLines[ipISO][nelem][ipHi][ipLo].ipCont-1;

						/* number of photons in the line has not been defined up til now,
						 * do so now.  this is redone in lines.  */
						EmisLines[ipISO][nelem][ipHi][ipLo].phots = 
							EmisLines[ipISO][nelem][ipHi][ipLo].Aul*
							EmisLines[ipISO][nelem][ipHi][ipLo].PopHi*
							EmisLines[ipISO][nelem][ipHi][ipLo].Pesc*
							dense.xIonDense[nelem][nelem+1-ipISO];

						/* inward fraction of line */
						xInWrd = EmisLines[ipISO][nelem][ipHi][ipLo].phots*
							EmisLines[ipISO][nelem][ipHi][ipLo].FracInwd;

						/* reflected part of line */
						rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn);

						/* inward beam that goes out since sphere set */
						/* in all this the ColOvTot term has been commented out,
						 * since this is not meaningful for something like the hydrogen atom,
						 * where most forms by recombination */
						rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]/*
						EmisLines[ipISO][nelem][ipHi][ipLo].ColOvTot*/);

						/* outward part of line */
						rfield.outlin[ip] += 
							(float)(EmisLines[ipISO][nelem][ipHi][ipLo].phots*
							(1. - EmisLines[ipISO][nelem][ipHi][ipLo].FracInwd)*
							radius.BeamOutOut*opac.tmn[ip]/*
							EmisLines[ipISO][nelem][ipHi][ipLo].ColOvTot*/);
					}
				}

				/*Iso treatment of two photon emission.  */
				/* NISO could in the future be increased, but we want this assert to blow
				 * so that it is understood this may not be correct for other iso sequences,
				 * probably should break since will not be present */
				ASSERT( ipISO <= ipHE_LIKE );
					
				/* upper limit to 2-phot is energy of 2s to ground */
				EnergyLimit = EmisLines[ipISO][nelem][TwoS][0].EnergyWN/RYD_INF;
				
				/* this could fail when pops very low and Pop2Ion is zero */
				ASSERT( iso.ipTwoPhoE[ipISO][nelem]>0 && EnergyLimit>0. );
				
				sum = 0.;
				/* two photon emission, iso.ipTwoPhoE[ipISO][nelem] is 
				 * continuum array index for Lya energy */
				for( nu=0; nu<iso.ipTwoPhoE[ipISO][nelem]; nu++ )
				{
					/* We do not assert rfield.anu[nu]<=EnergyLimit because the maximum 
					 * index could be set to point to the next highest bin.	*/
					ASSERT( rfield.anu[nu]/EnergyLimit < 1.01 || rfield.anu[nu-1]<EnergyLimit);

					/* iso.As2nu[ipISO][nelem][nu] is transition probability A per bin
					 * So sum is the total transition probability - this sum should
					 * add up to the A two photon */
					sum += iso.As2nu[ipISO][nelem][nu];

					/* flag saying whether to also do he triplets */
#					define lgDo2TripSToo	FALSE

					if( ipISO == ipHE_LIKE && lgDo2TripSToo )
					{
						fach = 2.f*dense.xIonDense[nelem][nelem+1-ipISO]*iso.As2nu[ipISO][nelem][nu]*
							((float)iso.Pop2Ion[ipISO][nelem][TwoS]+0.0001f*(float)iso.Pop2Ion[ipISO][nelem][ipHe2s3S]);
					}
					else
					{
						/* iso.Pop2Ion[ipISO][nelem][TwoS] is dimensionless and represents the population 
						 * of the TwoS level relative to that of the ion.  dense.xIonDense[nelem][nelem+1-ipISO]
						 * is the density of the current ion (cm^-3).  The factor of 2 is for two photons per 
						 * transition. Thus fach has dimension photons cm-3 s-1.	*/
						fach = 2.f*(float)iso.Pop2Ion[ipISO][nelem][TwoS]*dense.xIonDense[nelem][nelem+1-ipISO]
							*iso.As2nu[ipISO][nelem][nu];
					}

					/* >>chng 03 mar 26, only do if induced processes turned on,
					 * otherwise inconsistent with rate solver treatment.	*/
					/* >>chng 02 aug 14, add induced two photon emission */
					/* product of occupation numbers for induced emission */
					if( iso.lgInd2nu_On)
					{
						/* this is the total rate */
						fach *= (1.f + rfield.SummedOcc[nu]) *
							(1.f+rfield.SummedOcc[iso.ipSym2nu[ipISO][nelem][nu]-1]);
					}

					rfield.TotDiff2Pht[nu] += fach;

					/* total local diffuse emission */
					rfield.ConEmitLocal[nu] += fach;

					/* this is escaping part, as determined from optical depth to illuminated face */
					DiffuseEscape[nu] += fach*opac.ExpmTau[nu];

					/* save locally destroyed OTS two-photon continuum */
					rfield.ConOTS_local_photons[nu] += fach*(1.f - opac.ExpmTau[nu]);
				}

				/* 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 - (float)sum/EmisLines[ipISO][nelem][TwoS][0].Aul ) < 0.01f );
			}

			/* option to print hydrogen and helium two-photon emission coefficients?	*/
			if( lgPrt2PhtEmisCoef )
			{
				long yTimes20;
				double y, E2nu;

				lgPrt2PhtEmisCoef = FALSE;

				fprintf( ioQQQ, "\ny\tGammaNot(2q)\n");

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

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

					E2nu = EmisLines[0][0][1][0].EnergyWN/RYD_INF;
					i = ipoint(y*E2nu);
					fprintf( ioQQQ, "%.3e\t", 
						8./3.*HPLANCK*iso.Pop2Ion[0][0][1]/dense.eden
						*y*iso.As2nu[0][0][i]*E2nu/rfield.widflx[i] );
					
					E2nu = EmisLines[1][1][2][0].EnergyWN/RYD_INF;
					i = ipoint(y*E2nu);
					fprintf( ioQQQ, "%.3e\n", 
						8./3.*HPLANCK*iso.Pop2Ion[1][1][2]/dense.eden
						*y*iso.As2nu[1][1][i]*E2nu/rfield.widflx[i] );

					/*
					fprintf( ioQQQ, "%.3e\t%.3e\n", 
						rfield.TotDiff2Pht[i]*rfield.anu[i]*EN1RYD
						/E2nu/rfield.widflx[i]/dense.eden/dense.xIonDense[nelem][nelem+1-ipISO]/FR1RYD,
						8./3.*HPLANCK*iso.Pop2Ion[ipISO][nelem][TwoS]/dense.eden
						*y*iso.As2nu[ipISO][nelem][i]*E2nu/rfield.widflx[i]);	*/
				}
				fprintf( ioQQQ, "eden%.3e\n", dense.eden );
			}
		}
	}

	/* add recombination continua for elements heavier than those done with iso seq */
	for( nelem=NISO; nelem < LIMELM; nelem++ )
	{
		/* do not include species with iso-sequence in following */
		/* >>chng 03 sep 09, upper bound was wrong, did not include NISO */
		for( ion=dense.IonLow[nelem]; ion < nelem-NISO+1; ion++ )
		{
			if( dense.xIonDense[nelem][ion+1] > 0. )
			{
				long int ns, nshell,igRec , igIon,
					iplow , iphi , ipop;

				ip = Heavy.ipHeavy[nelem][ion]-1;

				/* nflux was reset upward in ConvInitSolution to encompass all
				 * possible line and continuum emission.  this test should not
				 * possibly fail.  It could if the ionization were to increase with depth.
				 * This is important because the nflux cell in ConInterOut is used to carry out the
				 * unit integration, and if it gets clobbered by diffuse emission the code
				 * will declare insanity in PrtComment */
				ASSERT( ip >= 0 && ip < rfield.nflux );

				/* get shell number, stat weights for this species */
				atmdat_outer_shell( nelem+1 , nelem+1-ion , &nshell, &igRec , &igIon );
				gn = (double)igRec;
				gion = (double)igIon;

				/* shell number */
				ns = Heavy.nsShells[nelem][ion]-1;
				ASSERT( ns == (nshell-1) );

				/* lower and upper energies, and offset for opacity stack */
				iplow = opac.ipElement[nelem][ion][ns][0]-1;
				iphi = opac.ipElement[nelem][ion][ns][1];
				iphi = MIN2( iphi , rfield.nflux );
				ipop = opac.ipElement[nelem][ion][ns][2];

				/* now convert ipop to the offset in the opacity stack from threshold */
				ipop = ipop - iplow;

				gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte*dense.eden*dense.xIonDense[nelem][ion+1];
				
				/* this is ground state continuum from stored opacities */
				if( rfield.ContBoltz[iplow] > SMALLFLOAT )
				{
					for( nu=iplow; nu < iphi; ++nu )
					{
						photon = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[iplow]*
							rfield.widflx[nu]*opac.OpacStack[nu+ipop]*rfield.anu2[nu];
						/* add heavy rec to ground in active beam,
						 * TODO	2	should use ConEmitLocal for all continua, but not followed
						 * by DiffuseEscape - put that at the end.  Once continua all
						 * bundled this way, it will be easy to save them as a function
						 * of depth and then do exact rt */
						rfield.ConEmitLocal[nu] += (float)photon;
						/*DiffuseEscape[i] += (float)photon;*/
						/*>>chng 03 sep 08, index had been i, should have been nu */
						DiffuseEscape[nu] += (float)photon*opac.ExpmTau[nu];
					}
				}

				/* now do the recombination Lya */
				ipla = Heavy.ipLyHeavy[nelem][ion]-1;
				ASSERT( ipla >= 0 );
				/* xLyaHeavy is set to a fraction of the total rad rec in ion_recomb, includes eden */
				difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
				/* >>chng 03 jul 10, here and below, use outlin_noplot */
				rfield.outlin_noplot[ipla] += (float)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);

				/* now do the recombination Balmer photons */
				ipla = Heavy.ipBalHeavy[nelem][ion]-1;
				ASSERT( ipla >= 0 );
				/* xLyaHeavy is set to fraction of total rad rec in ion_recomb, includes eden */
				difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
				rfield.outlin_noplot[ipla] += (float)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
			}
		}
	}

	/* free-free free free brems emission for all ions */
	limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
	for( nu=0; nu < limit; nu++ )
	{
		double TotBremsAllIons = 0., BremsThisIon;

		/* do hydrogen first, before main loop since want to add on H- brems */
		nelem = ipHYDROGEN;
		ion = 1;
		BremsThisIon = POW2( (float)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
		ASSERT( BremsThisIon >= 0. );
		
		/* for case of hydrogen, add H- brems - OpacStack contains the ratio
		 * of the H- to H brems cross section - multiply by this and the fraction in ground */
		BremsThisIon *= (1.+opac.OpacStack[nu-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
		TotBremsAllIons = BremsThisIon;

		/* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
		 * using gaunt factors from rfield.gff.	*/
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			/* MAX2 occurs because we want to start at first ion (or above)
			 * and do not want atom */
			for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
			{
				/* eff. charge is ion, so first rfield.gff argument must be "ion".	*/
				BremsThisIon = POW2( (float)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
				/*ASSERT( BremsThisIon >= 0. );*/

				TotBremsAllIons += BremsThisIon;
			}
		}

		/*TODO	2	Replace this constant with the appropriate macro, if any */
		TotBremsAllIons *= dense.eden*1.032e-11*rfield.widflx[nu]*rfield.ContBoltz[nu]/rfield.anu[nu]/phycon.sqrte;
		ASSERT( TotBremsAllIons >= 0.);
		
		/* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
		 * to outward beam - ConLocNoInter array removed as result
		 * if problems develop with very dense blr clouds, this may be reason */
		/*rfield.ConLocNoInter[nu] += (float)fac;*/
		/*rfield.ConEmitLocal[nu] += (float)TotBremsAllIons;*/

		if( nu >= rfield.ipEnergyBremsThin )
		{
			/* >>chng 05 feb 20, move into this test on brems opacity - should not be
			 * needed since would use expmtau to limit outward beam */
			/* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
			* to outward beam - ConLocNoInter array removed as result
			* if problems develop with very dense blr clouds, this may be reason */
			/*rfield.ConLocNoInter[nu] += (float)fac;*/
			rfield.ConEmitLocal[nu] += (float)TotBremsAllIons;
			
			/* do not add optically thick part to outward beam since self absorbed
			* >>chng 96 feb 27, put back into outward beam since do not integrate
			* over it anyway. */
			/* >>chng 99 may 28, take back out of beam since DO integrate over it
			* in very dense BLR clouds */
			/* >>chng 01 jul 10, add here, in only one loop, where optically thin */
			DiffuseEscape[nu] += (float)TotBremsAllIons;
		}
	}

	/* grain dust emission */
	/* >>chng 01 nov 22, moved calculation of grain flux to qheat.c, PvH */
	if( gv.lgDustOn && gv.lgGrainPhysicsOn )
	{
		/* this calculates diffuse emission from grains,
		 * and stores the result in gv.GrainEmission */
		GrainMakeDiffuse();

		for( nu=0; nu < rfield.nflux; nu++ )
		{
			rfield.ConEmitLocal[nu] += gv.GrainEmission[nu];
			DiffuseEscape[nu] += gv.GrainEmission[nu];
		}
	}

	/* hminus emission */
	fac = dense.eden*(double)dense.xIonDense[ipHYDROGEN][0];
	gn = 1.;
	gion = 2.;
	gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
	/* >>chng 00 dec 15 change limit to -1 of H edge */
	limit = MIN2(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,rfield.nflux);

	if( rfield.ContBoltz[hmi.iphmin-1] > 0. )
	{
		for( nu=hmi.iphmin-1; nu < limit; nu++ )
		{
			/* H- flux photons cm-3 s-1 
			 * ContBoltz is ratio of Boltz fac for each freq */
			factor = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[hmi.iphmin-1]*rfield.widflx[nu]*
			  opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
			  rfield.anu2[nu]*fac;
			rfield.ConEmitLocal[nu] += (float)factor;
			DiffuseEscape[nu] += (float)factor;
		}
	}
	else
	{
		for( nu=hmi.iphmin-1; nu < limit; nu++ )
		{
			arg = MAX2(0.,TE1RYD*(rfield.anu[nu]-0.05544)/phycon.te);
			/* this is the limit sexp normally uses */
			if( arg > SEXP_LIMIT ) 
				break;
			/* flux is in photons per sec per ryd */
			factor = gamma*exp(-arg)*rfield.widflx[nu]*
				opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
			  rfield.anu2[nu]*fac;
			rfield.ConEmitLocal[nu] += (float)factor;
			DiffuseEscape[nu] += (float)factor;
		}
	}

	/* this dilution is needed to conserve volume in spherical models.  tests such
	 * as parispn.in will fault if this is removed */
	Dilution = (float)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[rfield.nflux] = 1.e-10f * Dilution;
	DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;

	/* opacity should be zero at this energy */
	ASSERT( opac.opacity_abs[rfield.nflux] == 0. );

	/* many tmn added to conserve energy in very high Z models
	 * rerun highZ qso model if any tmn ever deleted here
	 *
	 * tmn set in ZoneStart and includes both opacity and dilution 
	 *
	 * use duffuse lines and continuum to add flux to outward beam
	 *
	 * NB!!!  this routine adds flux to the outward beam
	 *  it can only be called once per zone
	 */
	 
	/* >>chng 96 nov 19, do not consider energies below plasma freq
	 * ipPlasmaFreq points to lowest energy thin to brems and plas freq */

	/* add DiffuseEscape continuum to ConInterOut,
	 * lower limit of rfield.ipPlasma-1 since continuum is zero below rfield.ipPlasma-1 
	 * due to plasma frequency
	 * note that upper range of sum is <= nflux, which contains the unit
	 * verification token in the highest cell*/
	for( nu=rfield.ipPlasma-1; nu <= rfield.nflux; nu++ )
	{
		/* ConInterOut is the interactive continuum
		 * DiffuseEscape contains all radiation thrown into outward beam */
		/* radius.dVolOutwrd has factor or 1 or 0.5 for outward part */
		/* >>chng 05 feb 23, activate this statement */
		rfield.ConInterOut[nu] += DiffuseEscape[nu]*(float)radius.dVolOutwrd;
		/*rfield.ConInterOut[nu] += rfield.ConEmitLocal[nu]*opac.ExpmTau[nu]*(float)radius.dVolOutwrd;*/
		ASSERT( rfield.ConInterOut[nu] >= 0.);
	}

	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			fprintf(ioQQQ,"rtdiffusedebugg %li\t%.2e\t%.2e\t%.2e\n", 
				nzone, rfield.ConInterOut[1158] , DiffuseEscape[1158],radius.dVolOutwrd);
		}
	}

	/* outward level 1 line photons, 0 is dummy line */
	for( i=1; i <= nLevel1; i++ )
	{
		outline( &TauLines[i] );
	}

	/* outward level 2 line photons */
	for( i=0; i < nWindLine; i++ )
	{
		/* must not also do lines that were already done as part
		 * of the isoelectronic sequences */
		if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
		{
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC /*&& nzone > 10*/ && i==4821 )
				{
					/* set up to dump the Fe 9 169A line */
					fprintf(ioQQQ,"DEBUG dump lev2 line %li\n", i );
					DumpLine( &TauLine2[i] );/**/
					fprintf(ioQQQ,"DEBUG dump %.3e %.3e %.3e\n",
						rfield.outlin[TauLine2[i].ipCont-1],
						TauLine2[i].phots*TauLine2[i].FracInwd*radius.BeamInOut*opac.tmn[i]*TauLine2[i].ColOvTot,
						TauLine2[i].phots*(1. - TauLine2[i].FracInwd)*radius.BeamOutOut* TauLine2[i].ColOvTot );
				}
			}
			outline( &TauLine2[i] );
			/*if( i==2576 ) fprintf(ioQQQ,"DEBUG dump %.3e %.3e \n",
				rfield.outlin[TauLine2[i].ipCont-1] , rfield.outlin_noplot[TauLine2[i].ipCont-1]);*/
		}
	}

	/* outward hyperfine structure line photons */
	for( i=0; i < nHFLines; i++ )
	{
		outline( &HFLines[i] );
	}

	/* carbon monoxide */
	for( i=0; i < nCORotate; i++ )
	{
		outline( &C12O16Rotate[i] );
		outline( &C13O16Rotate[i] );
	}

	/* H2 emission */
	H2_RT_diffuse();

	/* do outward parts of FeII lines, if large atom is turned on */
	FeII_RTOut();
	/*TODO	2	add fegrain to outward beams, but within main formalism by including grains
	 * in all x-ray processes */

	if( trace.lgTrace )
		fprintf( ioQQQ, " RT_diffuse returns.\n" );

	/* >>chng 02 jul 25, zero out all light below plasma freq */
	for( nu=0; nu<rfield.ipPlasma-1; nu++ )
	{
		rfield.flux[nu] = 0.;
		rfield.ConEmitLocal[nu] = 0.;
		rfield.otscon[nu] =0.;
		rfield.otslin[nu] =0.;
		rfield.outlin[nu] =0.;
		rfield.outlin_noplot[nu] =0.;
		rfield.reflin[nu] = 0.;
		rfield.TotDiff2Pht[nu] = 0.;
		rfield.ConOTS_local_photons[nu] = 0.;
		rfield.ConInterOut[nu] = 0.;
	}

	/* find occupation number, also assert that no continua are negative */
	for( nu=0; nu < rfield.nflux; nu++ )
	{
		/* >>chng 00 oct 03, add diffuse continua */
		/* local diffuse continua */
		rfield.OccNumbDiffCont[nu] = rfield.ConEmitLocal[nu]*rfield.convoc[nu];

		/* confirm that all are non-negative */
		ASSERT( rfield.flux[nu] >=0.);
		ASSERT( rfield.otscon[nu] >=0.) ;
		ASSERT( rfield.otslin[nu] >=0.) ;
		ASSERT( rfield.ConInterOut[nu] >=0.) ;
		ASSERT( rfield.outlin[nu] >=0.) ;
	}

	/* option to kill outward lines with no outward lines command*/
	if( rfield.lgKillOutLine )
	{
		for( nu=0; nu < rfield.nflux; nu++ )
		{
			rfield.outlin[nu] = 0.;
			rfield.outlin_noplot[nu] = 0.;
		}
	}

	/* option to kill outward continua with no outward continua command*/
	if( rfield.lgKillOutCont )
	{
		for( nu=0; nu < rfield.nflux; nu++ )
		{
			rfield.ConInterOut[nu] = 0.;
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->RT_diffuse()\n", debug_fp );
#	endif
	return;
}
#undef	TwoS
#undef lgDo2TripSToo

