/* 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 */
/*radius_increment do work associated with geometry increments of this zone, called before RT_tau_inc */
/*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "hydrogenic.h"
#include "dynamics.h"
#include "rfield.h"
#include "colden.h"
#include "geometry.h"
#include "opacity.h"
#include "thermal.h"
#include "doppvel.h"
#include "dense.h"
#include "mole.h"
#include "h2.h"
#include "timesc.h"
#include "neutrn.h"
#include "hmi.h"
#include "lines_service.h"
#include "taulines.h"
#include "trace.h"
#include "wind.h"
#include "phycon.h"
#include "pressure.h"
#include "grainvar.h"
#include "molcol.h"
#include "conv.h"
#include "hyperfine.h"
#include "mean.h"
#include "struc.h"
#include "rt.h"
#include "radius.h"
/*cmshft - so compton scattering shift of spectrum */
static void cmshft(void);

#if !defined(NDEBUG)
/*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */
/* function only used in debug mode */
static void pnegopc(void);
#endif

void radius_increment(void)
{
#	if !defined(NDEBUG)
	int lgFlxNeg;
#	endif

	long int nelem,
	  i, 
	  j ,
	  ion,
	  nd,
	  mol;
	double 
	  avWeight,
	  DilutionHere ,
	  escatc, 
	  fmol,
	  opfac, 
	  relec, 
	  rforce, 
	  t;

	double ajmass, 
		Error,
	  rjeans;

	float dradfac;

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

	/* when this sub is called RADIUS is the outer edge of zone */

	if( lgAbort )
	{
#		ifdef DEBUG_FUN
		fputs( " <->radius_increment()\n", debug_fp );
#		endif
		return;
	}

	/* following block only set of asserts to check that iso populations are sane */
#	if !defined(NDEBUG)
	if( !dynamics.lgAdvection )
	{
		long int ipISO ;
		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		{
			for( nelem=ipISO; nelem<LIMELM; ++nelem )
			{
				/* >>chng 05 feb 05, rom div by gas_phase for element, to eden, see explain
				 * for this date below */
				/*if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][nelem-ipISO]/dense.gas_phase[nelem]>1e-10 )*/
				if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10.)
				{
					/* we will check that grnd pops are within this factor of the xIonDense value
					* these are only converged down to some extent - probably this should be related
					* to the ionization convergence error */
					/* >>chng 04 dec 25, move here from conv_base since only valid check when fully converged */
					/* >>chng 05 feb 05, for very highly ionized sims, where N or O is He-like,
					 * the error can approach 1%, and depends on the convergence criteria.
					 * the code does not try to converge these quantities in an absolute sense, 
					 * only the quantites relative to the electron or hydrogen density.  so it is not
					 * safe to do the assert for small values 
					 * change is to only do this branch if abundance is large enough to be detected by
					 * the ionization convergence solvers */
#					define ERR_CHK	1.001
					double err_chk = ERR_CHK;
					/* >>chng 05 sep 02, when low-T solver used soln is approximate,
					 * and it must not matter (lot-T solver should not be used if it
					 * does matter - so use more liberal expectation */
					if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' )
						err_chk *= 10.;
					/* make sure that populations are valid, first check Pop2Ion  */
					if( iso.Pop2Ion[ipISO][nelem][0] >
						dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk )
					{
						fprintf(ioQQQ,
							" DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n", 
							fnzone,
							ipISO , nelem ,
							iso.Pop2Ion[ipISO][nelem][0]*dense.xIonDense[nelem][nelem+1-ipISO] ,
							dense.xIonDense[nelem][nelem-ipISO] ,
							iso.chTypeAtomUsed[ipISO][nelem]);
						fprintf(ioQQQ," level solver found pop_ion_ov_neut of %.5e",
							iso.pop_ion_ov_neut[ipISO][nelem] );
						fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n",
							dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) );
					}

					ASSERT( conv.lgSearch ||
						iso.Pop2Ion[ipISO][nelem][0] <=
						dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk/*(1.+10.*FLT_EPSILON)*/ );
					/* make sure that populations are valid, first check Pop2Ion  */

					if( EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] >
						dense.xIonDense[nelem][nelem-ipISO]*err_chk )
					{
						fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n", 
							fnzone,
							ipISO , nelem ,
							EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO],
							dense.xIonDense[nelem][nelem-ipISO] ,
							iso.chTypeAtomUsed[ipISO][nelem]);
					}
					ASSERT( (conv.lgSearch || !conv.nPres2Ioniz) || 
						EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] <=
						dense.xIonDense[nelem][nelem-ipISO]*err_chk/*(1.+10.*FLT_EPSILON)*/ );
#					undef ERR_CHK
				}
			}
		}
	}
#	endif

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n", 
		  radius.Radius, radius.rinner, radius.drad, radius.drNext, 
		  radius.router[iteration-1], radius.depth );
	}

	/* remember mean and largest errors on electron density */
	Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue);
	if( Error > conv.BigEdenError )
	{
		conv.BigEdenError = (float)Error;
		dense.nzEdenBad = nzone;
	}
	conv.AverEdenError += (float)Error;

	/* remember mean and largest errors between heating and cooling */
	Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot;
	conv.BigHeatCoolError = MAX2((float)Error , conv.BigHeatCoolError );
	conv.AverHeatCoolError += (float)Error;

	/* remember mean and largest pressure errors */
	Error = fabs(pressure.PresTotlCurr - pressure.PresTotlCorrect) / pressure.PresTotlCorrect;
	conv.BigPressError = MAX2((float)Error , conv.BigPressError );
	conv.AverPressError += (float)Error;

	/* integrate total mass over model */
	dense.xMassTotal += dense.xMassDensity * (float)(radius.drad_x_fillfac*radius.r1r0sq);

	/* check cooling time for this zone, remember longest */
	timesc.ttherm = MAX2(timesc.ttherm,1.5*dense.pden*BOLTZMANN*phycon.te/
	  thermal.ctot);

	/* fmol is fraction of hydrogen in form of any molecule or ion */
	/* if N_H_MOLEC != 8, we probably need to change this line... */
	assert (N_H_MOLEC == 8);
	fmol = (hmi.Hmolec[ipMHm] + 2.*(hmi.H2_total + hmi.Hmolec[ipMH2p]))/dense.gas_phase[ipHYDROGEN];

	/* remember the largest fraction that occurs in the model */
	hmi.BiggestH2 = MAX2(hmi.BiggestH2,(float)fmol);

	/* largest fraction of atoms in molecules */
	for( i=0; i<NUM_COMOLE_CALC; ++i )
	{
		if( dense.lgElmtOn[co.nelem_hevmol[i]] )
		{
			float frac = co.hevmol[i] / SDIV(dense.gas_phase[co.nelem_hevmol[i]]);
			/* these are special cases with two heavy elements in one molecule */
			if( co.nelem_hevmol[i]==ipOXYGEN && (i==ipO2 || i==ipO2P || i==ipNO2 || i==ipNO2P) )
				frac *= 2.f;
			else if( co.nelem_hevmol[i]==ipNITROGEN && ( i==ipN2O || i==ipN2 ) )
				frac *= 2.f;
			else if( co.nelem_hevmol[i]==ipSULPHUR && (i==ipS2||i==ipS2P ) )
				frac *= 2.f;			
			else if( co.nelem_hevmol[i]==ipCARBON && (i==ipC2||i==ipC2P ) )
				frac *= 2.f;
			co.xMoleFracMax[i] = MAX2(frac,co.xMoleFracMax[i]);
		}
	}

	/* H 21 cm equilibrium timescale, H21cm returns H (not e) collisional deexcit rate (not cs) */
	t = H21cm_H_atom( phycon.te )* dense.xIonDense[ipHELIUM][ipHYDROGEN] +
		/* >>chng 02 feb 14, add electron term as per discussion in */
		/* >>refer	H1	21cm	Liszt, H., 2001, A&A, 371, 698 */
		H21cm_electron( phycon.te )*dense.eden;

	if( t > SMALLFLOAT )
	{
		timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm );
	}

	/* remember longest CO timescale */
	if( dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0] > SMALLFLOAT )
	{
		double timemin;
		double product = SDIV(dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0]);
		timemin = MIN2(timesc.AgeCOMoleDest ,
			timesc.AgeCOMoleDest* co.hevmol[ipCO]/ product );
		/* this is rate CO is destroyed, equal to formation rate in equilibrium */
		timesc.BigCOMoleForm = MAX2(  timesc.BigCOMoleForm , timemin );
	}

	/* remember longest H2  destruction timescale timescale */
	timesc.time_H2_Dest_longest = MAX2(timesc.time_H2_Dest_longest, timesc.time_H2_Dest_here );

	/* remember longest H2 formation timescale timescale */
	timesc.time_H2_Form_longest = MAX2( timesc.time_H2_Form_longest , timesc.time_H2_Form_here );

	/* increment counter if this zone possibly thermally unstable
	 * this flag was set in ionte, deriv of heating and cooling negative */
	if( thermal.lgUnstable )
	{
		thermal.nUnstable += 1;
	}

	/* remember Stromgren radius - where hydrogen ionization falls below half */
	if( !rfield.lgUSphON && dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN] > 0.49 )
	{
		rfield.rstrom = (float)radius.Radius;
		rfield.lgUSphON = TRUE;
	}

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

	relec = 0.;
	rforce = 0.;

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

	/* diffuse continua factor
	 * if lgSphere then all diffuse radiation is outward (COVRT=1)
	 * lgSphere not set then COVRT=0.0 */

	/* begin sanity check */
	/* only do this test in debug mode */
#	if !defined(NDEBUG)
	lgFlxNeg = FALSE;
	for( i=0; i < rfield.nflux; i++ )
	{
		if( rfield.flux[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[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[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[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 punch 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();
		puts( "[Stop in radius_increment]" );
		cdEXIT(EXIT_FAILURE);
	}
	/*end sanity check*/
#	endif

	/* rfield.lgOpacityFine flag set false with no fine opacities command */
	/* only update this if results will be punched */
	/* tests show that always evaluating this changes fast run of
	 * parispn from 26.7 sec to 35.1 sec */
	/* >>chng 05 mar 01, must always update them since used for transmission */
	if( rfield.lgOpacityFine /*&& rfield.lgPunchOpacityFine*/ )
	{
		dradfac = (float)radius.drad_x_fillfac;
		/* increment the fine opacity array */
		for( i=0; i<rfield.nfine; ++i )
		{
			rfield.fine_opt_depth[i] += 
				rfield.fine_opac[i]*dradfac;
		}

		/* >>chng 05 feb 28, evaluate fine opacity transmission coefficient for transmission
		 * through this zone. 
		 * this is assumed to be a scattering opacity, only do this if including scattering */
		if( opac.lgScatON )
		{
			/* sum over coarse continuum */
			for( i=0; i < rfield.nflux-1; i++ )
			{
				/* find transmission coef if lower and upper bounds of coarse continuum is within
				 * boundaries of fine continuum 
				 * unity is default */
				if( rfield.ipnt_coarse_2_fine[i] && rfield.ipnt_coarse_2_fine[i+1] )
				{
					rfield.trans_coef[i] = 0.;
					for( j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j )
					{
						/* find transmission coef for this zone */
						/* here assume continuum transmission is exp(-tau) */
						/*rfield.trans_coef[i] += 
							(float)sexp( rfield.fine_opt_depth[j] );*/
						rfield.trans_coef[i] += 1.f / (1.f + rfield.fine_opac[j]*dradfac );
					}
					/* first branch is normal case, where fine continuum is finer than
					 * coarse continuum.  But, when end temp is very high, fine continuum is
					 * very coarse, so may be just one cell, and following will not pass */
					if( rfield.ipnt_coarse_2_fine[i+1]>rfield.ipnt_coarse_2_fine[i] )
					{
						rfield.trans_coef[i] /= (rfield.ipnt_coarse_2_fine[i+1]-rfield.ipnt_coarse_2_fine[i]);
					}
					else
					{
						/* in case where fine is coarser than coarse, just use firs cell */
						rfield.trans_coef[i] = 
							1.f / (1.f + rfield.fine_opac[rfield.ipnt_coarse_2_fine[i]]*dradfac );
					}
				}
				else
				{
					/* this branch, not within find continuum energy range - transmission unity */
					rfield.trans_coef[i] = 1.f;
				}
			}
		}
	}

	/* attenuate the flux array after doing radiative acceleration */
	escatc = 6.65e-25*dense.eden;

	/* this loop should not be to <= nflux since we not want to clobber the
	 * continuum unit integration */
	for( i=0; i < rfield.nflux; i++ )
	{
		/* sum total continuous optical depths */
		opac.TauAbsGeo[0][i] += (float)(opac.opacity_abs[i]*radius.drad_x_fillfac);
		opac.TauScatGeo[0][i] += (float)(opac.opacity_sct[i]*radius.drad_x_fillfac);

		/* following only optical depth to illuminated face */
		opac.TauAbsFace[i] += (float)(opac.opacity_abs[i]*radius.drad_x_fillfac);
		opac.TauScatFace[i] += (float)(opac.opacity_sct[i]*radius.drad_x_fillfac);

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

		/* TMN is array of scale factors which account for attenuation
		 * of continuum across the zone (necessary to conserve energy
		 * at the 1% - 5% level.) sphere effects in
		 * drNext was set by NEXTDR and will be next dr */
		/* compute both total and thompson scat rad acceleration */
		rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
			rfield.anu[i]*(opac.opacity_abs[i] + 
		  opac.opacity_sct[i]);

		relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
			escatc)*rfield.anu[i];

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

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

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

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

		rfield.flux[i] *= (float)opfac;
		/* >>chng 03 nov 08, update SummedCon here since flux changed */
		/* >>chng 05 mar 5, add transmission coefficient due to fine opacites */
		/*rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];*/
		rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];

		/* outward lines and continua */
		rfield.ConInterOut[i] *= (float)opfac;
		rfield.outlin[i] *= (float)opfac;
		rfield.outlin_noplot[i] *= (float)opfac;

		/* this simulates the behavior in C90 and before,
		 * first attenuate flux acculumated so far, then
		 * add on flux from this zone 
		rfield.ConEmitOut[i] *= (float)opfac;*/

		rfield.ConEmitOut[i] *= (float)opfac;
		rfield.ConEmitOut[i] += rfield.ConEmitLocal[i]*(float)radius.dVolOutwrd*opac.tmn[i]/**/;

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

		/* >>chng 00 oct 03, add diffuse continua */
		/* local diffuse continua */
		rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[i]*rfield.convoc[i];

		/* >>chng 05 feb 16, occupation number of outward diffuse continuum */
		rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[i]*rfield.convoc[i];

	}

	/* begin 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 */
	/* >>chng 02 dec 05, add test for small float, protect against models where we have
	 * gone below smallfloat, and so float is ragged */
	if( rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
		(rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/SDIV(rfield.FluxSave[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-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 tways - 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[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
			(double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
			(double)rfield.FluxSave[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);

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

		/* 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 ) / t > 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.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-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();
		}
	}
	/* end sanity check */

	/* do scattering opacity, intensity goes as 1/(1+tau) 
	 * scattering opacity is turned off when sphere is set */
	if( opac.lgScatON )
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			/* assume half forward scattering
			 * opfac = 1. / ( 1. + dReff*0.5 * scatop(i) )
			 * >>chng 97 apr 25, remove .5 to get agreement with
			 * Lightman and White equation 11 in small epsilon limit,
			 * >>refer	continuum	RT	Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */
			opfac = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
			ASSERT( opfac <= 1.0 );
			/* >>chng 05 mar 03, include transmission coefficient due to fine opacities 
			opfac *= rfield.trans_coef[i];*/
			rfield.flux[i] *= (float)opfac;
			rfield.ConInterOut[i] *= (float)opfac;
			rfield.ConEmitOut[i] *= (float)opfac;
			rfield.outlin[i] *= (float)opfac;
			rfield.outlin_noplot[i] *= (float)opfac;
		}
	}

	/* now do slight reshuffling of energy due to compton scattering */
	cmshft();

	relec *= EN1RYD;

	/* radiative acceleration; xMassDensity is gm per cc, eval when PTOT called */
	t = 1./SPEEDLIGHT/dense.xMassDensity;
	/* >>chng 03 feb 01, lgLineRadPresOn was missing */
	wind.AccelLine = (float)(RT_line_driving()*t)*pressure.lgLineRadPresOn;
	/* >>chng 01 aug 03, add lgConPress to kill AccelCont here */
	wind.AccelCont = (float)(rforce*EN1RYD*t)*pressure.lgContRadPresOn;
	/* this is numerically unstable */
	wind.AccelPres = 0.;
	/* total acceleration */
	wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;
	/* remember the largest value */
	wind.AccelMax = (float)MAX2(wind.AccelMax,wind.AccelTot);
	/* fprintf(ioQQQ," winddebug  %.2f %e  %e %e %e\n", 
		fnzone, wind.AccelTot , wind.AccelMax,	wind.AccelLine ,
		wind.AccelCont); */

	/* force multiplier; RELEC*T can be zero for very low densities */
	if( relec*t > SMALLFLOAT )
	{
		wind.fmul = (float)((wind.AccelLine + wind.AccelCont)/(relec*t));
	}
	else
	{
		wind.fmul = 0.;
	}

	/* following is integral of radiative force */
	pressure.pinzon = (float)(wind.AccelTot*dense.xMassDensity*radius.drad_x_fillfac*geometry.AngleIllum);
	/*fprintf(ioQQQ," debuggg pinzon %.2f %.2e %.2e %.2e\n", 
		fnzone,pressure.pinzon,dense.xMassDensity,wind.AccelTot);*/
	pressure.PresInteg += pressure.pinzon;

	/* sound is sound travel time, sqrt term is sound speed */
	timesc.sound_speed_isothermal = sqrt(pressure.PresGasCurr/dense.xMassDensity);
	/* adiabatic sound speed assuming monatomic gas - gamnma is 5/3*/
	timesc.sound_speed_adiabatic = sqrt(5.*pressure.PresGasCurr/(3.*dense.xMassDensity) );
	timesc.sound += radius.drad_x_fillfac / timesc.sound_speed_isothermal;

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

	/* following radiation factors are extinguished by 1/r**2 and e- opacity
	 * also bound electrons */
	
	/* do all emergent spectrum from illuminated face if model is NOT spherical */
	if( !geometry.lgSphere )
	{
		double Reflec_Diffuse_Cont;

		/* >>chng 01 jul 14, from lower limit of 0 to plasma frequency -
		 * never should have added diffuse emission from below the plasma frequency */
		for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
		{
			if( opac.TauAbsGeo[0][i] < 30. )
			{
				/* ConEmitLocal is diffuse emission per unit vol, fill fac
				 * the 1/2 comes from isotropic emission
				 * scatsv(i) = (flux(i)+ConInterOut(i))*0.5*scatop(i)
				 * >>chng 97 apr 25, remove fac of 2 to
				 * get agreement wth lightman&white
				 * Reflec_Diffuse_Cont = (ConEmitLocal(i)+scatsv(i)) * dReff *e2TauAbs(i)
				 * >>chng 97 may 08, needed 1/2 in front of ConEmitLocal 
				 * had been there, but also with scatsv, and lost both
				 * when /2 removed to get agreement with lightman and white */
				/* >>chng 99 dec 21, add ConLocNoInter to account for rec continua
				 * >>chng 01 jul 01, remove ConLocNoInter from code - only had thick brems */
				Reflec_Diffuse_Cont = rfield.ConEmitLocal[i]/2.*
					radius.drad_x_fillfac * opac.e2TauAbs[i]*radius.r1r0sq;

				/* Reflec_Diffuse_Cont = (ConEmitLocal(i)+scatsv(i)) * dReff/2. *e2TauAbs(i)
				 * ConEmitReflec - reflected diffuse continuum */
				rfield.ConEmitReflec[i] += (float)(Reflec_Diffuse_Cont);

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

				/* >>chng 97 apr 25, remove fac of 2 to get
				 * agreement wth lightman&white */
				rfield.reflin[i] += (float)(rfield.outlin[i]*opac.opacity_sct[i]*
				  radius.drad_x_fillfac*opac.e2TauAbs[i]*radius.r1r0sq);
			}
		}
	}

	/* following is general method to find means weighted by various functions
	 * called in IterStart to initialize to zero, called here to put in numbers
	 * results will be weighted by radius and volume
	 * this is the only place things must be entered to create averages 
	 * code is in mean.c */
	aver("zone",1.,1.,"          ");
	aver("doit",phycon.te,1.,"    Te    ");
	aver("doit",phycon.te,dense.eden,"  Te(Ne)  ");
	aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHYDROGEN][1]," Te(NeNp) ");
	aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][1]  ," Te(NeHe+)");
	aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][2]  ,"Te(NeHe2+)");
	aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][1]  ," Te(NeO+) " );
	aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][2]  ," Te(NeO2+)");
	/*aver("doit",phycon.te,hmi.Hmolec[ipMH2g],"  Te(H2)  ");*/
	aver("doit",phycon.te,hmi.H2_total,"  Te(H2)  ");
	aver("doit",dense.gas_phase[ipHYDROGEN],1.,"   N(H)   ");
	aver("doit",dense.eden,dense.xIonDense[ipOXYGEN][2],"  Ne(O2+) ");
	aver("doit",dense.eden,dense.xIonDense[ipHYDROGEN][1],"  Ne(Np)  ");

	/* save information about structure of model, now used to get t^2 */
	/*nd = MIN2(nzone , NZLIM)-1;*/
	/* max because if program aborts during search phase, will get to here
	 * with nzone = -1 */
	nd = MAX2( 0, nzone-1 );
	/* this is number of struc.xx zones with valid data */
	struc.nzone = nd+1;
	ASSERT(nd>=0 && nd < struc.nzlim );
	struc.testr[nd] = (float)phycon.te;
	/* number of particles per unit vol */
	struc.DenParticles[nd] = dense.pden;
	/* >>chng 02 May 2001 rjrw: add hden for dilution */
	struc.hden[nd] = (float)dense.gas_phase[ipHYDROGEN];
	/* total grams per unit vol */
	struc.DenMass[nd] = dense.xMassDensity;
	struc.heatstr[nd] = thermal.htot;
	struc.coolstr[nd] = thermal.ctot;
	struc.volstr[nd] = (float)radius.dVeff;
	struc.drad[nd] = (float)radius.drad;
	struc.drad_x_fillfac[nd] = (float)radius.drad_x_fillfac;
	struc.histr[nd] = dense.xIonDense[ipHYDROGEN][0];
	struc.hiistr[nd] = dense.xIonDense[ipHYDROGEN][1];
	struc.ednstr[nd] = (float)dense.eden;
	struc.o3str[nd] = dense.xIonDense[ipOXYGEN][2];
	struc.pressure[nd] = (float)pressure.PresTotlCurr;
	struc.PresRadCurr[nd] = (float)pressure.PresRadCurr;
	struc.GasPressure[nd] = (float)pressure.PresGasCurr;
	struc.depth[nd] = (float)radius.depth;
	/* save absorption optical depth from illuminated face to current position */
	struc.xLyman_depth[nd] = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]];
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
	{
		struc.gas_phase[nelem][nd] = dense.gas_phase[nelem];
		for( ion=0; ion<nelem+2 ; ++ion )
		{
			struc.xIonDense[nelem][ion][nd] = dense.xIonDense[nelem][ion];
		}
	}

	/* the hydrogen molecules */
	for (mol=0;mol<N_H_MOLEC;mol++) 
	{
		struc.Hmolec[mol][nd] = hmi.Hmolec[mol];
	}

	/* the heavy element molecules */
	for (mol=0;mol<NUM_COMOLE_CALC;mol++) 
	{
		struc.COmolec[mol][nd] = co.hevmol[mol];
	}

	colden.dlnenp += dense.eden*(double)(dense.xIonDense[ipHYDROGEN][1])*radius.drad_x_fillfac;
	colden.dlnenHep += dense.eden*(double)(dense.xIonDense[ipHELIUM][1])*radius.drad_x_fillfac;
	colden.dlnenHepp += dense.eden*(double)(dense.xIonDense[ipHELIUM][2])*radius.drad_x_fillfac;

	/* >>chng 05 jan 09, add integral of n(H0) / Tspin */
	colden.H0_ov_Tspin += (double)(dense.xIonDense[ipHYDROGEN][0]) / hyperfine.Tspin21cm*radius.drad_x_fillfac;

	/* >>chng 05 Mar 07, add integral of n(OH) / Tspin */
	colden.OH_ov_Tspin += (double)(co.hevmol[ipOH]) / phycon.te*radius.drad_x_fillfac;

	/* this is Lya excitation temperature, moved here from pressuretotal */
	hydro.TexcLya = (float)TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );

	/* count number of times Lya excit temp hotter than gas */
	if( hydro.TexcLya > phycon.te )
	{
		hydro.nLyaHot += 1;
		if( hydro.TexcLya > hydro.TLyaMax )
		{
			hydro.TLyaMax = hydro.TexcLya;
			hydro.TeLyaMax = (float)phycon.te;
			hydro.nZTLaMax = nzone;
		}
	}

	/* column densities in various species */
	colden.colden[ipCOL_HTOT] += (float)(dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac);
	ASSERT( colden.colden[ipCOL_HTOT] >SMALLFLOAT );
	colden.colden[ipCOL_HMIN] += hmi.Hmolec[ipMHm]*(float)radius.drad_x_fillfac;
	/* >>chng 02 sep 20, from htwo to H2_total */
	/* >>chng 05 mar 14, rather than H2_total, give H2g and H2s */
	colden.colden[ipCOL_H2g] += hmi.Hmolec[ipMH2g]*(float)radius.drad_x_fillfac;
	colden.colden[ipCOL_H2s] += hmi.Hmolec[ipMH2s]*(float)radius.drad_x_fillfac;
	/* this is a special form of column density - should be proportional to total shielding */
	colden.coldenH2_ov_vel += hmi.H2_total*(float)radius.drad_x_fillfac / DoppVel.doppler[LIMELM+2];
	colden.colden[ipCOL_HeHp] += hmi.Hmolec[ipMHeHp]*(float)radius.drad_x_fillfac;
	colden.colden[ipCOL_H2p] += hmi.Hmolec[ipMH2p]*(float)radius.drad_x_fillfac;
	colden.colden[ipCOL_H3p] += hmi.Hmolec[ipMH3p]*(float)radius.drad_x_fillfac;
	colden.colden[ipCOL_Hp] += dense.xIonDense[ipHYDROGEN][1]*(float)radius.drad_x_fillfac;
	colden.colden[ipCOL_H0] += dense.xIonDense[ipHYDROGEN][0]*(float)radius.drad_x_fillfac;

	colden.colden[ipCOL_elec] += (float)(dense.eden*radius.drad_x_fillfac);
	/* He I t2S column density*/
	colden.He123S += dense.xIonDense[ipHELIUM][1]*
		iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe2s3S]*(float)radius.drad_x_fillfac;

	/* the ortho and para column densities */
	h2.ortho_colden += h2.ortho_density*radius.drad_x_fillfac;
	h2.para_colden += h2.para_density*radius.drad_x_fillfac;
	ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / SDIV( hmi.H2_total ) < 1e-4 );
	/*fprintf(ioQQQ,"DEBUG ortho para\t%.3e\t%.3e\ttot\t%.3e\t or pa colden\t%.3e\t%.3e\n", 
		h2.ortho_density, h2.para_density,hmi.H2_total,
		h2.ortho_colden , h2.para_colden);*/

	/* the CII and SiII atoms are resolved */
	for( i=0; i < 5; i++ )
	{
		/* pops and column density for C II atom */
		colden.C2Colden[i] += colden.C2Pops[i]*(float)radius.drad_x_fillfac;
		/* pops and column density for SiII atom */
		colden.Si2Colden[i] += colden.Si2Pops[i]*(float)radius.drad_x_fillfac;
	}
	for( i=0; i < 3; i++ )
	{
		/* pops and column density for CI atom */
		colden.C1Colden[i] += colden.C1Pops[i]*(float)radius.drad_x_fillfac;
		/* pops and column density for OI atom */
		colden.O1Colden[i] += colden.O1Pops[i]*(float)radius.drad_x_fillfac;
	}
	for( i=0; i < 4; i++ )
	{
		/* pops and column density for CIII atom */
		colden.C3Colden[i] += colden.C3Pops[i]*(float)radius.drad_x_fillfac;
	}

	/* now add total molecular column densities */
	molcol("ADD ",ioQQQ);

	/* increment forming the mean ionizatin and temperature */
	MeanInc();

	/*-----------------------------------------------------------------------*/

	/* calculate average atomic weight per hydrogen of the plasma */
	avWeight = 0.;
	for( nelem=0; nelem < LIMELM; nelem++ ) 
	{
		avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
	}
	avWeight /= dense.gas_phase[ipHYDROGEN];

	/* compute some average grain properties */
	rfield.opac_mag_B_point = 0.;
	rfield.opac_mag_V_point = 0.;
	rfield.opac_mag_B_extended = 0.;
	rfield.opac_mag_V_extended = 0.;
	if( gv.lgGrainPhysicsOn )
	{
		for( nd=0; nd < gv.nBin; nd++ )
		{
			gv.bin[nd]->avdust += gv.bin[nd]->tedust*(float)radius.drad_x_fillfac;
			gv.bin[nd]->avdft += gv.bin[nd]->DustDftVel*(float)radius.drad_x_fillfac;
			gv.bin[nd]->avdpot += (float)(gv.bin[nd]->dstpot*EVRYD*radius.drad_x_fillfac);
			gv.bin[nd]->avDGRatio += (float)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]*
				gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight*
				radius.drad_x_fillfac);

			/* this is total extinction in magnitudes at V and B, for a point source 
			 * total absorption and scattering,
			 * does not discount forward scattering to be similar to stellar extinction
			 * measurements made within ism */
			rfield.extin_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
				radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
			rfield.opac_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
				dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

			rfield.extin_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
				radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
			rfield.opac_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
				dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

			/* this is total extinction in magnitudes at V and B, for an extended source 
			 * total absorption and scattering,
			 * DOES discount forward scattering to apply for extended source like Orion */
			rfield.extin_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
				radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
			rfield.opac_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
				dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;

			rfield.extin_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
				radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
			rfield.opac_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
				gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
				dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
		}
	}

	/* there are some quantities needed to calculation the Jeans mass and radius */
	colden.TotMassColl += dense.xMassDensity*(float)radius.drad_x_fillfac;
	colden.tmas += (float)phycon.te*dense.xMassDensity*(float)radius.drad_x_fillfac;
	colden.wmas += dense.wmole*dense.xMassDensity*(float)radius.drad_x_fillfac;

	/* now find minimum Jeans length and mass; length in cm */
	rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity*
	  geometry.FillFac))/2.;

	/* minimum Jeans mass in gm */
	ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity*
	  geometry.FillFac);

	/* now remember smallest */
	colden.rjnmin = MIN2(colden.rjnmin,(float)rjeans);
	colden.ajmmin = MIN2(colden.ajmmin,(float)ajmass);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " radius_increment returns\n" );
	}

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

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

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

	if( opac.lgNegOpacIO )
	{
		/* option to punch negative opacities */
		if( NULL==(ioFile=fopen("negopc.txt","w")) )
		{
			fprintf( ioQQQ,"pnegopc could not open negopc.txt for writing\n");
			puts( "[Stop in opac0]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( i=0; i < rfield.nflux; i++ )
		{
			fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i], 
			  opac.opacity_abs[i] );
		}
		fclose( ioFile);
	}


#	ifdef DEBUG_FUN
	fputs( " <->pnegopc()\n", debug_fp );
#	endif
	return;
}
#endif
/* Note - when this file is compiled in fast optimized mode,
 * a good compiler will complain that the routine pnegopc is a
 * static, local routine, but that it has not been used.  It is
 * indeed used , but only when NDEBUG is not set, so it is only
 * used when the code is compiled in debug mode.  So this is
 * not a problem. */


/*cmshft - so compton scattering shift of spectrum */
static void cmshft(void)
{
	long int i;

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

	/* first check whether compton scattering is in as heat/cool */
	if( rfield.comoff == 0. )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->cmshft()\n", debug_fp );
#		endif
		return;
	}

	if( rfield.comoff != 0. )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->cmshft()\n", debug_fp );
#		endif
		return;
	}

	/* do reshuffle */
	for( i=0; i < rfield.nflux; i++ )
	{
		continue;
		/* watch this space for some really great code!!!!
		 * COMUP needs factor of TE to be compton cooling */
	}

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

