/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*RT_OTS compute diffuse fields due to H, He atoms, ion, triplets, metal recombination,
 * called by ConvBase  */
/*RT_OTS_AddLine add local destruction of lines to ots field */
/*RT_OTS_AddCont add local destruction of continuum to ots field */
/*RT_OTS_Update sum flux, otscon, otslin, ConInterOut, outlin, to form SummeDif, SummedCon SummedOcc */
/*RT_OTS_Zero - zero out some vectors - 
 * this is only called when code initialized by ContSetIntensity */
/*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
#include "cddefines.h"
#include "taulines.h"
#include "opacity.h"
#include "dense.h"
#include "iso.h"
#include "hmi.h"
#include "h2.h"
#include "rfield.h"
#include "converge.h"
#include "rt.h"
#include "co.h"
#include "atomfeii.h"
#include "heavy.h"
#include "he.h"
#include "trace.h"

/* this flag may be used for debugging ots rate changes */
static int nOTS_Line_type = -1;
static int nOTS1=-1 , nOTS2=-1;
/*add local destruction of continuum to ots field */
static void RT_OTS_AddCont(
	/* the ots rate itself */
	float ots, 
	/* pointer to continuum cell for ots, on f scale */
	long int ip);

/* =================================================================== */

void RT_OTS(void)
{
	long int
		ipla,
		ipISO ,
		nelem,
		n;
	float 
		difflya,
		esc,
		ots;

	/* the Bowen HeII yield */
#	define BOWEN 0.3f
	long int ipHi, 
	  ipLo;

	float bwnfac ,
	  ots660 ,
	  cont_phot_destroyed,
	  save_lya_dest,
	  save_he2lya_dest;

	double save_he2rec_dest;

	/*static long int nCall=0;
	++nCall;
	fprintf(ioQQQ,"debugggtos enter %li\n", nCall );*/

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

	/**************************************************************************
	 *
	 * the bowen HeII - OIII fluorescense problem
	 *
	 **************************************************************************/
	nOTS_Line_type = 0;
	nelem = ipHELIUM;
	if( dense.lgElmtOn[nelem] )
	{
		/* conversion per unit atom to OIII, at end of sub we divide by it,
		 * to fix lines back to proper esc/dest probs */
		bwnfac = BOWEN * MAX2(0.f,1.f- EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pesc - 
			EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pelec_esc );

		/* the last factor accounts for fact that two photons are produced,
		 * and the branching ratio */
		ots660 = (float)(EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Aul*
		  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopHi*dense.xIonDense[nelem][nelem+1]*
		  EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest *BOWEN*2.0 * 0.8 );

		/* now add this to the ots field */
		if( ots660 > SMALLFLOAT ) 
			RT_OTS_AddLine(ots660 , he.ip660 );

		/* decrease the destruction prob by the amount we will add elsewhere,
		 * ok since dest probs updated on every iteration*/
		EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest *= bwnfac;
		ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest >= 0. );
	}

	else
	{
		bwnfac = 1.;
	}

	/* save Lya loss rates so we can reset at end */
	save_lya_dest = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest;

	/* >>chng 03 may 28, hydro.dstfe2lya should not go into the ots field since it was
	 * actually lost in exciting FeII 
	 * this is not ready to be used - need to recover ots with small feii atom
	fprintf(ioQQQ,"feii debug\t%.3e\t%.3e\n", 
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , hydro.dstfe2lya);
	EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = MAX2(0.f,
		EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest-hydro.dstfe2lya);*/

	/* this is option to kill Lya ots rates, 
	 * rfield.lgLyaOTS is usually TRUE (1), and set FALSE (0) with
	 * no lya ots command */
	EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest *= rfield.lgLyaOTS;

	/* option to kill heii lya and rec continua ots */
	save_he2lya_dest = EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest;
	EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest *= rfield.lgHeIIOTS;

	/* option to kill heii lya and rec continua ots */
	save_he2rec_dest = iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad];
	iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] *= rfield.lgHeIIOTS;

	nOTS_Line_type = 1;

	/* make ots fields due to lines and continua of species treated with unified 
	 * isoelectronic sequence */
	/* loop over all elements */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			nOTS1 = ipISO;
			nOTS2 = nelem;
			/* do if this stage exists */
			if( (dense.IonHigh[nelem] >= nelem+1-ipISO )  )
			{
				/* generate line ots rates */
				/* now loop over all possible levels, but cannot include two photon
				* since there is no pointer to this continuum */
				for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ )
				{
					for( ipLo=0; ipLo < ipHi; ipLo++ )
					{
						/* this signifies a fake line */
						/* >>chng 03 may 19, DEST0 is the smallest possible
						 * dest prob - not a real number, don't add to ots field */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont< 1 ||
							EmisLines[ipISO][nelem][ipHi][ipLo].Pdest<= DEST0 )
							continue;

						/* ots rates, the destp prob was set in hydropesc */
						EmisLines[ipISO][nelem][ipHi][ipLo].ots = (float)(
							EmisLines[ipISO][nelem][ipHi][ipLo].Aul*
							EmisLines[ipISO][nelem][ipHi][ipLo].PopHi*dense.xIonDense[nelem][nelem+1-ipISO]*
							EmisLines[ipISO][nelem][ipHi][ipLo].Pdest);

						ASSERT( EmisLines[ipISO][nelem][ipHi][ipLo].ots >= 0. );
						/* way to kill lyalpha ots rates
						if( nelem==ipHYDROGEN && ipHi==ipH2p && ipLo==ipH1s )
						{
							EmisLines[ipISO][nelem][ipHi][ipLo].ots = 0.;
						} */

						/* finally dump the ots rate into the stack */
						if( EmisLines[ipISO][nelem][ipHi][ipLo].ots > SMALLFLOAT ) 
							RT_OTS_AddLine(EmisLines[ipISO][nelem][ipHi][ipLo].ots,
								EmisLines[ipISO][nelem][ipHi][ipLo].ipCont );
					}
				}
				{
					/* debugging code for line oscillation problems */
					/*@-redef@*/
					enum {DEBUG_LOC=FALSE};
					/*@+redef@*/
					if( DEBUG_LOC )
					{
						long ip;
						if( ipISO==0 && nelem==0  && nzone>90  ) 
						{
							ipHi = 2;
							ipLo = 0;
							ip = EmisLines[ipISO][nelem][ipHi][ipLo].ipCont;
							fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
								fnzone, 
								EmisLines[ipISO][nelem][ipHi][ipLo].ots,
								opac.opacity_abs[ip-1],
								EmisLines[ipISO][nelem][ipHi][ipLo].PopHi*dense.xIonDense[nelem][nelem+1-ipISO],
								EmisLines[ipISO][nelem][ipHi][ipLo].PopHi,
								dense.xIonDense[nelem][nelem+1-ipISO],
								EmisLines[ipISO][nelem][ipHi][ipLo].Pdest,
								rfield.otslinNew[ip-1],
								rfield.otslin[ip-1]);
						}
					}
				}

				/**************************************************************************
				*
				* ots recombination bound-free b-f continua continuum
				*
				**************************************************************************/

				/* put in OTS continuum */
				for( n=0; n < iso.numLevels[ipISO][nelem]; n++ )
				{
					cont_phot_destroyed = (float)(iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
						(1. - iso.RadRecomb[ipISO][nelem][n][ipRecEsc])*dense.eden*
						dense.xIonDense[nelem][nelem+1]);
					ASSERT( cont_phot_destroyed >= 0. );

					/* continuum energy index used in this routine is decremented by one there */
					RT_OTS_AddCont(cont_phot_destroyed,iso.ipIsoLevNIonCon[ipISO][nelem][n]);
					/* debugging code for rec continua */
					{
						/*@-redef@*/
						enum {DEBUG_LOC=FALSE};
						/*@+redef@*/
						if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 )
						{
							long ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;
							fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
								fnzone, 
								n ,
								EmisLines[ipISO][ipHYDROGEN][2][0].PopHi*dense.xIonDense[ipHYDROGEN][1],
								cont_phot_destroyed,
								cont_phot_destroyed/opac.opacity_abs[ip],
								rfield.otsconNew[ip] ,
								rfield.otscon[ip] ,
								opac.opacity_abs[ip] ,
								hmi.Hmolec[ipMHm] ,
								hmi.HMinus_photo_rate);
						}
					}
				}
			}
		}
	}
	/* more debugging code for rec continua */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			nelem = 0;
			fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
				nzone, 
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][0]-1],
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][1]-1],
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][3]-1],
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][4]-1],
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][5]-1],
				rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1],
				opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1]);
		}
	}

	/* now reset Lya dest prob in case is was clobbered */
	EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = save_lya_dest;
	EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest = save_he2lya_dest;
	iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] = save_he2rec_dest;

	nelem = ipHELIUM;
	if( dense.lgElmtOn[nelem] && bwnfac > 0. )
	{
		/* increase the destruction prob by the amount we decreased it above */
		EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest /= bwnfac;
	}

	if( trace.lgTrace )
	{
		fprintf(ioQQQ,"     RT_OTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n",
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ots , 
			rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1]  ,
			opac.opacity_abs[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1]
			);
	}

	nOTS_Line_type = 2;
	/* recombination Lya for all elements not yet converted into std isoelectronc form */
	for( nelem=NISO ; nelem < LIMELM; nelem++ )
	{
		long int ion;
		/* do not include species treated in iso-electronic fashon in the following,
		 * these were treated above */
		for( ion=0; ion < nelem+1-NISO; ion++ )
		{
			if( dense.xIonDense[nelem][ion+1] > 0. )
			{
				nOTS1 = nelem;
				nOTS2 = ion;
				/* now do the recombination Lya */
				ipla = Heavy.ipLyHeavy[nelem][ion];
				esc = opac.ExpmTau[ipla-1];
				/* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
				difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
				/* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
				 * setting the ots rates to zero */
				ots = difflya*MAX2(0.f,1.f-esc);
				/*if( nelem==6 && ion==2 )
					fprintf(ioQQQ," debugggnly\t %.2f\t%.2e\n",fnzone, ots );*/
				ASSERT( ots >= 0.);
				/*if( iteration == 2 && nzone>290 && ipla== 2339 )
					fprintf(ioQQQ,"recdebugg1 %.2e %li %li %.2e %.2e \n", 
					ots, nelem, ion,
					esc , dense.xIonDense[nelem][ion+1]);*/
				if( ots > SMALLFLOAT ) 
					RT_OTS_AddLine(ots,ipla);

				/* now do the recombination balmer lines */
				ipla = Heavy.ipBalHeavy[nelem][ion];
				esc = opac.ExpmTau[ipla-1];
				/* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
				difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
				/* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
				 * setting the ots rates to zero */
				ots = difflya*MAX2(0.f,1.f-esc);
				ASSERT( ots >= 0.);
				/*if( iteration == 2 &&nzone==294 && ipla== 2339 )
					fprintf(ioQQQ,"recdebugg2 %.2e %li %li\n", ots, nelem, ion );*/
				if( ots > SMALLFLOAT ) 
					RT_OTS_AddLine(ots,ipla);
			}
		}
	}

	nOTS_Line_type = 3;
	/* do OTS and outward parts of FeII lines, if large atom is turned on */
	FeII_OTS();

	nOTS_Line_type = 4;
	/* do the main set of level1 lines */
#	if 0
	{
#	include "lines_service.h"
	if( nzone>290 ) DumpLine( &TauLines[74] );
	}
#	endif
/*#	define	FRACNEW	1.0f
#	undef FRACNEW*/
	for( nOTS1=1; nOTS1 < nLevel1; nOTS1++ )
	{
		/*float otsold = TauLines[nOTS1].ots;*/
		TauLines[nOTS1].ots = (float)TauLines[nOTS1].PopHi * TauLines[nOTS1].Aul * TauLines[nOTS1].Pdest;
		/* * TauLines[nOTS1].ColOvTot;*/
		/*if( nzone )
			TauLines[nOTS1].ots = TauLines[nOTS1].ots*FRACNEW + otsold*(1.f-FRACNEW);*/
		if( TauLines[nOTS1].ots > SMALLFLOAT ) 
			RT_OTS_AddLine( TauLines[nOTS1].ots , TauLines[nOTS1].ipCont);
	}

	nOTS_Line_type = 5;
	/* do the level2 level 2 lines */
	for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ )
	{
		if( TauLine2[nOTS1].IonStg < TauLine2[nOTS1].nelem+1-NISO )
		{
			TauLine2[nOTS1].ots = (float)(TauLine2[nOTS1].PopHi * TauLine2[nOTS1].Aul * TauLine2[nOTS1].Pdest);
			if( TauLine2[nOTS1].ots > SMALLFLOAT ) 
				RT_OTS_AddLine( TauLine2[nOTS1].ots , TauLine2[nOTS1].ipCont);
		}
	}
	/*fprintf(ioQQQ,"debuggne1 \t%.2e\t%.2e\t%.2e\t%.2e\n", 
		TauLine2[743].ots/SDIV(opac.opacity_abs[743]), 
		TauLine2[744].ots/SDIV(opac.opacity_abs[744]) ,
		TauLine2[743].ots,opac.opacity_abs[743]);*/
	/*CO_OTS - add CO lines to ots fields, called by RT_OTS */

	nOTS_Line_type = 6;
	CO_OTS();

	/* the large H2 molecule */

	nOTS_Line_type = 7;
	H2_RT_OTS();
	/*fprintf(ioQQQ,"debugggtos exit\n");*/

#	ifdef DEBUG_FUN
	fputs( " <->RT_OTS()\n", debug_fp );
#	endif

	return;
}

/* =================================================================== */

void RT_OTS_AddLine(float ots, 
	/* pointer on the f scale */
  long int ip )
{

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

	/* add ots due to line destruction to radiation field */

	/* return if outside bounds of this continuum source, ip > rfield.nflux
	 * first case ip==0 happens when called with dummy line */
	if( ip==0 || ip > rfield.nflux )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->RT_OTS_AddLine()\n", debug_fp );
#		endif
		return;
	}

	/*the local ots rate must be non-negative */
	ASSERT( ots >= 0. );
	/* continuum pointer must be positive */
	ASSERT( ip > 0 );

	/* add locally destroyed flux of photons to line OTS array */
	/* check whether local gas opacity  (units cm-1) is positive, if so
	 * convert line destruction rate into ots rate by dividing by it */
	if( opac.opacity_abs[ip-1] > 0. )
	{
		rfield.otslinNew[ip-1] += (float)(ots/opac.opacity_abs[ip-1]);
	}
	/* first iteration is 1, second is two */
	{
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && /*iteration == 2 && nzone>294 &&*/ ip== 2363 /*&& ots > 1e16*/)
		{
			fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t",
				ots ,
				opac.opacity_abs[ip-1],
				ots/opac.opacity_abs[ip-1] );
			fprintf(ioQQQ,"iteration %li type %i %i %i \n", 
				iteration, 
				nOTS_Line_type,
				nOTS1,nOTS2 );
		}
	}

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

/* =================================================================== */

/*add local destruction of continuum to ots field */
static void RT_OTS_AddCont(
					/* the ots rate itself */
					float ots, 
					/* pointer to continuum cell for ots, on f scale */
					long int ip)
{

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

	/* 
	 * routine called to add ots due to continuum destruction to
	 * radiation field
	 */
	 
	/* check if outside bounds of this continuum source */
	if( ip > rfield.nflux )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->RT_OTS_AddCont()\n", debug_fp );
#		endif
		return;
	}

	ASSERT( ip> 0 );

	ASSERT( ots >= 0. );

	ASSERT( ip <= rfield.nupper );

	/* add locally destroyed flux of photons to continuum OTS array */
	/* check whether local gas opacity  (units cm-1) is positive, if so
	 * convert continuum destruction rate into ots rate by dividing by it */
	if( opac.opacity_abs[ip-1] > 0. )
	{
		rfield.otsconNew[ip-1] += (float)(ots/opac.opacity_abs[ip-1]);
		/*if( ots>0. ) fprintf(ioQQQ,
			"buggg ip %li ots %.2e opac %.2e cont_phot_destroyed %.2e\n",
			ip, ots , opac.opacity_abs[ip-1] ,rfield.otsconNew[ip-1]);*/
	}

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

#if 0
/* check whether Lya OTS rate is oscillating */
static int lgLyaOTS_oscil( float ots_new )
{
	static float old , older;
	int lgReturn = FALSE;

	if( !conv.nTotalIoniz || conv.lgSearch )
	{
		/* this is very first call on this iteration, zero everything out */
		old = 0.;
		older = 0.;
	}

	/* is sign of change oscillating? */
	if( conv.nPres2Ioniz>1 && (ots_new-old) * (old-older) < 0.)
		lgReturn = TRUE;

	older = old;
	old = ots_new;
	return( lgReturn );
}

	/* >>chng 04 jan 26, check whether Lya is oscillating */
	if( lgLyaOTS_oscil( rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ) )
	{
		fprintf(ioQQQ,"DEBUG lya ots %.2f %.3e\n", 
			fnzone, 
			rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
	}
#endif
/* =================================================================== */

/*RT_OTS_Update update ots rates, called in ConvBase */
void RT_OTS_Update(
	/* summed ots rates */
	double *SumOTS , 
	/* array index for constituent that changed the most
	long int *ipOTSchange,  */
	/* limit on how large a change to allow, 0 for no limit */
	double BigFrac)
{
	long int i;

	static int lgNeedSpace=TRUE;
	static double *old_OTS_line_x_opac , *old_OTS_cont_x_opac;
	float FacBig , FacSml;

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

	if( BigFrac <= 0. )
		BigFrac = 10.;
	FacBig = (float)(1. + BigFrac);
	FacSml = (float)(1./FacBig);

	/* create space to save arrays if needed, one time per coreload */
	if( lgNeedSpace )
	{
		if( (old_OTS_line_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ) ) == NULL )
			BadMalloc();
		if( (old_OTS_cont_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) ) ) == NULL )
			BadMalloc();
	}
	lgNeedSpace = FALSE;

	*SumOTS = 0.;

	/* option to kill ots rates with no ots lines command */
	if( rfield.lgKillOTSLine )
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			rfield.otslinNew[i] = 0.;
		}
	}

	/* during search phase set ots to current values */
	if( nzone==0 )
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			old_OTS_line_x_opac[i] = rfield.otslinNew[i] * opac.opacity_abs[i];
			old_OTS_cont_x_opac[i] = rfield.otsconNew[i] * opac.opacity_abs[i];
		}
	}

	/* remember largest change in ots rates */
	*SumOTS = 0.;
	/* now update new ots rates */
	for( i=0; i < rfield.nflux; ++i )
	{
		double OTS_line_x_opac = rfield.otslinNew[i] * opac.opacity_abs[i];
		double OTS_cont_x_opac = rfield.otsconNew[i] * opac.opacity_abs[i];
		double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] ) ;

		/* get new ots line rates for each cell, but don't let change be too big */
		if( OTS_line_x_opac > old_OTS_line_x_opac[i]*FacBig )
		{
			rfield.otslin[i] = (float)(rfield.otslin[i]*FacBig);
		}
		else if( OTS_line_x_opac < old_OTS_line_x_opac[i]*FacSml )
		{
			rfield.otslin[i] = (float)(rfield.otslin[i]*FacSml);
		}
		else 
		{
			rfield.otslin[i] = (float)(rfield.otslinNew[i]);
		}

		/* get new ots continuum rates for each cell */
		if( OTS_cont_x_opac > old_OTS_cont_x_opac[i]*FacBig )
		{
			rfield.otscon[i] = (float)(rfield.otscon[i]*FacBig);
		}
		else if( OTS_cont_x_opac < old_OTS_cont_x_opac[i]*FacSml )
		{
			rfield.otscon[i] = (float)(rfield.otscon[i]*FacSml);
		}
		else 
		{
			rfield.otscon[i] = (float)(rfield.otsconNew[i]);
		}

		/* zero out the new otscon vector*/
		rfield.otsconNew[i] = 0.;

		/* zero out the new otslin vector*/
		rfield.otslinNew[i] = 0.;

		/* now update to new values */
		old_OTS_line_x_opac[i] = rfield.otslin[i] * opac.opacity_abs[i];
		old_OTS_cont_x_opac[i] = rfield.otscon[i] * opac.opacity_abs[i];

		/* this is local ots continuum created by destroyed diffuse continua, 
		 * currently only two-photon */
		rfield.ConOTS_local_OTS_rate[i] = (float)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);

		/* remember sum of ots rates for convergence criteria */
		*SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];

		rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 
		  rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
		  rfield.ConOTS_local_OTS_rate[i];

		rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
		rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];

	}


	/* >>chng 02 oct 18, add this */
	/* sum of accumulated flux from particular frequency to infinity */
	rfield.flux_accum[rfield.nflux-1] = 0;
	for( i=1; i < rfield.nflux; i++ )
	{
		rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
			rfield.SummedCon[rfield.nflux-i-1];
	}
	/* this has to be positive since is sum of all photons in SummedCon */
	ASSERT( rfield.flux_accum[0]> 0. );

	/* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
	 * between plasma freqnecy and energy where brems is thin */
	ASSERT(rfield.ipPlasma>0 );

	/* >>chng 02 jul 25, set all radiation fields to zero below plasma frequence */
	for( i=0; i < rfield.ipPlasma-1; i++ )
	{
		rfield.otscon[i] = 0.;
		rfield.ConOTS_local_OTS_rate[i] = 0.;
		rfield.ConOTS_local_photons[i] = 0.;
		rfield.otslin[i] = 0.;
		rfield.SummedDif[i] = 0.;
		rfield.OccNumbBremsCont[i] = 0.;
		rfield.SummedCon[i] = 0.;
		rfield.SummedOcc[i] = 0.;
		rfield.ConInterOut[i] = 0.;
	}
	/* this loop evaluates occupation number for brems continuum,
	 * only used for induced two photon emission */
	if( rfield.ipEnergyBremsThin > 0 )
	{
		for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
		{
			/* this corrects for opacity / optical depth in brems - brems opacity goes as
			 * energy squared. */
			/* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
			 * so need this max */
			float factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);

			/* occupation number for black body is 1/ (exp hn/kT) -1) */
			rfield.OccNumbBremsCont[i] = (float)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
		}
	}

#	ifdef DEBUG_FUN
	fputs( " <->RT_OTS_Update()\n", debug_fp );
#	endif
	return;

#	if 0
	OTSOld = 0.;
	OTSNew = 0.;
	BigOTSNew = 0.;
	/* find summed ots rates, old and new */
	for( i=0; i < rfield.nflux; i++ )
	{
		OTSOld += rfield.otscon[i] + rfield.otslin[i];
		OTSNew += rfield.otsconNew[i] + rfield.otslinNew[i];
		if( BigOTSNew < rfield.otsconNew[i] + rfield.otslinNew[i] )
		{
			BigOTSNew = rfield.otsconNew[i] + rfield.otslinNew[i];
		}
	}
	
	/* we now have old and new rates, what is the ratio, and by how much will
	 * we allow this to change? */
	if( BigFrac == 0. || conv.lgSearch || OTSOld < SMALLFLOAT )
	{
		/* this branch, want a clear update of ots rates, either because
		 * requested with BigFrac == 0, or we are in search phase */
		FracOld = 0.;
	}
	else
	{
		if( OTSNew > OTSOld )
		{
			chng = fabs(1. - OTSOld/SDIV(OTSNew) );
		}
		else
		{
			chng = fabs(1. - OTSNew/SDIV(OTSOld) );
		}

		if( chng < BigFrac )
		{
			/* this branch, old and new do not differ by much, so use new */
			FracOld = 0.;
		}
		else
		{
			/* this branch, too large a difference between old and new, cap it to BigFrac */
			FracOld = (1. - BigFrac / chng);
			ASSERT( FracOld >= 0. );
			FracOld = MIN2( 0.25 , FracOld );
		}
	}

	/* fraction old and new ots rates */
	FracNew = 1. - FracOld ;
	fprintf(ioQQQ," DEBUG FracNew\t%.2e\n", FracNew);

	/* remember largest change in ots rates */
	changeOTS = 0.;
	/*fprintf(ioQQQ," sumcontinuum zone%li 1168=%e\n", nzone,rfield.otslin[1167]);*/

	for( i=0; i < rfield.nflux; i++ )
	{
		/* >>chng 01 feb 01, define inverse opacity in safe manner */
		double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] ) ;

		/* remember which energy had largest change in ots rates */
		if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS)
		{
			changeOTS = 
			  fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i]);
		}

		/* >>chng 01 apr 10, only change ots with means if FracOld not zero,
		 * this is to avoid taking means when this routine called by StartEndIter 
		 * to reset sums of rates */
		if( BigFrac > 0. && conv.nTotalIoniz > 0 )
		{
			/* the New vectors are the ones updated by the AddOTS routines,
			 * and have the current OTS rates.  The otscon and otslin vectors
			 * have the rates from the previous refresh of the New vectors*/
			/* here is the refresh.  all were initially set to zero in call to 
			 * RT_OTS_Zero (below) from zero */
			rfield.otscon[i] = (float)((rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew));

			/* this is local ots continuum created by destroyed diffuse continua, 
			 * currently only two-photon */
			/* >>chng 00 dec 27, define this continuum and do not count 2-photon in old var */
			rfield.ConOTS_local_OTS_rate[i] = (float)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);

			/* current ots will be average of old and new */
			rfield.otslin[i] = (float)((rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew));
		}

		/* zero out the new otscon vector*/
		rfield.otsconNew[i] = 0.;

		/* zero out the new otslin vector*/
		rfield.otslinNew[i] = 0.;

		/* remember sum of ots rates for convergence criteria */
		*SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
		{
			/* following should be set true to print strongest ots contributors */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC && (nzone==378)/**/ )
			{
				if( conv.nPres2Ioniz > 3 ) 
					exit(1);
				fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
					conv.nPres2Ioniz,
					rfield.anu[i],
					opac.opacity_abs[i],
					rfield.otscon[i],
					rfield.otslin[i]);
			}
		}

		/* >>chng 00 dec 27, include ConOTS_local_OTS_rate */
		/* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
		rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 
		  rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
		  rfield.ConOTS_local_OTS_rate[i];

		rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
		rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
	}

	ASSERT(*SumOTS >= 0. );

	/* >>chng 02 oct 18, add this */
	/* sum of accumulated flux from particular frequency to infinity */
	rfield.flux_accum[rfield.nflux-1] = 0;
	for( i=1; i < rfield.nflux; i++ )
	{
		rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
			rfield.SummedCon[rfield.nflux-i-1];
	}
	/* this has to be positive since is sum of all photons in SummedCon */
	ASSERT( rfield.flux_accum[0]> 0. );

	/* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
	 * between plasma freqnecy and energy where brems is thin */
	ASSERT(rfield.ipPlasma>0 );

	/* >>chng 02 jul 25, set all radiation fields to zero below plasma frequence */
	for( i=0; i < rfield.ipPlasma-1; i++ )
	{
		rfield.otscon[i] = 0.;
		rfield.ConOTS_local_OTS_rate[i] = 0.;
		rfield.ConOTS_local_photons[i] = 0.;
		rfield.otslin[i] = 0.;
		rfield.SummedDif[i] = 0.;
		rfield.OccNumbBremsCont[i] = 0.;
		rfield.SummedCon[i] = 0.;
		rfield.SummedOcc[i] = 0.;
		rfield.ConInterOut[i] = 0.;
	}
	/* this loop evaluates occupation number for brems continuum,
	 * only used for induced two photon emission */
	if( rfield.ipEnergyBremsThin > 0 )
	{
		for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
		{
			/* this corrects for opacity / optical depth in brems - brems opacity goes as
			 * energy squared. */
			/* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
			 * so need this max */
			float factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);

			/* occupation number for black body is 1/ (exp hn/kT) -1) */
			rfield.OccNumbBremsCont[i] = (float)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
		}
	}

	{
		/* following should be set true to print strongest ots contributors */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>0)/**/ )
		{
			double BigOTS;
			int ipOTS=0;
			BigOTS = -1.;
			/* find the biggest ots contributor */
			/*for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]; i < rfield.nflux; i++ )*/
			for( i=0; i < rfield.nflux; i++ )
			{
				if( (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i] > BigOTS )
				{
					BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
					ipOTS = (int)i;
				}
			}
			fprintf(ioQQQ,
				" sumcontinuunots zone %li SumOTS %.2e %.2eRyd opc:%.2e OTS%.2e lin%.2e con:%.2e %s %s cell%i FracNew %.2f \n",
				nzone,
				*SumOTS,
				rfield.anu[ipOTS] , 
				opac.opacity_abs[ipOTS],
				BigOTS ,
				rfield.otslin[ipOTS],
				rfield.otscon[ipOTS] , 
				/* line	label */
				rfield.chLineLabel[ipOTS] ,
				/* cont label*/
				rfield.chContLabel[ipOTS] ,
				ipOTS,
				FracNew);
		}
	}

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

/* =================================================================== */

/*RT_OTS_Zero zero out some vectors - this is only called when code initialized by ContSetIntensity */
void RT_OTS_Zero( void )
{
	long int i;

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

	/* this loop goes up to nflux itself (<=) since the highest cell
	 * will be used to pass unity through the code to verify integrations */
	for( i=0; i <= rfield.nflux; i++ )
	{
		rfield.SummedDif[i] = 0.;
		/* the four main ots vectors */
		rfield.otscon[i] = 0.;
		rfield.otslin[i] = 0.;
		rfield.otslinNew[i] = 0.;
		rfield.otsconNew[i] = 0.;

		rfield.ConInterOut[i] = 0.;
		rfield.outlin[i] = 0.;
		rfield.outlin_noplot[i] = 0.;
		rfield.SummedDif[i] = 0.;
		/* "zero" for the summed con will be just the incident radiation field */
		rfield.SummedCon[i] = rfield.flux[i];
		rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
		rfield.ConOTS_local_photons[i] = 0.;
		rfield.ConOTS_local_OTS_rate[i] = 0.;
	}

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

/* =================================================================== */

/*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
void RT_OTS_ChkSum(long int ipPnt)
{
	static long int nInsane=0;
	long int i;
	double phisig;
#	define LIM_INSAME_PRT	30

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

	/* this entire sub is a sanity check */
	/* >>chng 02 jul 23, lower bound from 0 to rfield.ipEnergyBremsThin - since now
	 * set radiation field to black body below this energy */
	for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ )
	{
		phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly + 
		  rfield.outlin[i]+
		  rfield.outlin_noplot[i]+
		  rfield.ConOTS_local_OTS_rate[i];
		/* >>chng 02 sep 19, add sec test on SummedDif since it can be zero whild 
		 * phisig is just above small float */
		if( phisig > 0. && rfield.SummedDif[i]> 0.)
		{
			if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 )
			{
				++nInsane;
				/* limit amount of printout - in many failures would print entire
				 * continuum array */
				if( nInsane < LIM_INSAME_PRT )
				{
					fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n", 
					rfield.anu[i], rfield.SummedDif[i]/phisig - 1., i );
					fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n", 
					rfield.SummedDif[i], phisig );
					fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n", 
					rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 
					rfield.outlin[i]+rfield.outlin_noplot[i] );
					fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n", 
					rfield.chLineLabel[i], rfield.chContLabel[i] );
				}

			}
		}
		phisig += rfield.flux[i];
		/* >>chng 02 sep 19, add sec test on SummedDif since it can be zero when 
		 * phisig is just above small float */
		if( phisig > 0. && rfield.SummedDif[i]> 0. )
		{
			if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 )
			{
				++nInsane;
				/* limit amount of printout - in many failures would print entire
				 * continuum array */
				if( nInsane < LIM_INSAME_PRT )
				{
					fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum%3ld, insane SummedCon at energy%.5e error=%.2e i=%ld\n", 
					ipPnt, rfield.anu[i], rfield.SummedCon[i]/phisig - 1., i );
					fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n", 
					rfield.SummedCon[i], phisig );
					fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n", 
					rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 
					rfield.outlin[i]+rfield.outlin_noplot[i], rfield.flux[i] );
					fprintf( ioQQQ, " line continuum here are %s %s\n", 
					rfield.chLineLabel[i], rfield.chContLabel[i]
					);
				}
			}
		}
	}

	if( nInsane > 0 )
	{
		fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n");
		/* TotalInsanity exits after announcing a problem */
		TotalInsanity();
	}

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

/* =================================================================== */

/*RT_OTS_PrtRate print continuum and line ots rates when trace ots is on */
void RT_OTS_PrtRate(
	  /* weakest rate to print */
	  double weak ,
	  /* flag, 'c' continuum, 'l' line, 'b' both */
	  int chFlag )
{
	long int i;

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

	/* arg must be one of these three */
	ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' );

	/* 
	 * both printouts have cell number (on C array scale)
	 * energy in ryd
	 * the actual value of the ots rate 
	 * the ots rate relative to the continuum at that energy 
	 * rate times opacity 
	 * all are only printed if greater than weak 
	 */

	/*===================================================================*/
	/* first print ots continua                                          */
	/*===================================================================*/
	if( chFlag == 'c' || chFlag == 'b' )
	{
		fprintf( ioQQQ, "     DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n",
			weak,fnzone ,TorF(conv.lgConvIoniz) );

		for( i=0; i < rfield.nupper; i++ )
		{
			if( rfield.otscon[i]*opac.opacity_abs[i] > weak )
			{
				fprintf( ioQQQ, "     %4ld%12.4e%12.4e%12.4e%12.4e %s \n", 
				  i, 
				  rfield.anu[i], 
				  rfield.otscon[i], 
				  opac.opacity_abs[i], 
				  rfield.otscon[i]*opac.opacity_abs[i], 
				  rfield.chContLabel[i]);

			}
		}
	}

	/*===================================================================*/
	/* second print ots line rates                                       */
	/*===================================================================*/
	if( chFlag == 'l' || chFlag == 'b' )
	{
		fprintf( ioQQQ, "     DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n",
			weak,fnzone,TorF(conv.lgConvIoniz) );

		for( i=0; i < rfield.nupper; i++ )
		{
			if( rfield.otslin[i]*opac.opacity_abs[i] > weak )
			{
				fprintf( ioQQQ, "     %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n", 
				  i, 
				  rfield.anu[i], 
				  rfield.otslin[i], 
				  opac.opacity_abs[i],
				  rfield.otslin[i]*opac.opacity_abs[i],
				  rfield.chLineLabel[i] ,
				  rfield.line_count[i] );
			}
		}
	}

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