/*RadInc do work associated with geometry increments of this zone, called before TauInc */
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "opacpoint.h"
#include "hydrogenic.h"
#include "rfield.h"
#include "showme.h"
#include "colden.h"
#include "exptau.h"
#include "e2tau.h"
#include "linelabl.h"
#include "opac.h"
#include "deptau.h"
#include "timesc.h"
#include "csphot.h"
#include "flxsav.h"
#include "crsnut.h"
#include "hlmax.h"
#include "hcbr.h"
#include "uspher.h"
#include "neutrn.h"
#include "ionfracs.h"
#include "hevmolec.h"
#include "hmi.h"
#include "trace.h"
#include "sphere.h"
#include "wind.h"
#include "zonecnt.h"
#include "h.h"
#include "phycon.h"
#include "radius.h"
#include "densty.h"
#include "nh.h"
#include "facexp.h"
#include "itercnt.h"
#include "pressure.h"
#include "angllum.h"
#include "tcool.h"
#include "nustbl.h"
#include "hphoto.h"
#include "insane.h"
#include "e2.h"
#include "sexp.h"
#include "pnegopc.h"
#include "forlin.h"
#include "cmshft.h"
#include "radinc.h"

void RadInc(void)
{
#	if !defined(NDEBUG)
	int lgFlxNeg;
#endif
	long int i, 
	  limit, 
	  n;
	double crs, 
	  DilutionHere ,
	  escatc, 
	  flx, 
	  hlrate, 
	  opfac, 
	  relec, 
	  rforce, 
	  scatsv[NCELL], 
	  t;

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

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" RadInc 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[IterCnt.iter-1], radius.depth );
	}

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

	/* remember longest CO timescale */
	if( xIonFracs[1][5]*xIonFracs[1][7] > SMALLFLOAT )
	{
		double timemin;
		timemin = MIN2(timesc.AgeCOMoleDest ,
			timesc.AgeCOMoleDest* hevmolec.hevmol[IPCO-1]/
			MAX2(SMALLFLOAT,xIonFracs[1][5]*xIonFracs[1][7]) );
		/* this is rate CO is destroyed, equal to formation rate in equilibrium */
		timesc.BigCOMoleForm = MAX2(  timesc.BigCOMoleForm , timemin );
			/*timesc.AgeCOMoleDest* hevmolec.hevmol[IPCO-1]/
			(xIonFracs[1][5]*xIonFracs[1][2]) );*/
	}

	/* remember longest H2 timescale */
	if( (double)xIonFracs[1][0]*(double)xIonFracs[1][0] > SMALLFLOAT )
	{
		double timemin;
		timemin = MIN2(timesc.AgeH2MoleDest ,
			timesc.AgeH2MoleDest*hmi.htwo /((double)xIonFracs[1][0]*(double)xIonFracs[1][0]));
		timesc.BigH2MoleForm = MAX2(timesc.BigH2MoleForm, timemin );
			/*timesc.AgeH2MoleDest*hmi.htwo /(xIonFracs[1][0]*xIonFracs[1][0])));*/
	}

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

	/* find Stromgren radius */
	if( !uspher.lgUSphON && h.hi/phycon.hden > 0.49 )
	{
		uspher.rstrom = (float)radius.Radius;
		uspher.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.outcon[i], 
			  rfield.otscon[i] + rfield.otslin[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, " RadInc 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, " RadInc 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, " RadInc finds negative tmn.\n" );
			fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
			  opac.tmn[i], rfield.anu[i], i, LineLabl.chLineLabel[i] );
			lgFlxNeg = TRUE;
		}
		if( rfield.otslin[i] < 0. )
		{
			fprintf( ioQQQ, " RadInc 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, LineLabl.chLineLabel[i]  );
			lgFlxNeg = TRUE;
		}
		if( rfield.outlin[i] < 0. )
		{
			fprintf( ioQQQ, " RadInc 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, LineLabl.chLineLabel[i]  );
			lgFlxNeg = TRUE;
		}
		if( rfield.outcon[i] < 0. )
		{
			fprintf( ioQQQ, " RadInc finds negative intensity in outcon.\n" );
			fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
			  rfield.outcon[i], rfield.anu[i], i, LineLabl.chContLabel[i]  );
			lgFlxNeg = TRUE;
		}
		if( opac.opac[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", 
		  ZoneCnt.nzone );
		ShowMe();
		puts( "[Stop in RadInc]" );
		exit(1);
	}
	/*end sanity check*/
#	endif

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

	for( i=0; i < rfield.nflux; i++ )
	{
		/* sum total continuous optical depths */
		opac.tauabs[0][i] += (float)(opac.opac[i]*radius.dReff);

		/* following only depth to illuminated face */
		deptau.depabs[i] += (float)(opac.opac[i]*radius.dReff);
		deptau.depsct[i] += (float)(opac.scatop[i]*radius.dReff);

		/* these are total in inward direction, large if sherical */
		opac.tausct[0][i] += (float)(opac.scatop[i]*radius.dReff);
		opac.TauTotal[0][i] = opac.tauabs[0][i] + opac.tausct[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.outcon[i] + 
		  rfield.ConOutNoInter[i])*rfield.anu[i]*(opac.opac[i] + 
		  opac.scatop[i]);

		relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outcon[i] + 
		  rfield.ConOutNoInter[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 */
		facexp.ExpZone[i] = sexp(opac.opac[i]*radius.dReff*Angllum.AngleIllum);

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

		/* e2(tau) in inward direction, up to illuminated face */
		e2tau.e2TauAbs[i] = (float)e2(deptau.depabs[i],exptau.ExpmTau[i]);
		assert( e2tau.e2TauAbs[i] <= 1. );

		/* DilutionHere is square of ratio of inner to outer radius */
		opfac = facexp.ExpZone[i]*DilutionHere;

		rfield.flux[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.ConOutNoInter[i] *= (float)opfac;
		rfield.ConOutNoInter[i] += rfield.ThrowOutNoInter[i]*(float)radius.dVolOutwrd*opac.tmn[i]/**/;

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

		/* set occupation numbers */
		rfield.occnum[i] = rfield.flux[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 */
	if( (rfield.flux[nh.ipHn[0][0]-1]/MAX2(SMALLFLOAT,flxsav.FluxSave[nh.ipHn[0][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 */
		flx = -log((double)rfield.flux[nh.ipHn[0][0]-1]/
			(double)opac.tmn[nh.ipHn[0][0]-1]/
			(double)flxsav.FluxSave[nh.ipHn[0][0]-1]);

		/* this is computed absorption optical depth to illuminated face */
		t = deptau.depabs[nh.ipHn[0][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( flx - t ) / t > 0.01 && 
			/* for very small inner optical depths, the tmn correction is major,
			 * and this test is not precise */
			fabs(flx-t)>MAX2(0.001,1.-opac.tmn[nh.ipHn[0][0]-1]) )
		{
			/* in print below must add extra HI col den since this will later be
			 * incremented in TauInc */
			fprintf( ioQQQ,
				" RadInc lyman continuum insanity zone %li, attenuat, atomic taus=%g %g H=%g\n",
				ZoneCnt.nzone, flx , t ,
				6.327e-18*(h.hi*radius.dReff+coldenCom.colden[IPCHI-1]) );
			insane();
		}
	}
	/* end sanity check */

	/* do scattering opacity, intensity goes as 1/(1+tau) */
	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 formula */
			opfac = 1./(1. + radius.dReff*opac.scatop[i]);
			rfield.flux[i] *= (float)opfac;
			rfield.ConOutNoInter[i] *= (float)opfac;
			rfield.outcon[i] *= (float)opfac;
			rfield.outlin[i] *= (float)opfac;
			/* dReff for following appears when it is used
			 * >>chng 96 aug 21, extra factor of two in below */
			scatsv[i] = (rfield.flux[i] + rfield.outcon[i])*opac.scatop[i];
		}
	}
	else
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			scatsv[i] = 0.;
		}
	}

	/* 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/densty.xMassDensity;
	wind.AccelLine = (float)(forlin()*t);
	wind.AccelCont = (float)(rforce*EN1RYD*t);
	wind.AccelTot = wind.AccelCont + wind.AccelLine;
	wind.AccelMax = (float)MAX2(wind.AccelMax,wind.AccelTot);
	/* this shows that AccelLine becomes unphysicsl, presumaby due to forlin
	 * fprintf(ioQQQ," wind te %e  %e %e %e\n", phycon.te,	wind.AccelLine ,
		wind.AccelCont ,wind.AccelTot );*/

	/* 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*densty.xMassDensity*radius.dReff*Angllum.AngleIllum);
	pressure.PresInteg += pressure.pinzon;

	/* sound is sound travel time, sqrt term is sound speed */
	timesc.sound += radius.dReff/sqrt(pressure.GasPres/densty.xMassDensity);

	/* attenuate neutrons if they are present */
	if( neutrn.lgNeutrnHeatOn )
	{
		/* correct for optical depth effects */
		neutrn.totneu *= (float)sexp(crsnut.CrsSecNeutron*phycon.hden*radius.dReff*Angllum.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 continua if model is NOT spherical */
	if( !sphere.lgSphere )
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			if( opac.tauabs[0][i] < 30. )
			{
				/* DiffLocal is diffuse emission per unit vol, fill fac
				 * the 1/2 comes from isotropic emission
				 * scatsv(i) = (flux(i)+outcon(i))*0.5*scatop(i)
				 * >>chng 97 apr 25, remove fac of 2 to
				 * get agreement wth lightman&white
				 * flx = (DiffLocal(i)+scatsv(i)) * dReff *e2TauAbs(i)
				 * >>chng 97 may 8, needed 1/2 in front of DiffLocal 
				 * 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 ThrowOutNoInter to account for rec continua*/
				flx = (rfield.DiffLocal[i]/2. + scatsv[i] +rfield.ThrowOutNoInter[i]/2.)*
					radius.dReff * e2tau.e2TauAbs[i];

				/* flx = (DiffLocal(i)+scatsv(i)) * dReff/2. *e2TauAbs(i)
				 * refcon - reflected diffuse plus back-scattered incident cont */
				rfield.refcon[i] += (float)(flx*radius.r1r0sq);

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

				rfield.reflux[i] += (float)(rfield.flux[i]*opac.scatop[i]*
				  radius.dReff*e2tau.e2TauAbs[i]*radius.r1r0sq);
			}
		}
	}

	/* remember largest line photoionization rate for COMNT */
	limit = MIN2(nhlevel-3,NHPLPHOT-2);
	for( n=3; n <= limit; n++ )
	{
		crs = csphot( HydroLines[0][n][n-1].ipCont ,
		  nh.ipHn[0][n+1] , OpacPoint.ipHOpac[0][n+1] );

		hlrate = 
			rfield.otslin[
			HydroLines[0][n][n-1].ipCont-1]*crs;

		hlmaxCom.hlmax[n] = (float)MAX2(hlmaxCom.hlmax[n],hlrate/MAX2(1e-30,
		  HPhoto.hgamnc[0][n]+phycon.eden*hcbr.hcoldn[0][n]+hcbr.hradn[0][n]));
	}

	hydro.FreeFreeTotHeat += hydro.FreeFreeHeat*(float)radius.dVeff;

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

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