/* 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 */
/* routines dealing with creation and evaluration of opacities */

/* set true when malloced, init to false */
extern int lgOpacMalloced;

/*OpacityCreateAll compute initial set of opacities for all species */
void OpacityCreateAll(void);

/*OpacityAdd1Subshell add opacity due to single species to main opacity array*/
void OpacityAdd1Subshell(
	/*ipOpac is opacity index within opac opacity offset for this species */
	long int ipOpac, 
	/* lower freq limit to opacity range on energy mesh */
	long int ipLowLim, 
	/* upper limit to opacity range on energy mesh */
	long int ipUpLim, 
	/* abundance, we bail if zero */
	float abundance,
	/* either static 's' or volitile 'v' */
	char chStat );

/*OpacityAddTotal derive total opacity for this position */
void OpacityAddTotal(void);

/*OpacityAdd1Element enter total photo cross section for all subshells 
 * of a single element into opacity array */
void OpacityAdd1Element(
		/* ipZ is 0 for H, 1 for He, etc */
		long int ipZ);

/* OpacityZero - called by OpacityAddTotal to zero opacity array after saving
 * it into oldopac */
void OpacityZero(void);

/* OpacityZeroOld - only set old opac to current value during search phase */
void OpacityZeroOld(void);

/*OpacityAdd1SubshellInduc add opacity of individual species, including stimulated emission */
void OpacityAdd1SubshellInduc(
  /* pointer to opacity offset with stack */
  long int ipOpac, 
  /* low energy limit to opacity bound */
  long int low, 
  /* high energy limit to opacity bound */
  long int ihi, 
  /* the abundance of the species in this level */
  double a, 
  /* the departure coefficient */
  double b,
  /* either 's' for static opacities, or 'v' for volitile */
  char chStat );

/* limit to number of opacity cells available in the opacity stack
#define	NOPSV	1300000L*/

/*lowest level where hydrogenic nu^-3 photo cross sect used*/
/*#define	NHPLPHOT	400*/

/* limit to number of cells in frequency array */
/* >>chng 00 oct 13, increase from 2000 to 2200 due to increase in FIR resolution */
/* >>chng 01 aug 19, increase from 2000 to 4000 due to increase in mid-range resolution */
/* >>chng 01 sep 24, get rid of NC ELL from code, totally MALLOC */
/*#define	NC ELL	4000L*/

EXTERN struct t_opacCom {

	/* [0] is optical depth for current position,
	 * [1] is total optical depth from previous iteration */

	/*TauAbsGeo, TauScatGeo, absorption, scattering, optical depths,
	 * at current position, includes through far side of slab*/
	float **TauAbsGeo;
	float **TauScatGeo;

	/* TauTotalGeo is total optical depth at each energy, at current position, 
	 * includes through far side of slab */
	float **TauTotalGeo;

	/* these are the integrated absorption and scattering optical depths
	 * to the illuminated face of the cloud */
	float *TauAbsFace, 
	  *TauScatFace;

	/*local absorption opacity, cm-1*/
	double *opacity_abs;
	 
	/* local scattering opacity, cm-1 */
	double *opacity_sct ;

	/*save previous opacity */
	double *OldOpacSave;

	/* albedo is local gas albedo*/
	double *albedo;

	/* initial opacities from zone 1, used to reset opacity at restart */
	/* and the saved value of local opacity*/
	double *opacity_abs_savzon1;
	/* opacity_sct_savzon1 save local opacity at start of calculations*/
	double *opacity_sct_savzon1;

	/* these static opacities are only evaluated 
	 * one time per zone*/
	double *OpacStatic;

	/* density/temp factors needed for free-free opacity */
	double *FreeFreeOpacity;

	/* exp(-dTau) for the current zone, evaluated in radinc */
	double *ExpZone;

	/* saves value of E2(tau), exponential integral,
	 * where tau is optical depth to illuminated face */
	float *e2TauAbs;

	/* exp(-tau) to illunated face */
	float *ExpmTau;

	/* factors that account for attenuation of light across this zone
	 * should be nearly unity */
	float *tmn;

	/* this flag is set true in zero, and set false with the
	 * no static opacities command.  When false always update
	 * all opacites */
	int lgOpacStatic;

	/* this flag is set true in OpacityZero 
	 * when the OpacStatic array is zeroed, 
	 * and is false if the array has been left alone. 
	 * all later opacities must be reevaluated when this
	 * flag is true */
	int lgRedoStatic;

	  /* this is the stack used to hold opacities - entered one time when code
	 * is initialized, in routine OpacityCreateAll */
	/* >>chng 01 mar 08, from float to double */
	double *OpacStack/*[NOPSV]*/;

	/* parameters for CaII pumping */
	float tpcah[2];

	/* taumin is the smallest optical depths allowed, */
	float taumin;

	/* tlamin is smallest Lya optical depth, 
	 * modified to large value if case b used */
	float tlamin;

	/* this flag says we are doing case b */
	int lgCaseB;

	/* this flag turns off n=2 collisional excitations for comparison 
	 * with their paper, keyword HUMMER on caseb command */
	int lgCaseB_HummerStorey;

	/* this flag turns off excited state photoionization, 
	 * keyword NO PHOTO on caseb command */
	int lgCaseB_no_photo;

	/* another case b option, turn off background opacities, no Pdest */
	int lgCaseB_no_pdest;

	/* optical depths to compton and H- */
	float telec, 
	  thmin;

	/* flag set if negative opacities occured */
	int lgOpacNeg;

	/* flag saying whether (true) or not (false) scattering opacity
	 * is enabled */
	int lgScatON;

	/* IO unit to punch negative opacities */
	int lgNegOpacIO;

	/* variable dealing with the option to generate and use an
	 * ancillary file containing a stored opacity set.  These
	 * are all initialized in zero1  */
	/*lgCompileOpac flag saying to compile opacity only */
	int lgCompileOpac;

	/*lgOpacExist flag saying that existing opacities used*/
	int lgOpacExist;

	/* set FALSE if no file opacity command entered, option to
	 * ignore opacity file */
	int lgUseFileOpac ;

	/* indicates whether (true) or not the outward optical
	 * depths have been predicted yet */
	int lgTauOutOn;

	/* which opacity to punch out  */
	char chOpcTyp[5];

	/* total number of opacity cells used in opacity stack
	 * in OpacityCreateAll used as a counter to remember where
	 * next opacity goes */
	long int nOpacTot;

	/*
	 *
	 * NBNBNBthis must exactly parallel the read/write statements
	 * in OpacityCreateAll
	 */

	/* ipRayScat opacity pointer for rayleigh scattering*/
	long int ipRayScat, 

	/* iopcom compton scatterin, total recoil*/
	iopcom, 

	/* ippr is threshold for pair production, ioppr is opacity offset*/
	ippr, 
	ioppr, 

	/* ipBrems opacity offset pointer for brems (free-free)*/
	ipBrems, 

	/* iphmra ratio of h minus to neut h ff brems opacity*/
	iphmra, 

	/* iphmop H- bound free opacity*/
	iphmop, 

	/* ih2pnt lower, upper limits to bound, ih2pof, opacity offset*/
	ih2pnt[2], 
	ih2pof, 

	/* iophe1 points for photo from singlet levels*/
	iophe1[9], 

	/* lowest levels of he triplets*/
	ioptri, 

	/*  
	 * ipElement[nelem][ioniz][shell][purpose],
	 * array index on fortran scale, to energy range are set in ipShells, called by ContCreatePointers
	 * pointers to opacity stack offset defined in OpacityCreate1Element where
	 * opacities are set
	 *
	 * first dim 
	 * [0] is pointer to low-energy threshold in energy array
	 * [1] is highest energy for shell, set by LimitSh
	 * [2] is opacity offset within large opacity stack
	 *
	 * if shell does not exist, set upper limit
	 * to less than lower limit so this never looped upon
	 * these are used as flags by LimitSh to check whether
	 * this is a real shell - 
	 * following code is in ipShells for this non-existant case 
	 * OpacPoint.ipElement[nelem][ion][nshell][0] = 2;
	 * OpacPoint.ipElement[nelem][ion][nshell][1] = 1;
	 * all routines must protect against this since opacities for these
	 * undefined shells is also undefined
	 *
	 * second dim is shell, 0 for k shell, up to 6,
	 * (Shell.chShell[ns] give the label for each of these shells)
	 * These quantities are only defined for the number of shells that exist. 
	 * The number of shells is given by Heavy.nsShells[nelem][ion]
	 *
	 * third dimension is ion stage, 0 for atom
	 * 
	 * last dim is atomic number of element, 0 for H
	 *
	 */
	ipElement[LIMELM][LIMELM][7][3], 

	/* in1 is [NI] excited state*/
	in1[3], 

	/* pointers to oxy excited states*/
	ipo3exc[3], 
	ipo3exc3[3], 
	ipo1exc[3], 

	/* photo to excited O+ levels*/
	iopo2d, 

	/* photoionization from upper lev of Mg II 2798*/
	ipmgex, 
	ipOpMgEx, 

	/* calcium excited states*/
	ica2ex[2], 
	ica2op;

	/* threshold for photoionization of H2 - 15.4 eV */
	long int ipH2_photo_thresh , ipH2_photo_opac_offset ;

	/* index of carbon k-shell threshold in energy array*/
	long int ipCKshell;

	/* smallest ots rate, set to 0 in scalar */
	float otsmin;

	/* max correction for stim emission in continuum opacities at Lyman and Balmer edges */
	float stimax[2];

	}	opac;


