/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*InitDefaultsPreparse initialization at start of simulation, called from cloudy
* before parser, sets initial values of quantities changed by parser 
* called for each point in a grid but one time in multi-iteration sims */
#include "cddefines.h" 
#include "physconst.h"
#include "phycon.h"
#include "radius.h"
#include "trace.h"
#include "dynamics.h"
#include "geometry.h"
#include "noexec.h"
#include "rfield.h"
#include "opacity.h"
#include "stopcalc.h"
#include "atomfeii.h"
#include "plot.h"
#include "rt.h"
#include "fudgec.h"
#include "abund.h"
#include "h2.h"
#include "ionbal.h"
#include "hextra.h"
#include "wind.h"
#include "atmdat.h"
#include "pressure.h"
#include "thermal.h"
#include "continuum.h"
#include "optimize.h"
#include "save.h"
#include "dense.h"
#include "hcmap.h"
#include "prt.h"
#include "init.h" 
#include "taulines.h"
#include "lines_service.h"
#include "grid.h"

/*InitDefaultsPreparse initialization at start of simulation, called from cloudy
 * before parser, will now be called one time per sim in grid but long term
 * goal is to call parser only one time in a grid.  This would be called 
 * first to set defaults that may be changed by the parser.  This routine
 * should only set defaults for variables that change in the parser.
 * It does not initialize variables for the start of a calculation.  That
 * is done by InitSimPostparse 
 * called one time in multi iteration sim */
void InitDefaultsPreparse( void )
{
	long int i, 
		ipISO,
		nelem;

	DEBUG_ENTRY( "InitDefaultsPreparse()" );

	/* init vars before parsing commands for each sim */

	/* option to turn off collisional ionization with "no collisional ionization" cmmnd */
	atmdat.lgCollIonOn = true;

	// default, Chianti not on, but if on use hybrid, do not print
	atmdat.lgChiantiOn = true;
	atmdat.lgChiantiHybrid = true;
	atmdat.lgChiantiPrint = false;
	//Use gbar to fill in dBase transitions if they lack collision strengths
	atmdat.lgGbarOn = true;
	//Tells Cloudy to exclusively use experimental energies in Chianti.
	atmdat.lgChiantiExp = true;
	// Set the default number of Chianti energy levels to use for Fe for photoionization case
	atmdat.nChiantiMaxLevelsFe = atmdat.nChiantiPhotoLevelsFe;
	// Set the default number of Chianti energy levels to use for all other elements
	atmdat.nChiantiMaxLevels = atmdat.nChiantiPhotoLevels;
	// ChiantiLevelsSet is false until the user specifies the number of chianti levels to use.
	atmdat.lgChiantiLevelsSet = false;
	/* nChiantiMaxLevelsFe and nChiantiMaxLevels are defaulted to 100 and 50 respectively if
	 * the coronal command is used. See parse_coronal.cpp */
	/* Set the default number Stout max levels
	 * For now use the default Chianti Fe levels = 100 */
	atmdat.nStoutMaxLevels = atmdat.nChiantiMaxLevelsFe;

	/* Stout on by default, optional printout off */
	atmdat.lgStoutOn = true;
	atmdat.lgStoutPrint = false;
	atmdat.lgStoutHybrid = true;

	// By default use Dima (Voronov97) data for collisional ionization rate coefficients
	// possible options are DIMA or HYBRID
	atmdat.CIRCData = t_atmdat::HYBRID;

	/** Set the default collision strength for dBase transitions when there is no radiative data */
	atmdat.collstrDefault = 0.01;

	atmdat.lgLamdaOn = true;
#	ifdef USE_CDMS
	atmdat.lgCalpgmOn = true;
#	else
	atmdat.lgCalpgmOn = false;
#	endif
	strcpy(atmdat.chCloudyChiantiFile, "CloudyChianti.ini");
	strcpy(atmdat.chStoutFile, "Stout.ini");

	/* drChange was reset to get orion flux in h-beta correct
	 * drChange is really tau of current zone */
	radius.drChange = 0.15f;

	radius.glbdst = 0.;
	radius.glbrad = 0.;

	// energy conservation check is slow, turn on with set check energy every zone
	continuum.lgCheckEnergyEveryZone = false;

	/* option to read in all input quantities and NOT execute the actual model
	* only check on input parameters - set by calling cdNoExec */
	noexec.lgNoExec = false;

	/* constant for the extinguish command */
	rfield.ExtinguishColumnDensity = 0.;
	rfield.ExtinguishLeakage = 0.;
	rfield.ExtinguishLowEnergyLimit = 1.;

	/* parameters having to do with thermal map */
	hcmap.RangeMap[0] = 10.f;
	hcmap.RangeMap[1] = .99e10f;
	/* zone where map is to be done */
	hcmap.MapZone = -1;
	hcmap.nMapStep = 20;

	thermal.ConstGrainTemp = 0.;
	thermal.lgTemperatureConstant = false;
	thermal.lgTemperatureConstantCommandParsed = false;
	thermal.ConstTemp = 0.;
	thermal.lgTeHigh = false;
	thermal.lgTeBD96 = false;
	thermal.lgTLaw = false;
	thermal.lgTeSN99 = false;
	/* try toe predict next zone's temperature in constant density models,
	* as done in ZoneStart.  Turned off with no tepred command */
	thermal.lgPredNextTe = true;

	/* turbulent heating - set with hextra command */
	hextra.TurbHeat = 0.;
	/* set true with TIME option on hextra command for time dependent sims */
	hextra.lgTurbHeatVaryTime = false;
	/* options for for extra heating to depend on scale radius */
	hextra.lgHextraDepth = false;
	/* options for for extra heating to depend on density */
	hextra.lgHextraDensity = false;
	/* options for alpha disk model heating */
	hextra.lgHextraSS = false;

	/* options set by cosmic ray command */
	hextra.cryden = 0.;
	hextra.cryden_ov_background = 0.;
	hextra.lg_CR_B_equipartition = false;
	hextra.crtemp = 0.;
	hextra.crpowr = 0.;
	hextra.cr_energydensity = 0;

	/* set with coronal equilibrium init time command */
	dynamics.lg_coronal_time_init = false;
	dynamics.lgTracePrint = false;

	/* parameters to do with wind */
	wind.lgWindOK = true;
	wind.DiskRadius = 0;
	wind.lgDisk = false;
	wind.windv0 = 0.;
	wind.setStatic();
	wind.comass = 0.;
	wind.windv = 0.;
	wind.dvdr = 0.;
	wind.emdot = 0.;
	wind.AccelAver = 0.;
	wind.acldr = 0.;
	wind.AccelGravity = 0.;
	wind.AccelTotalOutward = 0.;
	wind.AccelCont = 0.;
	wind.AccelElectron = 0.;
	wind.AccelLine = 0.;
	wind.AccelMax = 0.;
	wind.fmul = 0.;
	wind.lgVelPos = true;

	/* argument on ELEMENT LIMIT OFF XXX command, lowest abundance 
	 * element to include in the calculation */
	dense.AbundanceLimit = 0.;

	/* controls density fluctuations, when true variations are due to density changing,
	* when false are due to abundances changing if changes are on */
	dense.lgDenFlucOn = true;
	dense.lgDenFlucRadius = true;
	dense.flong = 0.;
	dense.flcPhase = 0.;

	/* this says keep initial density constant, 
	* so pressure from iter to iter not really const */
	dense.lgDenseInitConstant = true;

	// pressure does not vary with time by default
	dense.lgPressureVaryTime = false;
	//  required number is timescale for time variation
	dense.PressureVaryTimeTimescale = -1.;
	//  optional number is index for time variation
	dense.PressureVaryTimeIndex = 0.;

	/* extra electron density, set with eden command */
	dense.EdenExtra = 0.;

	/* forced electron density, set with set eden command */
	dense.EdenSet = 0.;

	/* option to set electron fraction, n_e/n_H */
	dense.EdenFraction = 0.;

	/* individual terms for the pressure equation of state */
	/* >>chng 05 jan 01, all three are set true at start *
	* code default is constant density, so all this is ignored
	* is turned to constant pressure then these flags must be adjusted so
	* that we get the pressure we expect.  with all true, all three
	* contributors to pressure will be counted - with constant gas pressure
	* these are turned off */
	pressure.lgPres_radiation_ON = true;
	pressure.lgPres_magnetic_ON = true;
	pressure.lgPres_ram_ON = true;
	/* flag for constant pressure, include continuum too */
	pressure.lgContRadPresOn = true;
	/* constant density is the default */
	strcpy( dense.chDenseLaw, "CDEN" );
	/* number on line is log of nT - option to specify initial pressure */
	pressure.lgPressureInitialSpecified = false;
	/* this is log of nT product - if not present then set zero */
	pressure.PressureInitialSpecified = 0;

	/* select certain atomic data, changed with set atomic data command 
	* this says to use Zeippen 1982 [OII] transition probabilities */
	dense.lgAsChoose[ipOXYGEN][1] = false;

	abund.lgAbnSolar = false;

	/* the ipSolar array is a set of pointers used for reading
	* in abundances with the "abundances" command
	* the hydrogen abundance of unity is not read in */
	abund.npSolar = LIMELM - 1;
	for( i=0; i < abund.npSolar; i++ )
	{
		abund.ipSolar[i] = i + 2;
	}

	/* option to turn off an element */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* set of scale factors for changing abundances with elements command */
		abund.ScaleElement[nelem] = 1.;
		abund.solar[nelem] = abund.SolarSave[nelem];

		// default scale factors for SET DIELECTRONIC RECOMBINATION KLUDGE SCALE
		ionbal.DR_mean_scale[nelem] = 1.;
	}

	abund.lgAbTaON = false;

	/* option to turn off an element */
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		/* option to have abundances from table */
		abund.lgAbunTabl[nelem] = false;
	}

	/* threshold for faintest heating cooling to save with save heating or 
	 * save cooling commands, reset with PUNCH WEAKHEATCOOL command */
	save.WeakHeatCool = 0.05f;

	/* set of variables used to control save command */
	save.nsave = 0;
	save.lgPunContinuum = false;
	save.lgDRPLst = false;
	save.lgDRHash = true;
	save.lgTraceConvergeBaseHash = true;
	/* this is the string that will appear after each model in the save output,
	 * reset with the "set save hash" command */
	strcpy( save.chHashString , "###########################"  );
	/* save every one continuum point, set skip save command */
	save.ncSaveSkip = 1;
	/* flush file after every iteration - set save flush */
	save.lgFLUSH = false;
	for( i=0; i<LIMPUN; ++i )
	{
		save.lgHashEndIter[i] = true;
		/* set false for time dependent calculations*/
		save.lg_separate_iterations[i] = true;
		save.lgSaveEveryZone[i] = false;
		save.nSaveEveryZone[i] = -1;
		save.emisfreq[i].set( -1. );
		save.ipEmisFreq[i] = -1;
	}

	/* default is to conserve energy, reset with
	 * set save line width / resolution command */
	save.Resolution = realnum(-1.);
	save.ResolutionAbs = realnum(-1.);

	/* punch dominant rates variables */
	for( long i=0; i<LIMPUN; i++ )
		strcpy( save.chSpeciesDominantRates[i], "" );

	/* default no printing of optical depths, TooFaint is .1 */
	prt.lgPrtTau = false;
	prt.PrtTauFnt = 0.1f;
	prt.lgPrtShort = false;
	prt.TooFaint = 1e-3f;
	prt.lgFaintOn = true;
	prt.lgFntSet = false;
	prt.lgPrnLineCell = false;
	prt.nPrnLineCell = -1000;
	prt.lgPrtCitations = false;

	/* change angle of illumination 
	* this is angle away from the normal, so 0 is a normal ray, the default*/
	geometry.DirectionalCosin = 1.;
	geometry.size = 1.f;
	geometry.lgSizeSet = false;
	geometry.covaper = -1.f;

	/* if true then print main block of lines as array,
	 * set false with print lines column, will then
	 * do a single column of lines */
	prt.lgPrtLineArray = true;

	/* when printing a column this is option to print linear rather than log */
	prt.lgPrtLineLog = true;

	/* print ages */
	prt.lgPrnAges = false;

	/* print column densities */
	prt.lgPrintColumns = true;

	/* option to sort lines by wavelength, print sort command */
	prt.lgSortLines = false;

	prt.lgPrtMaser = false;
	prt.lgPrintTime = true;

	prt.lgPrtContIndices = false;
	prt.lgPrnPump = false;
	prt.lgPrnInwd = false;
	prt.lgPrnColl = false;
	prt.lgPrnHeat = false;
	/* >>chng 00 dec 08, these determine the standard items included in "nFnu", PvH */
	prt.lgSourceReflected = true;
	prt.lgSourceTransmitted = false;
	prt.lgDiffuseInward = true;
	prt.lgDiffuseOutward = true;
	prt.lgPrtLastIt = false;
	prt.lgOnlyZone = false;
	prt.lgOnlyHead = false;
	prt.lgPrtStart = false;
	prt.nstart = 0;
	/* print predictions from collapsed levels of iso sequences */
	prt.lgPrnIsoCollapsed = false;

	/* turn off printing of heating agents */
	prt.lgPrintHeating = false;

	/* flag saying to print all matrix elements for ionization balance 
	 * set with PRINT ARRAYS command */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		prt.lgPrtArry[nelem] = false;
	}

	/* print line flux at earth */
	prt.lgPrintFluxEarth = false;

	/* print line surface brightness, def sr, option arcsec */
	prt.lgSurfaceBrightness = false;
	prt.lgSurfaceBrightness_SR = true;

	/* print line cumulative sets true, print integrated line intensity over
	 * time in temp dependent simulation */
	prt.lgPrintLineCumulative = false;

	prt.nzdump = -100;

	trace.lgSecIon = false;
	trace.lgTrOvrd = true;
	trace.lgOpacBug = false;
	trace.nTrConvg = 0;
	trace.lgTr8446 = false;
	trace.lgTrLevN = false;
	trace.lgTrGant = false;
	trace.lgOptcBug = false;
	trace.lgTrace3Bod = false;
	trace.lgOTSBug = false;
	trace.lgESOURCE = false;
	trace.lgTraceMole = false;
	trace.lgHeatBug = false;
	trace.lgHeavyBug = false;
	trace.lgBug2nu = false;
	trace.lgDrBug = false;
	trace.lgWind = false;
	trace.lgPtrace = false;
	trace.lgDrv_cdLine = false;
	trace.lgDustBug = false;
	trace.lgComBug = false;
	trace.lgHeBug = false;
	trace.lgCarBug = false;
	trace.lgCalBug = false;
	trace.lgConBug = false;
	trace.lgNeBug = false;
	trace.lgFeBug = false;
	trace.lgHBug = false;
	trace.lgTrLine = false;
	trace.nznbug = 10000;
	trace.npsbug = 10000;
	trace.lgTrace = false;
	trace.lgPointBug = false;
	trace.lgNeonBug = false;
	trace.lgCoolTr = false;
	trace.lgTrDiff = false;
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		trace.lgIsoTraceFull[ipISO] = false;

	/* variables used in stop ... command */

	/* various criteria for stopping model */
	/* >>chng 04 dec 21, remove from here and init to 1e30 in zero */
	/*StopCalc.tauend = 0.;*/
	StopCalc.tauend = 1e30f;

	/* >>chng 05 nov 22 - NPA.  Stop calculation when fraction of oxygen frozen
	* out on grains ices gets too high - formation of ices */
	/*StopCalc.StopDepleteFrac = 0.99f;*/
	/* >>chng 05 dec 16, with revised ion solver logic, code should be able to
	* converge away from situation where ices have disturbed the chemistry and
	* net negative atomic abundances result.  now we say solution not converged and
	* soldier on 
	* this test should not be necessary */
	StopCalc.StopDepleteFrac = 1.02f;

	StopCalc.xMass = 0.;
	StopCalc.taunu = 0.;
	StopCalc.iptnu = -1;
	/* stopping AV */
	StopCalc.AV_extended = 1e30f;
	StopCalc.AV_point = 1e30f;
	/* highest allowed temperature */
	StopCalc.TempHiStopZone = (realnum)phycon.TEMP_LIMIT_HIGH;
	StopCalc.TempHiStopIteration = (realnum)phycon.TEMP_LIMIT_HIGH;

	/* the floor sets a limit to the temperature in the calculation -
	* if te falls below this, we do a constant temperature cloud at
	* this temperature */
	StopCalc.TeFloor = 0.;

	/* stop zone calculations when Te falls below this,
	 * TEMP_STOP_DEFAULT in cddefines.h and is 4000 */
	StopCalc.TempLoStopZone = (realnum)phycon.TEMP_STOP_DEFAULT;
	/* stop iterations, used to stop time dependent command */
	StopCalc.TempLoStopIteration = -1.;

	/* ending column densities */
	StopCalc.HColStop = COLUMN_INIT;
	StopCalc.colpls = COLUMN_INIT;
	StopCalc.colnut = COLUMN_INIT;
	StopCalc.col_h2 = COLUMN_INIT;
	StopCalc.col_h2_nut = COLUMN_INIT;
	StopCalc.col_H0_ov_Tspin = COLUMN_INIT;
	StopCalc.col_monoxco = COLUMN_INIT;

	StopCalc.lgStopSpeciesColumn = false;
	StopCalc.chSpeciesColumn[0] = '\0';
	StopCalc.col_species = COLUMN_INIT;

	/* stopping electron density */
	StopCalc.StopElecDensity = -COLUMN_INIT;

	/* stopping electron and molecular fractions */
	StopCalc.StopElecFrac = -FLT_MAX;
	StopCalc.StopHPlusFrac = -FLT_MAX;
	/* stopping molecular fraction has opposite sign - want to stop when 2H_2/NH gt this */
	StopCalc.StopH2MoleFrac = FLT_MAX;
	/* this flag says that 21cm line optical depth is the stop quantity */
	StopCalc.lgStop21cm = false;
	/* debugging aid - stop when conv.nTotalIoniz reaches this value */
	StopCalc.nTotalIonizStop = 0;
	/* stop when absolute value of velocity falls below this */
	StopCalc.StopVelocity = 0.;
	/* number of stop line commands entered */
	StopCalc.nstpl = 0;

	/* initialize some variables for the optimizer */
	optimize.nIterOptim = 400;
	optimize.OptGlobalErr = 0.10f;
	optimize.nEmergent.clear();
	optimize.chLineLabel.clear();
	optimize.wavelength.clear();
	optimize.errorwave.clear();
	optimize.ipobs.clear();
	optimize.xLineInt_Obs.clear(); 
	optimize.xLineInt_error.clear();	
	optimize.chTempLab.clear();
	optimize.ionTemp.clear();
	optimize.temp_obs.clear(); 
	optimize.temp_error.clear();
	optimize.chTempWeight.clear();
	optimize.chColDen_label.clear();
	optimize.ion_ColDen.clear();
	optimize.ColDen_Obs.clear();
	optimize.ColDen_error.clear();
	optimize.ContIndex.clear();
	optimize.ContEner.clear();
	optimize.ContNFnu.clear();
	optimize.ContNFnuErr.clear();
	optimize.nRangeSet = 0;
	strcpy( optimize.chOptRtn, "PHYM" );

	/* flags says what is to be matched */
	optimize.lgOptLum = false;
	optimize.lgOptDiam = false;
	optimize.lgOptimize = false;

	/* trace flag for optimization process */
	optimize.lgTrOpt = false;

	optimize.lgOptimFlow = false;
	optimize.optint = 0.;
	optimize.optier = 0.;
#	if defined(__unix) || defined(__APPLE__)
	optimize.lgParallel = !cpu.i().lgMPISingleRankMode();
#	else
	optimize.lgParallel = false;
#	endif
	if( optimize.lgParallel )
		optimize.useCPU = cpu.i().nCPU();
	else
		optimize.useCPU = 1;
	optimize.lgOptCont = false;

	/* the fudge factors command */
	fudgec.nfudge = 0;
	fudgec.lgFudgeUsed = false;
	for( i=0; i < NFUDGC; i++ )
		fudgec.fudgea[i] = 0.;

	EmLineZero( DummyEmis );
	TauZero( DummyEmis );
	DummyEmis.iRedisFun() = 0;
	DummyEmis.ipFine() = -1;
	DummyEmis.gf() = 0.;
	DummyEmis.damp() = 0.;
	DummyEmis.dampXvel() = 0.;
	DummyEmis.opacity() = 0.;
	DummyEmis.Aul() = 1e-30f;

	/* parameters dealing with printer plots */
	for( i=0; i < NDPLOT; i++ )
	{
		plotCom.lgPltTrace[i] = false;
	}

	/* this says what types of printer plots we will make */
	for( i=0; i < NDPLOT; i++ )
	{
		strcpy( plotCom.chPType[i], "NONE" );
	}
	plotCom.lgPlotON = false;

	/* following were block data logic */
	rt.lgStarkON = true;

	/* by default use Federman form of shielding function */
	rt.nLineContShield = LINE_CONT_SHIELD_FEDERMAN;

	/* parameters set with Case A and Case B commands */
	/* this is flag for turning on case b */
	opac.lgCaseB = false;

	/* this is separate flag for turning off collisions from n=2 */
	opac.lgCaseB_HummerStorey = false;

	/* this is separate flag for turning off excited state photoionization */
	opac.lgCaseB_no_photo = false;
	/* another case b option, turn off background opacities, no Pdest */
	opac.lgCaseB_no_pdest = false;

	/* smallest allowed line and Lya optical depths, reset with
	* Case B command */
	opac.tlamin = 1e-20f;

	/* taumin command minimum optical depths for lines default 1e-20 */
	opac.taumin = 1e-20f;

	/* set false with no induced processes */
	rfield.lgInducProcess = true;

	/* this flag says that CMB has been set */
	rfield.lgCMB_set = false;
	
	rfield.lgComptonOn = true;

	for( i=0; i < LIMSPC; i++ )
	{
		/* this is set true if particular continuum source can vary with time 
		* set true if TIME appears on intensity / luminosity command line */
		rfield.lgTimeVary[i] = false;
		/* most continua enter as a beam rather than isotropic */
		rfield.lgBeamed[i] = true;
		// default is radiation from the "illuminated" face
		rfield.Illumination[i] = Illuminate::FORWARD;
		// optical depth = normal optical depth * this scale factor,
		// is 1 / cos theta
		rfield.OpticalDepthScaleFactor[i] = 1.;
		/* default energy range is H-ionizing radiation */
		rfield.range[i][0] = HIONPOT;
		rfield.range[i][1] = rfield.egamry;
		rfield.RSFCheck[i] = 0.;
	}

	/* line overlap opacity, turn off with no fine opacity command */
	rfield.lgOpacityFine = true;

	/* lower and upper wavelength bounds, Angstroms, for the FeII continuum */
	FeII.feconwlLo = 1000.;
	FeII.feconwlHi = 7000.;
	FeII.nfe2con = 1000;

	/* this is the faintest the high-energy tail of the continuum be */
	rfield.FluxFaint = 1e-10f;

	/* >>chng 01 jul 26, moved next statement from below loop to avoid bug in gcc 2.95.3, PvH */
	/* default diffuse fields is outward only */
	strcpy( rfield.chDffTrns, "OU2" );
	rfield.lgOutOnly = true;

	/* flags for whether continuum is defined over all energies */
	rfield.lgMMok = true;
	rfield.lgHPhtOK = true;
	rfield.lgXRayOK = true;
	rfield.lgGamrOK = true;

	for( rfield.ipSpec=0; rfield.ipSpec < LIMSPC; rfield.ipSpec++ )
	{
		if( rfield.lgContMalloc[rfield.ipSpec] )
		{
			rfield.tNu[rfield.ipSpec].clear();
			rfield.tslop[rfield.ipSpec].clear();
			rfield.tFluxLog[rfield.ipSpec].clear();
			rfield.lgContMalloc[rfield.ipSpec] = false;
		}
	}

	grid.lgNegativeIncrements = false;
	grid.lgSaveXspec = false;

	return;
}