/* 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_first derive thickness of first zone, called after conditions in first zone 
 * are established, sets 
radius.drad_x_fillfac
radius.drad
 */
#include "cddefines.h"
#define	Z	1.0001
#include "wind.h"
#include "stopcalc.h"
#include "thermal.h"
#include "trace.h"
#include "punch.h"
#include "pressure.h"
#include "iso.h"
#include "h2.h"
#include "rt.h"
#include "rfield.h"
#include "dense.h"
#include "hmi.h"
#include "geometry.h"
#include "opacity.h"
#include "ipoint.h"
#include "tabden.h"
#include "radius.h"

void radius_first(void)
{
	long int i ,
		ip;

	int lgDoPun;

	double accel, 
	  BigOpacity, 
	  BigOpacityAnu,
	  change,
	  dr912, 
	  drH2 ,
	  drContPres ,
	  drOpacity ,
	  drStromgren, /* used for Stromgren length */
	  drTabDen, 
	  dradf, 
	  drcol, 
	  drthrm, 
	  factor, 
	  winddr;

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

	/***********************************************************************
	 *
	 * wind model, use acceleration length                          
	 *
	 ***********************************************************************/

	if( wind.windv > 0. )
	{
		/* evaluate total pressure, although value not used (so stuffed into dr912) */
		/* >>chng 01 nov 02, remove call, confirm vals defined with assert */
		ASSERT( dense.pden > 0. && dense.wmole > 0. );
		accel = 1.3e-10*dense.pden*dense.wmole;
		winddr = POW2(wind.windv)/25./accel;
	}
	else
	{
		winddr = 1e30;
	}

	/* key off of Lyman continuum optical depth */
	if( StopCalc.taunu > 0.99 && StopCalc.taunu < 3. )
	{
		dr912 = StopCalc.tauend/6.3e-18/(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac)*Z/50.;
	}
	else
	{
		dr912 = 1e30;
	}

	/***********************************************************************
	 *
	 * key off of column density; total, neutral, or ionized                          
	 *
	 ***********************************************************************/

	if( StopCalc.HColStop < 5e29 )
	{
		/* this is useful for very thin columns, normally larger than 1st DR */
		drcol = log10(StopCalc.HColStop) - log10(dense.gas_phase[ipHYDROGEN]*geometry.FillFac* 20.);
	}
	else if( StopCalc.colpls < 5e29 )
	{
		/* ionized column density */
		drcol = log10(StopCalc.colpls) - log10(dense.xIonDense[ipHYDROGEN][1]*geometry.FillFac* 20.);
	}
	else if( StopCalc.colnut < 5e29 )
	{
		/* neutral column denisty */
		drcol = log10(StopCalc.colnut) - log10(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac*50.);
	}
	else
	{
		/* not used */
		drcol = 30.;
	}
	/* finally convert the drived column density to linear scale */
	drcol = pow(10.,MIN2(35.,drcol));

	/***********************************************************************
	 *
	 * key off of density or abundance fluctuations, must be small part of wavelength                          
	 *
	 ***********************************************************************/

	if( dense.flong != 0. )
	{
		/* flong set => density fluctuations */
		dradf = 6.283/dense.flong/10.;
		dradf = vfmin(dradf,radius.router[iteration-1]*
		  Z,drcol,dr912,FEND);
	}
	else
	{
		dradf = FLT_MAX;
	}

	/* >>>chng 99 nov 18, add check on stromgren length */
	/* estimate Stromgren length, but only if there are ionizing photons
	 * and not constant temperature model */
	if( (rfield.qhtot>0.) && (rfield.qhtot> rfield.qbal*0.01) && (rfield.uh>1e-10) )
	{
		/* >>chng 99 dec 23, double to allow lte.in to work on alphas */
		/* >>chng 03 mar 15, density to double to avoid overflow, PvH */
		drStromgren = (double)(rfield.qhtot)/iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]/
			POW2((double)dense.gas_phase[ipHYDROGEN]);

		/* different logic if this is a sphere */
		if( drStromgren/radius.rinner > 1. )
		{
			/* >>chng 03 mar 15, to double to avoid FP overflow, PvH */
			drStromgren = (double)rfield.qhtot*3./(radius.rinner*
			    iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*POW2((double)dense.gas_phase[ipHYDROGEN]) );
			drStromgren += 1.;
			/* this results in r_out / r_in */
			drStromgren = pow( drStromgren , 0.33333);
			/* make it a physics thickness in cm */
			drStromgren *= radius.rinner;
		}

		/* take one hundredth of this */
		drStromgren /= 100.;
	}
	else
	{
		drStromgren = FLT_MAX;
	}

	/***********************************************************************
	 *
	 * find largest opacity, to keep the first zone optical depth 1
	 * this is usually the physics that sets the first zone thickness
	 *
	 ***********************************************************************/

	/* >>>chng 99 jun 25, this is to simulate behavior of code before extension
	 * of continuum arrray to 1e-8 ryd */
	ip = ipoint(1e-5);

	/* find largest opacity */
	BigOpacity = 0.;
	BigOpacityAnu = -1.;
	for( i=ip; i < rfield.nflux; i++ )
	{
		/* remember largest opacity, and energy where this happened,
		 * make sure flux at energy is gt 0, can be zero for case where
		 * nflux increased to include emission from some ions */
		if( rfield.flux[i]>0. && opac.opacity_abs[i] > BigOpacity )
		{
			BigOpacity = opac.opacity_abs[i];
			BigOpacityAnu = rfield.anu[i];
		}
	}
	/* BigOpacity may be zero on very first call */

	/* drChange set with set didz command, is only number set with this command,
	 * default in zerologic is 0.15 
	 * set drad to small part of*/
	if( BigOpacity > SMALLFLOAT )
	{
		drOpacity = (radius.drChange/100.)/BigOpacity/geometry.FillFac;
	}
	else
	{
		drOpacity = 1e30;
	}

	/***********************************************************************
	 *
	 * thermalization length of typical lines
	 *
	 ***********************************************************************/

	drthrm = 1.5e31/MAX2(1.,POW2((double)dense.gas_phase[ipHYDROGEN]));

	/***********************************************************************
	 *
	 * make sure we resolve initial structure in tabden command
	 * if interpolated table we need to make sure that we resolve the 
	 * initial changes in the structure
	 *
	 ***********************************************************************/

	if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
	{
		drTabDen = 1.;
		i = 1;
		factor = 0.;
		while( i < 100 && factor < 0.05 )
		{
			/* check densities at ever larger dr's, until factor becomes more than 5% */
			factor = dense.gas_phase[ipHYDROGEN]/
				tabden(radius.Radius+drTabDen, drTabDen );
			/* density change can be positive or negative sign */
			factor = fabs(factor-1.);
			drTabDen *= 2.;
			i += 1;
		}
		drTabDen /= 2.;
	}
	else
	{
		drTabDen = 1e30;
	}

	/* >>chng 03 mar 20, add check on lyman band optical depth - want first zone
	 * to be thin in H2 bands */
	/* some tests are fully molecular with solomon process turned off,
	 * do not sense this when already almost fully molecular */
	if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 )
	{
		change = 0.1;
	}
	else
	{
		/* >>chng 04 mar 14, this branch, H is quite molecular,
		 * still do not want large changes in solomon rate since linearization
		 * would not work in hmole network, bu do not need such fine steps */
		change = 1.;
	}

	/* >>chng 04 mar 14 go back to original logic since molecular
	 * pdr's had big jump in conditions from
	 * first to second zon even when most H in H2
	change = 0.1; */
	/* >>chng 04 apr 18, change from 0.1 to 0.001, inital zones too large in 
	 * leiden test case f1 */
	change = 0.001;
	/* >>chng 04 mar 13, not too large when big H2 is on */
	if( h2.lgH2ON  && hmi.lgBigH2_evaluated )
	{
		if( fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > 0.05 )
		{
			/* changes in H2 heating caused by changes in solomon rate
			 * would drive temperature failures */
			/* >>chng 04 apr 18, change from 0.001 to 0.0001, inital zones too large in 
			 * leiden test case f1 */
			change = 0.0001;
		}
		else
		{
			/* >>chng 04 apr 18, change from 0.01 to 0.001, inital zones too large in 
			 * leiden test case f1 */
			change = 0.001;
		}
	}
	drH2 = change / SDIV( 
		hmi.H2_total * geometry.FillFac * hmi.H2Opacity );

	/* >>chng 06 feb 01, very high U ulirg models had dramatic increase in
	 * cont pre in first few zones,
	 * in constant total pressure case, don't want acceleration across first zone to
	 * be large compared with current gas pressure */
	if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn )
	{
		RT_radiative_acceleration();/**/
		drContPres = 0.05 * pressure.PresTotlInit / 
			(wind.AccelTot*dense.xMassDensity*geometry.FillFac*geometry.AngleIllum);
	}
	else
		drContPres = 1e30;

	/***********************************************************************
	 *
	 * we have now generated range of estimates of first thickness,
	 * now choose smallest of the group
	 *
	 ***********************************************************************/

	/* radius div by 20, to prevent big change in iron ionization for high ioniz gas,
	 * this is also the ONLY place that sphericity comes in */
	radius.drad = vfmin( drOpacity  , radius.Radius/20.,drStromgren ,
		radius.router[iteration-1]/10.   ,   drcol,dr912,
		drthrm  ,  winddr  ,  dradf  ,  drTabDen  ,  drH2 , drContPres , FEND);

	/* option to set lower limit to zone thickness, with set drmin command*/
	radius.drad = MAX2( radius.drad , radius.sdrmin );

	/* set flag for comment if the previous line forced a larger dr than
	 * would otherwise have been chosen.  will cause comment to be generated
	 * in PrtComment if set true*/
	/*lint -e777 float test equality */
	if( radius.drad == radius.sdrmin )
	{
		radius.lgDR2Big = TRUE;
	}
	else
	{
		radius.lgDR2Big = FALSE;
	}

	/* this min had been in the big min set above, but caused a false alarm
	 * on the lgDR2Big test above since the set dr command sets both in and max */
	radius.drad = MIN2( radius.sdrmax  ,  radius.drad );
	radius.drad_x_fillfac = radius.drad * geometry.FillFac;

	/* drMinimum is smallest acceptable DRAD, and is 1/100 OF DRAD(1) */
	/* this can be turned off by GLOB command */
	if( radius.lgDrMnOn )
	{
		/* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
		 * avoid false trigger across thermal fronts 
		 * add * dense.gas_phase */
		radius.drMinimum = (float)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7);
	}
	else
	{
		radius.drMinimum = 0.;
	}

	/* if set drmin is used, make sure drMinimum (which will cause an abort) is
	 * smaller than drmin */
	if( radius.lgSMinON )
	{
		/* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth
		 * avoid false trigger across thermal fronts 
		 * add * dense.gas_phase */
		radius.drMinimum = MIN2(radius.drMinimum * dense.gas_phase[ipHYDROGEN],(float)(radius.sdrmin/10.f) );
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n", 
		  radius.drad, radius.drMinimum/ dense.gas_phase[ipHYDROGEN], radius.sdrmin, radius.sdrmax );
	}

	if( radius.drad < SMALLFLOAT*1.1 )
	{
		fprintf( ioQQQ, 
			" PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad);
		fprintf( ioQQQ, 
			" radius_first: calculation continuing but crash is likely. \n");
		/* this sets flag that insanity has occurred */
		TotalInsanity();
	}
	
	/* all this is to only punch on last iteration
	 * the punch dr command is not really a punch command, making this necessary
	 * lgDRon is set true if "punch dr" entered */
	if( punch.lgDROn )
	{
		lgDoPun = TRUE;
	}
	else
	{
		lgDoPun = FALSE;
	}

	/* punch what we decided up? */
	if( lgDoPun )
	{
		/* create hash marks on second and later iterations */
		if( iteration > 1 && punch.lgHashEndIter )
		{
			static int iter_punch=-1;
			/*fprintf( punch.ipPnunit[ipPun], " ###########################\n" );*/
			if( iteration !=iter_punch )
				fprintf( punch.ipDRout, "%s\n",punch.chHashString );
			iter_punch = iteration;
		}
		/* this is common part of each line, the zone count, depth, chosen dr, and depth2go */
		/* >>chng 05 aug 15, had printed drNext, always zero, rather the drad, which is set here */
		fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go );

		if( radius.drad == drOpacity  )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n", 
			   BigOpacity , BigOpacityAnu );
		}
		else if( radius.drad == radius.Radius/20. )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from radius.Radius\n" );
		}
		else if( radius.drad == drStromgren )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from drStromgren\n");
		}
		else if( radius.drad == radius.router[iteration-1]/10. )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from radius.router[iteration-1]\n");
		}
		else if( radius.drad == drcol )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from drcol\n");
		}
		else if( radius.drad == radius.sdrmin )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from radius.sdrmin\n");
		}
		else if( radius.drad == dr912 )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from dr912\n");
		}
		else if( radius.drad == radius.sdrmax )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from radius.sdrmax\n");
		}
		else if( radius.drad == drthrm )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from drthrm\n");
		}
		else if( radius.drad == winddr )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from winddr\n");
		}
		else if( radius.drad == drH2 )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from H2 lyman lines\n");
		}
		else if( radius.drad == dradf )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from dradf\n");
		}
		else if( radius.drad == drTabDen )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from drTabDen\n");
		}
		else if( radius.drad == drContPres )
		{
			fprintf( punch.ipDRout, 
				"radius_first keys from radiative acceleration across zone\n");
		}
		else
		{
			fprintf( punch.ipDRout,  "radius_first insanity\n" );
			fprintf( ioQQQ, "radius_first insanity, radius is %e\n" ,
				radius.drad);
			ShowMe();
		}
	}
	/*lint +e777 float test equality */

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

