/* 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 */
/*Zerologic set initial values for logical variables 
 * this is called before any commands are parsed */
#ifdef _MSC_VER
#	pragma warning( disable : 4305 )/* disable const double to float warning in MS VS - 
									   very large number of warns result */
#endif
#include "cddefines.h"
#include "kshllenr.h"
#include "map.h"
#include "optimize.h"
#include "prt.h"
#include "extinc.h"
#include "plot.h"
/*  */
/* following used to enter total predicted continuum into emission line array 
 * actually entered into line array in lineset1
 */
#include "predcont.h"
/* energies where diffuse continuuum is to be entered into line array 
 * variables are declared in predcont.h, where array dimensioned 
 * NB - if these numbers change, the wavelength in the printout will change 
 * too, since the wavelength is derived form this */
/* >>>chng 99 mar 23, adjusted energies so that wavelength line list is
 * the same as it was in C90 - small changes were caused by going over
 * to proper rydberg constant */
float EnrPredCont[NPREDCONT] = {
	/* this is radio continuum at 3.4 cm */
	  2.680e-06
	, 7.445e-04 
	, 1.498e-03 
	, 2.211e-03 
	, 2.952e-03 
	, 3.677e-03 
	, 3.7501e-03 /* Ney-Allen */
	, 3.9915e-03 /* Ney-Allen */
	, 4.2543e-03 /* Ney-Allen */
	, 4.314e-03 
	, 4.6446e-03 /* Ney-Allen */
	, 5.162e-03 
	, 5.2462e-03 /* Ney-Allen */
	, 5.8079e-03 /* Ney-Allen */
	, 6.240e-03 
/* 	, 7.320e-03  */
	, 7.3312e-03 /* Ney-Allen */
	, 7.9936e-03 /* Ney-Allen */
	, 8.7119e-03 /* Ney-Allen */
	, 9.6125e-03 /* Ney-Allen */
	, 9.77243e-03
	, 1.1099e-02 /* Ney-Allen */
	, 1.2022e-02 /* Ney-Allen */
	, 1.29253e-02 
	, 2.2152e-02 
	, 3.92044e-02 
	, 5.54593e-02 
	/* next two either side of n=4 edge of hydrogen, set to 1.5% off either direction*/
	/* >>chng 00 sep 18, had been too close in energy */
	, 6.1563e-02 
	, 6.3437e-02 
	, 8.1460e-02 
	/* >>chng 00 sep 14, changed energies of paschen jump to be farther away as
	 * per note on BJ */
	, 0.1094 
	, 0.1128 
	, 0.14675 
	, 0.18653 
	/* >>chng 00 sep 14, next two energies changed since they were too close to BJ
	 * and so both ended up shortward of limit*/
	, 0.246   /* these two are the Balmer jump, below and above. */
	, 0.254  /* continuum binning not much better than 1% so offset energyies by more */
	, 0.375  /* peak on two photon continuum */
	, 0.38096 
	, 0.43994 
	, 0.44394 
	, 0.50811 
	, 0.57489
	, 0.62487 
	, 0.67155 
	, 0.70244 
	, 0.72163 
	, 0.74812 
	, 0.76172 
	, 0.77551 
	, 0.79681 
	, 0.81859 
	, 0.8260 
	, 0.84859 
	, 0.85618 
	, 0.87967 
	, 0.911267 /*exactly 1000A */ 
	/* points on either side of Lyman jump,
	 * energies changed to be robust when energy grid changes,
	 * grid resolution is about 1%, so change from 0.99783 and 1.000
	 * to 1 +/- 1.5%
	 * >>chng 00 sep 23 change wavelength points for next two */
	, 0.985 
	, 1.015 
	, 1.199 
	, 1.299 
	, 1.4984 
	, 1.58441 
	/* points on either side of Lyman jump,
	 * energies changed to be robust when energy grid changes,
	 * grid resolution is about 1%, so change from 1.80433 and 1.809
	 * to 1.807 +/- 1.5%
	 * >>chng 00 sep 23 change wavelength points for next two */
	, 1.780 
	, 1.834 
	, 2.283};

long int ipPredCont[NPREDCONT];
#include "heavy.h"
#include "fudgec.h"
#include "iterations.h"
#include "noexec.h"
#include "reason.h"
#include "rt.h"
#include "trace.h"
#include "magnetic.h"
#include "hextra.h"
#include "zero.h"
/*
 * long names of elements 
 * ok here since never change */
/*#include "elementnames.h"*/
struct /*t_elementnames*/ {

/* definition of this is in zerologic, where it is also initialized
 * THIS FILE MUST BE KEPT PARALLEL WITH THE DEFINITIOND IN zerologic.c */
	/* names are LIMELM+2 since c12o and c13o are last two */
	char chElementName[LIMELM+3][11] ;

	char chElementNameShort[LIMELM+3][5] ;

	/* symbols for elements, used in various printouts */
	char chElementSym[LIMELM+3][3] ;

	/* ionization stages */
	char chIonStage[LIMELM+4][3] ;

	/* string giving ionization stage as roman numerals */
	char chIonRoman[LIMELM+4][7];

} elementnames = {

	/*char chElementName[LIMELM+2][11] = */
	{
	 "Hydrogen  " ,
	 "Helium    " ,
	 "Lithium   " ,
	 "Beryllium " ,
	 "Boron     " ,
	 "Carbon    " ,
	 "Nitrogen  " ,
	 "Oxygen    " ,
	 "Fluorine  " ,
	 "Neon      " ,
	 "Sodium    " ,
	 "Magnesium " ,
	 "Aluminium " ,
	 "Silicon   " ,
	 "Phosphorus" ,
	 "Sulphur   " ,
	 "Chlorine  " ,
	 "Argon     " ,
	 "Potassium " ,
	 "Calcium   " ,
	 "Scandium  " ,
	 "Titanium  " ,
	 "Vanadium  " ,
	 "Chromium  " ,
	 "Manganese " ,
	 "Iron      " ,
	 "Cobalt    " ,
	 "Nickel    " ,
	 "Copper    " ,
	 "Zinc      " ,
	 "C12O18    " ,
	 "C13O18    " ,
	 "H2        " 
	},

	 /*char chElementNameShort[LIMELM][5] = */
	 {
 	 "HYDR" ,
	 "HELI" ,
	 "LITH" ,
	 "BERY" ,
	 "BORO" ,
	 "CARB" ,
	 "NITR" ,
	 "OXYG" ,
	 "FLUO" ,
	 "NEON" ,
	 "SODI" ,
	 "MAGN" ,
	 "ALUM" ,
	 "SILI" ,
	 "PHOS" ,
	 "SULP" ,
	 "CHLO" ,
	 "ARGO" ,
	 "POTA" ,
	 "CALC" ,
	 "SCAN" ,
	 "TITA" ,
	 "VANA" ,
	 "CHRO" ,
	 "MANG" ,
	 "IRON" ,
	 "COBA" ,
	 "NICK" ,
	 "COPP" ,
	 "ZINC" ,
	 "12CO" ,
	 "13CO" ,
	 "H2  " 
	 },

	/* symbols for elements, used in various printouts */
	 /*char elementnames.chElementSym[LIMELM][3] = */
	 { 
		"H " ,
		"He" ,
		"Li" ,
		"Be" ,
		"B " ,
		"C " ,
		"N " ,
		"O " ,
		"F " ,
		"Ne" ,
		"Na" ,
		"Mg" ,
		"Al" ,
		"Si" ,
		"P " ,
		"S " ,
		"Cl" ,
		"Ar" ,
		"K " ,
		"Ca" ,
		"Sc" ,
		"Ti" ,
		"V " ,
		"Cr" ,
		"Mn" ,
		"Fe" ,
		"Co" ,
		"Ni" ,
		"Cu" ,
		"Zn" ,
		"12" ,
	    "13" ,
	    "H2"} ,

	/* ionization stages */
	 /*char chIonStage[LIMELM+1][3] =*/
		{
		" 1" ,
		" 2" ,
		" 3" ,
		" 4" ,
		" 5" ,
		" 6" ,
		" 7" ,
		" 8" ,
		" 9" ,
		"10" ,
		"11" ,
		"12" ,
		"13" ,
		"14" ,
		"15" ,
		"16" ,
		"17" ,
		"18" ,
		"19" ,
		"20" ,
		"21" ,
		"22" ,
		"23" ,
		"24" ,
		"25" ,
		"26" ,
		"27" ,
		"28" ,
		"29" ,
		"30" ,
		"31" ,
		/* this is special for molecule */
		"CO",
		"CO",
		"H2"
		} ,

	 /* ionization stages as roman numeral */
	 /*char chIonRoman[LIMELM+1][7] =*/
		{
		"I" ,
		"II" ,
		"III" ,
		"IV" ,
		"V" ,
		"VI" ,
		"VII" ,
		"VIII" ,
		"IX" ,
		"X" ,
		"XI" ,
		"XII" ,
		"XIII" ,
		"XIV" ,
		"XV" ,
		"XVI" ,
		"XVII" ,
		"XVIII" ,
		"XIX" ,
		"XX" ,
		"XXI" ,
		"XXII" ,
		"XXIII" ,
		"XXIV" ,
		"XXV" ,
		"XXVI" ,
		"XXVII" ,
		"XXVIII" ,
		"XXIX" ,
		"XXX" ,
		"XXXI" ,
		/* this is special for molecule */
		"  ",
		"  ",
		"  "
		} 
	};

void Zerologic(void)
{
	long int i;

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

	/* parameters having to do with thermal map */
	map.RangeMap[0] = 10.;
	map.RangeMap[1] = .99e10;
	/* >>chng 04 jun 16, had been 9999, was accidentally hit once,
	 * reset to -1 to indicate not to be used */
	/*map.MapZone = 9999;*/
	map.MapZone = -1;
	map.lgMapOK = TRUE;
	map.lgMapDone = FALSE;
	map.nMapStep = 20;
	/* will be reset to positive number when map actually done */
	map.nMapAlloc = 0;
	map.nmap = 0;
	map.lgMapBeingDone = FALSE;

	/* upper limit to energies of inner shell opacities in ryd
	 * this is 1 mev by default */
	KshllEnr.EnergyKshell = 7.35e4;

	/* constant for the extinguish command */
	extinc.excolm = 0.;
	extinc.exleak = 0.;
	extinc.exlow = 1.;
	/* the constant that multiplies the column density to get optical depth at 1 Ryd */
	extinc.cnst_col2optdepth = 6.22e-18f ;
	/* the power on the energy */
	extinc.cnst_power = -2.43f ;

	/* list of shells, 1 to 7 */
	strcpy( Heavy.chShell[0], "1s" );
	strcpy( Heavy.chShell[1], "2s" );
	strcpy( Heavy.chShell[2], "2p" );
	strcpy( Heavy.chShell[3], "3s" );
	strcpy( Heavy.chShell[4], "3p" );
	strcpy( Heavy.chShell[5], "3d" );
	strcpy( Heavy.chShell[6], "4s" );

	/* zero out diffuse recombinaiton emission arrays */
	for( i=0; i < LIMELM; i++ )
	{
		int j;
		for( j=0; j < (LIMELM - 1); j++ )
		{
			Heavy.xLyaHeavy[i][j] = 0.;
		}
	}

	noexec.lgNoExec = FALSE;

	/* default no printing of optical depths, TooFaint is .1 */
	prt.lgPrtTau = FALSE;
	prt.PrtTauFnt = 0.1;
	prt.lgPrtShort = FALSE;
	prt.TooFaint = 1e-3;
	prt.lgFaintOn = TRUE;
	prt.lgFntSet = FALSE;

	/* 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.lgPrtCont = FALSE;
	prt.lgPrnDiff = 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;

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

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

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

	prt.nzdump = -100;

	fudgec.nfudge = 0;

	for( i=0; i < NFUDGC; i++ )
	{
		fudgec.fudgea[i] = 0.;
	}

	for( i=0; i < ITRDIM; i++ )
	{
		iterations.IterPrnt[i] = 10000;
	}
	iterations.itermx = 0;
	/* this implements set coverage command */
	iterations.lgConverge_set = FALSE;
	iteration = 0;

	hextra.TurbHeat = 0.;
	/* >>chng 01 feb 22, had set this to 1e30, change to 0 as sentinal that is not used */
	hextra.turrad = 0.;
	hextra.turback = 0.;
	hextra.cryden = 0.;
	hextra.cryden_ov_background = 0.;
	/* default cosmic ray background */
	/* >>chng 99 jun 24, slight change to value
	 * quoted by 
	 * >>refer	cosmic ray	ionization rate	McKee, C.M., 1999, astro-ph 9901370
	 * this will produce a total
	 * secondary ionization rate of 2.5e-17 s^-1, as tested in 
	 * tsuite cosmicray.in.  If each ionization produces 2.4 eV of heat,
	 * the background heating rate should be 9.6e-29 * n*/
	/* >>chng 04 jan 26, update cosmic ray ionization rate for H0 to
	 * >>refer	cosmicray	ionization	Williams, J.P., Bergin, E.A., Caselli, P., 
	 * >>refercon	Myers, P.C., & Plume, R. 1998, ApJ, 503, 689,
	 * H0 ionization rate of 2.5e-17 s-1 and a H2 ionization rate twice this
	 * >>chng 04 mar 15, comment said 2.5e-17 which is correct, but code produce 8e-17,
	 * fix back to correct value 
	 */
	/* NB - the rate is derived from the density.  these two are related by the secondary
	 * ionization efficiency problem.  background_rate is only here to provide the relationship
	 * for predominantly neutral gas.  the background_density is the real rate. 
	hextra.background_density = 1.99e-9f;*/
	/* >>chng 05 apr 16, to get proper ionization rate in ism_set_cr_rate, where
	 * H is forced to be fully atomic, no molecules, density from 1.99 to 2.15 */
	hextra.background_density = 2.15e-9f;
	hextra.background_rate = 2.5e-17f;
	hextra.crpowr = 0.;
	hextra.crtemp = 0.;

	/* parameters having to do with magnetic field */
	Magnetic_init();

	strcpy( reason.chReason, "reason not specified." );

	/* following were block data logic */
	rt.lgStarkON = TRUE;
	/*rt.nLineContShield = LINE_CONT_SHIELD_PESC;*/
	/* >>chng 04 jun 12, change from pesc to federman form of shielding function */
	rt.nLineContShield = LINE_CONT_SHIELD_FEDERMAN;

	optimize.lgVaryOn = FALSE;

	trace.lgSecIon = FALSE;
	trace.lgTrOvrd = TRUE;
	trace.lgTraceInput = FALSE;
	trace.lgOpacBug = FALSE;
	trace.lgTrConvg = 0;
	trace.lgTr8446 = FALSE;
	trace.lgTrLevN = FALSE;
	trace.lgTrGant = FALSE;
	trace.lgOptcBug = FALSE;
	trace.lgTrace3Bod = FALSE;
	trace.lgOTSBug = FALSE;
	trace.lgESOURCE = FALSE;
	trace.lgTr_CO_Mole = FALSE;
	trace.lgTr_H2_Mole = FALSE;
	trace.lgHeatBug = FALSE;
	trace.lgHeavyBug = FALSE;
	trace.lgBug2nu = FALSE;
	trace.lgDrBug = 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;
	trace.lgIsoTraceFull[ipH_LIKE] = FALSE;
	trace.lgIsoTraceFull[ipHE_LIKE] = FALSE;

	/* 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;

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