/*FirstDR derive thickness of first zone */
#include <float.h>
#include <string.h>
#include "cddefines.h"
#include "cdminmax.h"/* vfmin */
#define	Z	1.0001
#include "itercnt.h"
#include "wind.h"
#include "u.h"
#include "sdrmin.h"
#include "stopcalc.h"
#include "trace.h"
#include "insane.h"
#include "ipdr.h"
#include "qs.h"
#include "drmnon.h"
#include "rfield.h"
#include "drmin.h"
#include "showme.h"
#include "hcaseb.h"
#include "fluct.h"
#include "filfac.h"
#include "phycon.h"
#include "radius.h"
#include "opac.h"
#include "ipoint.h"
#include "h.h"
#include "didz.h"
#include "densty.h"
#include "pressure.h"
#include "tabden.h"
#include "firstdr.h"

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

	int lgDoPun;

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

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

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

	if( wind.windv > 0. )
	{
		/* evaluate total pressure, although value not used (so stuffed into dr912) */
		dr912 = TotalPressure();
		accel = 1.3e-10*densty.pden*densty.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/(h.hi*filfac.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(phycon.hden*filfac.FillFac* 20.);
	}
	else if( StopCalc.colpls < 5e29 )
	{
		/* ionized column density */
		drcol = log10(StopCalc.colpls) - log10(h.hii*filfac.FillFac* 20.);
	}
	else if( StopCalc.colnut < 5e29 )
	{
		/* neutral column denisty */
		drcol = log10(StopCalc.colnut) - log10(h.hi*filfac.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( fluct.flong != 0. )
	{
		/* flong set => density fluctuations */
		dradf = 6.283/fluct.flong/10.;
		dradf = vfmin(dradf,radius.router[IterCnt.iter-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( (qs.qhtot>0.) && (qs.qhtot> qs.qbal*0.01) && (u.uh>1e-10) )
	{
		/* >> chng 99 dec 23, double to allow lte.in to work on alphas */
		drStromgren = (double)(qs.qhtot)/(double)(hcaseb.HRecRad[0]) /POW2(phycon.hden);

		/* different logic if this is a sphere */
		if( drStromgren/radius.rinner > 1. )
		{
			drStromgren = qs.qhtot*3./(radius.rinner*hcaseb.HRecRad[0] *POW2(phycon.hden) );
			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.opac[i] > BigOpacity )
		{
			BigOpacity = opac.opac[i];
			BigOpacityAnu = rfield.anu[i];
		}
	}

	/* this cannot possibly happen */
	assert( BigOpacity > 0. );

	/* 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*/
	drOpacity = (didz.drChange/100.)/BigOpacity/filfac.FillFac;

	/* initialize something in following routine? vars not used here */
	/* ContRate(&freqm,&opacm); 98 jan 24, kill, see what happens */

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

	drthrm = 1.5e31/MAX2(1.,phycon.hden*phycon.hden);

	/***********************************************************************
	 *
	 * 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(pressure.chCPres,"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 = phycon.hden/
				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;
	}

	/***********************************************************************
	 *
	 * 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[IterCnt.iter-1]/10.   ,   drcol,dr912,
		drthrm  ,  winddr  ,  dradf  ,  drTabDen  ,  FEND);

	/* option to set lower limit to zone thickness, with set drmin command*/
	radius.drad = MAX2( radius.drad , sdrminCom.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*/
	if( radius.drad == sdrminCom.sdrmin )
	{
		sdrminCom.lgDR2Big = TRUE;
	}
	else
	{
		sdrminCom.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( sdrminCom.sdrmax  ,  radius.drad );

	/* drMinimum is smallest acceptable DRAD, and is 1/100 OF DRAD(1) */
	/* this can be turned off by GLOB command */
	if( drmnon.lgDrMnOn )
	{
		drmin.drMinimum = (float)(radius.drad/1e7);
	}
	else
	{
		drmin.drMinimum = 0.;
	}

	if( sdrminCom.lgSMinON )
	{
		drmin.drMinimum = sdrminCom.sdrmin;
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" FirstDR called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n", 
		  radius.drad, drmin.drMinimum, sdrminCom.sdrmin, sdrminCom.sdrmax );
	}

	if( radius.drad < SMALLFLOAT*1.1 )
	{
		fprintf( ioQQQ, 
			" FirstDR detected likely insanity, found dr=%13.5e \n", radius.drad);
		fprintf( ioQQQ, 
			" FirstDR: calculation continuing but crash is likely. \n");
		/* this sets flag that insanity has occurred */
		insane();
		ShowMe();
	}

	
	/* 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( ipdr.lgDROn )
	{
		lgDoPun = TRUE;
	}
	else
	{
		lgDoPun = FALSE;
	}

	/* punch what we decided up? */
	if( lgDoPun )
	{
		if( radius.drad == drOpacity  )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from drOpacity, opac was %.2e at %.2e Ryd\n", 
			   radius.drad , BigOpacity , BigOpacityAnu );
		}
		else if( radius.drad == radius.Radius/20. )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from radius.Radius\n", 
			   radius.drad );
		}
		else if( radius.drad == drStromgren )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from drStromgren\n", 
			   radius.drad );
		}
		else if( radius.drad == radius.router[IterCnt.iter-1]/10. )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from radius.router[IterCnt.iter-1]\n", 
			   radius.drad );
		}
		else if( radius.drad == drcol )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from drcol\n", 
			   radius.drad );
		}
		else if( radius.drad == dr912 )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from dr912\n", 
			   radius.drad );
		}
		else if( radius.drad == sdrminCom.sdrmax )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from sdrminCom.sdrmax\n", 
			   radius.drad );
		}
		else if( radius.drad == drthrm )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from drthrm\n", 
			   radius.drad );
		}
		else if( radius.drad == winddr )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from winddr\n", 
			   radius.drad );
		}
		else if( radius.drad == dradf )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from dradf\n", 
			   radius.drad );
		}
		else if( radius.drad == drTabDen )
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR keys from drTabDen\n", 
			   radius.drad );
		}
		else
		{
			fprintf( ipdr.ipDRout, 
				" %11.3e FirstDR insanity\n", 
			   radius.drad );
			ShowMe();
		}
	}

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