/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/* */
#include "cddefines.h"
/* cddefines.h set EXTERN,
 * now unset EXTERN so that everything is defined here */
#ifdef EXTERN
#	undef EXTERN
#endif
#define EXTERN 
/* */
/* some things in cddefines.h */
/* this is the output cloudy will talk to, default is stdout and is set below */
FILE *ioQQQ ;
FILE *ioStdin ;
/* these are options to print error output 
 * requested with print errors command */
FILE* ioPrnErr; 
/* this is declared in cddefines and is used for io
 * saying when subs are entered and return */
FILE *debug_fp;
/* used for punching map*/
FILE *ioMAP = NULL;
/* following set TRUE with print errors command */
int lgPrnErr;
/* set true with SET ASSERT FPEE */
int lgAssertFPE;
/* levels within hydrogenic atoms */
const int ipH1s = 0;
const int ipH2s = 1;
const int ipH2p = 2;

/* some levels for he-like species, defined in cddefines.c */

/* level 1 */
const int ipHe1s1S = 0;

/* level 2 */
const int ipHe2s3S = 1;
const int ipHe2s1S = 2;
const int ipHe2p3P0 = 3;
const int ipHe2p3P1 = 4;
const int ipHe2p3P2 = 5;
const int ipHe2p1P = 6;

/* level 3 */
const int ipHe3s3S = 7;
const int ipHe3s1S = 8;
const int ipHe3p3P = 9;
const int ipHe3d3D = 10;
const int ipHe3d1D = 11;
const int ipHe3p1P = 12;

/* level 4 */
const int ipHe4s3S = 13;
const int ipHe4s1S = 14;
const int ipHe4p3P = 15;
const int ipHe4d3D = 16;
const int ipHe4d1D = 17;
const int ipHe4f3F = 18;
const int ipHe4f1F = 19;
const int ipHe4p1P = 20;

/* level 5 */
const int ipHe5s3S = 21;
const int ipHe5s1S = 22;
const int ipHe5p3P = 23;
const int ipHe5d3D = 24;
const int ipHe5d1D = 25;
const int ipHe5f3F = 26;
const int ipHe5f1F = 27;
const int ipHe5g3G = 28;
const int ipHe5g1G = 29;
const int ipHe5p1P = 30;

/* these are array indices for isoelectronic sequences,
 * same as element but used for array addressing to make
 * context totally clear */
const int ipH_LIKE=0;
const int ipHE_LIKE=1;
const int ipLI_LIKE=2;
const int ipBE_LIKE=3;
const int ipB_LIKE=4;
const int ipC_LIKE=5;
const int ipN_LIKE=6;
const int ipO_LIKE=7;
const int ipF_LIKE=8;

/* these are array indices for some elements */
const int ipHYDROGEN=0;
const int ipHELIUM=1;
const int ipLITHIUM=2;
const int ipBERYLLIUM=3;
const int ipBORON=4;
const int ipCARBON=5;
const int ipNITROGEN=6;
const int ipOXYGEN=7;
const int ipFLUORINE=8;
const int ipNEON=9;
const int ipSODIUM=10;
const int ipMAGNESIUM=11;
const int ipALUMINIUM=12;
const int ipSILICON=13;
const int ipPHOSPHORUS=14;
const int ipSULPHUR=15;
const int ipCHLORINE=16;
const int ipARGON=17;
const int ipPOTASIUM=18;
const int ipCALCIUM=19;
const int ipSCANDIUM=20;
const int ipTITANIUM=21;
const int ipVANADIUM=22;
const int ipCHROMIUM=23;
const int ipMANGANESE=24;
const int ipIRON=25;
const int ipCOBALT=26;
const int ipNICKEL=27;
const int ipCOPPER=28;
const int ipZINC=29;

/* external ZeroNum used to div by zero 
 * ok here since never changed*/
const double ZeroNum = 0.;

/* nzone is zone counter, incremented in routine cloudy */
long int nzone;
/* this is nzone + conv.nPres2Ioniz/100 in ConvBase */
double fnzone;

/* the iteration counter, set and incremented in routine cloudy */
long int iteration;

#include "taulines.h"/* this must go here since it defines NTA needed for other lines*/

/* following is true extern in taulines.h */
long nWindLine = NWINDDIM;
float *cs1_flag_lev2;

/* set number of CO rotation levels to default, 20 as set in taulines.h */
long int nCORotate = NCOROTATE;

/* say that space not allocated yet */
int lgCORotateMalloc = FALSE;

#include "abund.h"
#include "atoms.h"
#include "dynamics.h"
#include "assertresults.h"
#include "bit32.h"
#include "broke.h"
#include "ca.h"
#include "called.h"
#include "cddrive.h"
/* this will be set true when cdInit is called.  The definition is in cdInit.
 * Other routines will check that this is TRUE when they are called, 
 * to verify that cdInit was called first */
int lgcdInitCalled=FALSE;
#include "atmdat.h"
#include "colden.h"
#include "co.h"
#include "converge.h"
#include "continuum.h"
#include "coolheavy.h"
#include "magnetic.h"
#include "version.h"
#include "dense.h"
#include "doppvel.h"
#include "dump.h"
#include "embesq.h"
#include "extinc.h"
#include "fe.h"
#include "fudgec.h"
#include "geometry.h"
#include "grainvar.h"
#include "h2_priv.h"
#include "h2.h"
int lgH2_READ_DATA = FALSE;
#include "helike.h" 
#include "he.h"
#include "heavy.h"
#include "hextra.h"
#include "hmi.h"
#include "hydrogenic.h"

/* this is set true once space malloced, then never change
 * number of levels again with hydrogenic command, 
 * it is used in gethydro to never set number of levels one set here,
 * also to make sure MALLOC only happens one time  */
int lgHydroMalloc = FALSE;

/* these specify the form of the line redistribution function */
/* partial redistribution with wings */
const int ipPRD = 1;
/* complete redistribution, core only, no wings, Hummer's K2 function */
const int ipCRD = -1;
/* complete redistribution with wings */
const int ipCRDW = 2;
/* redistribution function for Lya, calls Hummer routine for H-like series only */
const int ipLY_A = -2;
/* core function for K2 destruction */
const int ipDEST_K2 = 1;
/* core function for complete redist destruction */
const int ipDEST_INCOM = 2;
/* core function for simple destruction */
const int ipDEST_SIMPL = 3;

/* indices within recombination coefficient array */
/* ipRecEsc is state specific escape probability*/
const int ipRecEsc = 2;
/* the net escaping, including destruction by background and optical deepth*/
const int ipRecNetEsc = 1;
/* ipRecRad is state specific radiative recombination rate*/
const int ipRecRad = 0;
#include "hyperfine.h"
#include "input.h"
#include "intalk.h"
#include "ionbal.h"
#include "iso.h"
#include "iterations.h"
#include "jeans_mass.h"
#include "kshllenr.h"
#include "lines.h"
/* these are the definitions of the line save arrays in lines.h */
LinDstSv *LineDSv=NULL;
LinSv *LineSv=NULL;
#include "map.h"
#include "mean.h"
#include "mewecoef.h"
#include "mole.h" 
#include "negdrg.h"
#include "neutrn.h"
#include "noexec.h"
#include "numderiv.h"
#include "oxy.h"
#include "optimize.h"
#include "peimbt.h"
#include "phycon.h"
#include "physconst.h"
#include "plot.h"
#include "sil.h"
#include "carb.h"
#include "atomfeii.h"
/* this is set true when space is allocated for the FeII arrays,
 * once this happens nFeIILevel cannot be changed with the atom feii levels command */
int lgFeIIMalloc=FALSE;
/* */
#include "pressure.h"
#include "prt.h"
#include "punch.h"
#include "radius.h"
#include "reason.h"
#include "rfield.h"
/* set true when malloced, init to false */
int lgRfieldMalloced=FALSE;
#include "opacity.h"
int lgOpacMalloced=FALSE;
#include "rt.h"
#include "secondaries.h"
#include "stopcalc.h"
#include "struc.h"
#include "thermal.h"
#include "timesc.h"
#include "trace.h"
#define NAMELEN 80
char chOptimFileName[NAMELEN] = "optimal.in";
#include "warnings.h"
#include "wind.h"
#include "yield.h"

void cdDefines(void)
{
	long int nelem;

	debug_fp = stderr;
	/* this will become the sqrt of COLL_SONST2 in cddefines */
	COLL_CONST = sqrt(COLL_CONST2);
	MILNE_CONST = sqrt(MILNE_CONST2);

	/* only reset number of levels if space has not been allocated - in all single run
	 * models this will be false, so levels will be set - exception is grid of models,
	 * this will only happen in first call */
	if( !lgHydroMalloc )
	{

		/* topoff for he-like sequence */
		/* TODO	2	do something with this or set to LONG_MAX? 
		* after moving collapsed levels to iso.nCollapsed, this 
		* variable does nothing.	
		iso.nTopOff[ipHE_LIKE] = 0;*/

		/* variables for H-like sequence */
		/* default nnumber of levels for hydrogen atom */
		for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
		{
			iso.n_HighestResolved[ipH_LIKE][nelem] = 15;
			/* >>chng 02 jan 18, increase by one to include 2s, so that all loops will be to < nLevels */
			iso.numLevels[ipH_LIKE][nelem] = iso.n_HighestResolved[ipH_LIKE][nelem] + 1;
			/* >>chng 03 jun 06, put number of topoff to one less than number,
			 * since intend to have more levels */
			iso.nTopOff[ipH_LIKE][nelem] = iso.numLevels[ipH_LIKE][nelem]-1;

			/* number of levels to print */
			iso.numPrintLevels[ipH_LIKE][nelem] = iso.numLevels[ipH_LIKE][nelem]-1;
			/* number of collapsed levels */
			iso.nCollapsed[ipH_LIKE][nelem] = 0;
		}
		for( nelem=ipHYDROGEN; nelem <= ipHELIUM; ++nelem )
		{
			/* >>chng 03 jun 01 up number of levels for H and He from 15 to 26 */
			iso.n_HighestResolved[ipH_LIKE][nelem] = 25;
			/* >>chng 02 jan 18, increase by one to include 2s, so that all loops will be to < nLevels */
			iso.numLevels[ipH_LIKE][nelem] = iso.n_HighestResolved[ipH_LIKE][nelem] + 1;

			/* lowest level is 2s,
			* in versions 90 this was equal to 10, but was 1 in 91
			* >>chng 99 jan 16, changed back to 10*/
			/*iso.nTopOff[ipH_LIKE] = 1;*/
			/* >>chng 03 jun 06, put number of topoff to one less than number,
			 * since intend to have more levels */
			/*iso.nTopOff[ipH_LIKE][nelem] = iso.numLevels[ipH_LIKE][nelem]-1;*/
			/* >>chng 03 oct 16, from above to 22 to avoid maser near topoff */
			iso.nTopOff[ipH_LIKE][nelem] = 22;

			/* number of levels to print */
			iso.numPrintLevels[ipH_LIKE][nelem] = iso.numLevels[ipH_LIKE][nelem]-1;
		}
#	if 0
#	endif

		/* total number of He-like levels, hydrogen does not exist, of course */
		iso.n_HighestResolved[ipHE_LIKE][ipHYDROGEN] = 0;
		iso.numLevels[ipHE_LIKE][ipHYDROGEN] = 0;
		iso.numPrintLevels[ipHE_LIKE][ipHYDROGEN] = 0;
		iso.nCollapsed[ipHE_LIKE][ipHYDROGEN] = 0;

		/* these are the total number of l levels, including the resolved 2 3P */
		for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
		{
			/* put at least three resolved and 1 collapsed in every element for he-like */
			iso.n_HighestResolved[ipHE_LIKE][nelem] = 3;
			iso.nCollapsed[ipHE_LIKE][nelem] = 1;
		}

		/* Make n=7 default max resolved for Helium itself */
		iso.n_HighestResolved[ipHE_LIKE][ipHELIUM] = 7;
		/* And n=5 for carbon, nitrogen, oxygen, and iron because they are quite abundant */
		iso.n_HighestResolved[ipHE_LIKE][ipCARBON] = 5;
		iso.n_HighestResolved[ipHE_LIKE][ipNITROGEN] = 5;
		iso.n_HighestResolved[ipHE_LIKE][ipOXYGEN] = 5;
		iso.n_HighestResolved[ipHE_LIKE][ipIRON] = 5;

		for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
		{
			/* at the very least do n=3 for every element that is turned on.	*/
			iso.n_HighestResolved[ipHE_LIKE][nelem] = MAX2( iso.n_HighestResolved[ipHE_LIKE][nelem], 3 );
			/* Calculate total levels for each element.	*/
			iso.numLevels[ipHE_LIKE][nelem] = ( ( 1 + iso.n_HighestResolved[ipHE_LIKE][nelem] ) * iso.n_HighestResolved[ipHE_LIKE][nelem] + 1 );
		}
	}

	/* number of electrons in valence shell that can compton recoil ionize */
	{
		long int nCom[LIMELM] = 
		{
			1 , 2 ,								/* K 1s shell */
			1 , 2 ,								/* L 2s shell */
			1 , 2 , 3 , 4 , 5 , 6 ,				/* L 2p shell */
			1 , 2 ,								/* M 3s shell */
			1 , 2 , 3 , 4 , 5 , 6 ,				/* M 3p shell */
			1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 ,		/* M 3d shell */
			1 , 2 ,								/* N 4s shell */
			1 , 2								/* N 4p shell */
		};

		for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
		{
			ionbal.nCompRecoilElec[nelem] = nCom[nelem];
		}
	}
}

#undef EXTERN
#define EXTERN extern
