/* 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 */
/*ZoneEnd last routine called after all zone calculations, before iter_end_check,
 * upon exit radiation field is for outer edge of current zone */
/*ZoneStart set variables that change with each zone, like radius, depth,
 * upon exit flux will be flux at center of zone about to be computed */
#include "cddefines.h"
#include "phycon.h"
#include "opacity.h"
#include "rfield.h"
#include "struc.h"
#include "thermal.h"
#include "dense.h"
#include "h2.h"
#include "geometry.h"
#include "conv.h"
#include "dynamics.h"
#include "radius.h"
#include "zones.h"
/* this is number of zones to include in guess for next temperatures */
#define	IOFF	3

void ZoneStart(const char *chMode)
{
	int lgNeOscil,
		lgTeOscil;
	long int i;

	double dt1 , de1, 
	  /* this is correction factor for dilution of radiation
	   * across current zone, set in ZoneStart to correct
		* flux for center of current zone, and undone in ZoneDone */
	  DilutionCorrec ,
		drFac ,
		dTauThisZone,
	  outwrd, 
	  ratio, 
	  rm, 
	  rad_middle_zone,
	  vin, 
	  vout;

	static double DTaver , DEaver, 
	  dt2 , de2;

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

	/* this sub can be called in two ways, 'incr' means increment
	 * radius variables, while 'init' means only initialize rest */
	/* called at start of a zone to update all variables
	 * having to do with starting this zone */

	/* first establish current filling factor (normally 1) since used within
	 * following branches */
	geometry.FillFac = (float)(geometry.fiscal*
		pow( radius.Radius/radius.rinner ,(double)geometry.filpow));
	geometry.FillFac = (float)MIN2(1.,geometry.FillFac);

	if( strcmp(chMode,"init") == 0 )
	{
		/* initialize the variables at start of calculation, after this exits
		 * radius is equal to the initial radius, not outer edge of zone */

		/* depth to illuminated face */
		radius.depth = 1e-30;

		/* integral of depth times filling factor */
		radius.depth_x_fillfac = radius.depth*geometry.FillFac;
		/* effective radius */
		radius.drad_x_fillfac = radius.drad*geometry.FillFac;

		/* reset radius to inner radius, at this point equal to illuminated face radius */
		radius.Radius = radius.rinner;
		radius.Radius_mid_zone = radius.rinner + radius.drad/2.;

		/* thickness of next zone */
		radius.drNext = radius.drad;


		/* depth to illuminated face */
		radius.depth_mid_zone = radius.drad/2.;

		/* depth variable for globule */
		radius.glbdst = radius.glbrad;

		/* this is case where outer radius is set */
		if( radius.router[iteration-1] < 5e29 )
		{
			radius.Depth2Go = radius.router[iteration-1] ;
		}
		else if( iteration > 1 )
		{
			/* this case second or later iteration, use previous thickness */
			radius.Depth2Go = radius.router[iteration-2];
		}
		else
		{
			/* first iteration and depth not set, so we cannot estimate it */
			radius.Depth2Go = 0. ;
		}
	}
	else if( strcmp(chMode,"incr") == 0 )
	{
		/* update radius variables - called by cloudy at start of this zone's calcs */
		radius.drad = radius.drNext;
		ASSERT( radius.drad > 0. );
		radius.depth += radius.drad;
		radius.depth_mid_zone = radius.depth - radius.drad/2.;

		/* effective radius */
		radius.drad_x_fillfac = radius.drad*geometry.FillFac;

		/* integral of depth times filling factor */
		radius.depth_x_fillfac += radius.drad_x_fillfac;

		/* decrement Depth2Go but do not let become negative */
		radius.Depth2Go = MAX2( 0.,radius.Depth2Go - radius.drad );

		/* increment the radius, during the calculation Radius is the
		 * outer radius of the current zone*/
		radius.Radius += radius.drad*radius.dRadSign;

		/* Radius is now outer edge of this zone since incremented above, 
		 * so need to add drad to it */
		radius.Radius_mid_zone = radius.Radius - radius.drad/2.;

		/***********************************************************************
		 *
		 * attenuuate rfield to center of this zone
		 *
		 ***********************************************************************/

		 /* radius was just incremented above, so this resets continuum to
		  * flux at center of zone we are about to compute.  radius will be
		  * incremented immediately following this.  this correction is undone
		  * when ZoneDone called */

		/* this will be the optical thickness of the next zone
		 * AngleIllum is 1/COS(theta), is usually 1, reset with illuminate command,
		 * option for illumination of slab at an angle */

		/* radius.dRNeff is next dr with filling factor, this will only be
		 * used to get approximate correction for attenuation
		 * of radiation in this zone,
		 * equations drived in hazy, netzer&ferland 83, for factor accounting
		 * any changes here should be checked with both sphonly.in and pphonly*/
		/* honlyotspp seems most sensitive to the 1.35 */
		drFac = radius.drad*geometry.FillFac*geometry.AngleIllum*1.35;

		/* dilute radiation so that rfield is in center of zone, this
		 * goes into tmn at start of zone here, but is later divided out
		 * when ZoneEnd called */
		DilutionCorrec = 1./POW2(
			(radius.Radius-radius.drad/2.)/(radius.Radius-radius.drad) );

		/* note that this for loop is through <= nflux, to include the [nflux]
		 * unit integration verification token.  The token was set in
		 * in MadeDiffuse and propagated out in metdif.  the opacity at that energy is
		 * zero, so only the geometrical part of tmn will affect things.  The final
		 * sum is verified in PrtComment */
		for( i=0; i <= rfield.nflux; i++ )
		{
			dTauThisZone = opac.opacity_abs[i]*drFac;
			/* 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 */

			if( dTauThisZone < 1e-4 )
			{
				/* this small tau limit is the most likely branch, about 60% in parispn.in */
				opac.tmn[i] = 1.f;
			}
			else if( dTauThisZone < 5. )
			{
				/* this middle tau limit happens in the remaining 40% */
				opac.tmn[i] = (float)((1. - exp(-dTauThisZone))/(dTauThisZone));
			}
			else
			{
				/* this happens almost not at all */
				opac.tmn[i] = (float)(1./(dTauThisZone));
			}

			/* now add on correction for dilution across this zone */
			opac.tmn[i] *= (float)DilutionCorrec;

			/* actually do the corrections */
			rfield.flux[i] *= opac.tmn[i];

			/* >>chng 03 nov 08, update SummedCon here since flux changed */
			rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
		}

		/* following is distance to globule, usually doesnt matter */
		radius.glbdst -= (float)radius.drad;

		/* if gt 5th zone, and not constant pressure, and not oscillating te
		 * then guess next temperature : option to predict next temperature
		 * NZLIM is dim of structure variables saving temp, do data if nzone>NZLIM
		 * IOFF is number of zones to look over, set set to 3 in the define above */
		/* >>chng 01 mar 12, add option to not do this, set with no tepred command */
		if( (strcmp(dense.chDenseLaw,"CPRE") != 0) && 
			thermal.lgPredNextTe && 
			(nzone > IOFF + 1)  )
		{
			phycon.TeInit = phycon.te;
			phycon.EdenInit = dense.eden;
			lgTeOscil = FALSE;
			lgNeOscil = FALSE;
			dt1 = 0.;
			dt2 = 0.;
			de1 = 0.;
			de2 = 0.;
			DTaver = 0.;
			DEaver = 0.;
			for( i=nzone - IOFF-1; i < (nzone - 1); i++ )
			{
				dt1 = dt2;
				de1 = de2;
				/* this will get the average change in temperature for the 
				 * past IOFF zones */
				dt2 = struc.testr[i-1] - struc.testr[i];
				de2 = struc.ednstr[i-1] - struc.ednstr[i];
				DTaver += dt2;
				DEaver += de2;
				if( dt1*dt2 < 0. )
				{
					lgTeOscil = TRUE;
				}
				if( de1*de2 < 0. )
				{
					lgNeOscil = TRUE;
				}
			}

			/* option to guess next electron temperature if no oscillations */
			if( !lgTeOscil )
			{
				DTaver /= (double)(IOFF);
				/* don't want to overcorrect, use smaller of two */
				dt2 = fabs(dt2);
				rm = fabs(DTaver);
				DTaver = sign(MIN2(rm,dt2),DTaver);
				/* do not let it change by more than 5% of current temp */
				/* >>chng 03 mar 18, from 5% to 1% - convergence is much better
				 * now, and throwing the next Te off by 5% could easily disturb
				 * the solvers - primal.in was a good case */
				DTaver = sign(MIN2(rm,0.01*phycon.te),DTaver);
				/* this actually changes the temperature */
				phycon.te = phycon.te - DTaver;
				tfidle(TRUE);
			}
			else
			{
				/* temp was oscillating - too dangerous to do anything */
				DTaver = 0.;
			}
			/* this is used to remember the proposed temperature, for output
			 * in the punch predictor command */
			phycon.TeProp = phycon.te;

			/* option to guess next electron density if no oscillations */
			if( !lgNeOscil )
			{
				DEaver /= IOFF;
				de2 = fabs(de2);
				rm = fabs(DEaver);
				/* don't want to overcorrect, use smaller of two */
				DEaver = sign(MIN2(rm,de2),DEaver);
				/* do not let it change by more than 5% of current temp */
				DEaver = sign(MIN2(rm,0.05*dense.eden),DEaver);
				/* this actually changes the temperature */
				dense.eden -= DEaver;
			}
			else
			{
				/* temp was oscillating - too dangerous to do anything */
				DEaver = 0.;
			}
			/* this is used to remember the proposed temperature, for output
			 * in the punch predictor command */
			phycon.EdenProp = dense.eden;
			tfidle(FALSE);
		}
	}

	else
	{
		fprintf( ioQQQ, " PROBLEM ZoneStart called with insane argument, %4.4s\n", 
		  chMode );
		/* TotalInsanity exits after announcing a problem */
		TotalInsanity();
	}

	/* find the mean dr for the past few zones - this var is used
	 * in setting the optical scale in induced processes, and large
	 * changes in it will drive oscillations in the H+ fraction */
	radius.drad_x_fillfac_mean = radius.drad_x_fillfac;
	i = 1;
	while( i < 5 && nzone-i-1>=0 )
	{
		radius.drad_x_fillfac_mean += struc.drad_x_fillfac[nzone-i-1];
		++i;
	}
	radius.drad_x_fillfac_mean /= i;

	/* do advection if enabled */
	if( dynamics.lgAdvection )
		DynaStartZone();

	/* clear flag indicating the ionization convergence failures 
	 * have ever occured in current zone */
	conv.lgConvIonizThisZone = FALSE;

	/* this will say how many times the large H2 molecule has been called in this zone -
	 * if not called (due to low H2 abundance) then not need to update its line arrays */
	h2.nCallH2_this_zone = 0;

	/* check whether zone thickness is too small relative to radius */
	if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
	{
		ratio = radius.drad/(radius.glbdst + radius.drad);
		/* this only matters for globule command */
		if( ratio < 1e-14 )
		{
			radius.lgdR2Small = TRUE;
		}
		else
		{
			radius.lgdR2Small = FALSE;
		}
	}


	/* area factor, ratio of radius to out edge of this zone 
	 * relative to radius of illuminated face of cloud */
	/*radius.r1r0sq = (float)POW2(
		(radius.Radius - radius.drad*radius.dRadSign/2.)/radius.rinner);*/
	/*>>>chng 99 jan 25, definition of r1r0sq is now outer edge of current zone
	 * relative to inner edge of cloud */
	radius.r1r0sq = POW2( radius.Radius/radius.rinner );

	/* following only approximate, used for center of next zone */
	radius.dRNeff = radius.drNext*geometry.FillFac;

	/* this is the middle of the zone */
	rad_middle_zone = radius.Radius - radius.drad/2.;

	/* Radius is outer radius of this zone, so radius - drad is inner radius
	 * rinner is inner radius of nebula; dVeff is the effective vol element
	 * dVeff is vol rel to inner radius, so has units cm
	 * for plane parallel dVeff is dReff */
	if( radius.drad/radius.Radius > 0.01 )
	{
		vin = ((radius.Radius - radius.drad)/radius.rinner)*
			(MIN2((radius.Radius-radius.drad),radius.CylindHigh)/radius.rinner)*
			(radius.Radius - radius.drad);

		vout = (radius.Radius/radius.rinner)*
		  (MIN2(radius.Radius,radius.CylindHigh)/radius.rinner)*
		  radius.Radius;

		/* default of iEmissPower is 2, observing the full sphere */
		/* this is the usual case the full volume, just the difference in the two vols */
		/* this will become the true vol when multiplied by the square of the inner radius */
		radius.dVeff = (vout - vin)/3.*geometry.FillFac;
	}
	else
	{
		/* thin cell limit */
		/* rad_middle_zone is the middle of the zone;
		 * radius is always the OUTER edge of the current zone */
		radius.dVeff = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac*
			(MIN2(rad_middle_zone,radius.CylindHigh)/radius.rinner);
	}

	/* this is possible correction for slit projected onto resolved object -
	 * this only happens when aperture command is entered.
	 * NB not clear what happens when slit and cylendar are both used -
	 * this would be very orientation specific */

	/* default of iEmissPower is 2, set to 0 is just aperture beam, 
	 * and to 1 if aperture slit set */
	if( geometry.iEmissPower == 1 )
	{
		/* this is long slit option, remove one power of radius */
		radius.dVeff /= (rad_middle_zone/radius.rinner);
	}
	else if( geometry.iEmissPower == 0 )
	{
		/* this is short slit option, remove two powers of radius */
		radius.dVeff /= POW2(rad_middle_zone/radius.rinner);
	}

	/* covering factor, default is unity */
	radius.dVeff *= geometry.covgeo;

	/* these are needed for line intensities in outward beam
	 * lgSphere set, geometry.covrt usually 1, 0 when not lgSphere
	 * so outwrd is 1 when lgSphere set 1/2 when not set, */
	outwrd = (1. + geometry.covrt)/2.;

	/*>>>chng 99 apr 23, from above to below */
	/*radius.dVolOutwrd = outwrd*POW2( (radius.Radius-radius.drad_x_fillfac/2.)/radius.Radius) * 
	  radius.drad;*/
	/* this includes covering fact, the effective vol,, and 1/r^2 dilution,
	 * dReff includes filling factor.  the covering factor part is 1 for sphere,
	 * 1/2 for open */
	/*radius.dVolOutwrd = outwrd*radius.Radius*radius.drad_x_fillfac/(radius.Radius + 
	  2.*radius.drad);*/
	/* dReff from above, includes filling factor */
	radius.dVolOutwrd = outwrd*radius.drad_x_fillfac;
	ASSERT( radius.dVolOutwrd > 0. );

	/* following will be used to "uncorrect" for this in lines when predicting continua
	radius.GeoDil = radius.Radius/(radius.Radius + 2.*radius.drad);*/

	/* this should multiply the line to become the net inward.  geo part is 1/2 for
	 * open geomery, 0 for closed.  for case of isotropic emitter only this and dVolOutwrd
	 * above are used */
	radius.dVolReflec = (1. - outwrd)*radius.drad_x_fillfac*radius.r1r0sq;

	if( geometry.lgSphere )
	{
		/* inward beams do not go in when lgSphere set since geometry symmetric */
		radius.BeamInIn = 0.;
		radius.BeamInOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius + 
		  2.*radius.drad);
	}
	else
	{
		radius.BeamInIn = radius.drad_x_fillfac*radius.r1r0sq;

		/* inward beams do not go out */
		radius.BeamInOut = 0.;
	}
	/* factor for outwardly directed part of line */
	radius.BeamOutOut = radius.Radius*radius.drad_x_fillfac/(radius.Radius + 
	  2.*radius.drad);


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

void ZoneEnd(void)
{
	long i;

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

	/***********************************************************************
	 *
	 *  correct rfield for attenuation from center of zone to inner edge
	 *
	 ***********************************************************************/

	 /* radius is outer radius of this zone, this resets continuum to
	  * flux at illuminated face of zone we have already computed. */

	/* opac.tmn defined in ZoneStart */
	/* undilute radiation so that rfield is at face of zone */
	/* NB, upper limit of sum includes [nflux] to test unit verification cell */
	for( i=0; i <= rfield.nflux; i++ )
	{
		rfield.flux[i] /= opac.tmn[i];
		/* >>chng 03 nov 08, update SummedCon here since flux changed */
		rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
	}

	/* do advection if enabled */
	if( dynamics.lgAdvection )
		DynaEndZone();

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