/* 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 */
/*OpacityCreateAll compute initial set of opacities for all species */
/*OpacityCreate1Element generate ionic subshell opacities by calling atmdat_phfit */
/*opacity_more_memory allocate more memory for opacity stack */
/*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
/*hmiopc derive total H- H minus opacity */
/*rayleh compute Rayleigh scattering cross section for Lya */
/*OpacityValenceRescale routine to rescale non-OP valence shell cross sections */
/*Yan_H2_CS - cross sections for the photoionization of H2 */
/******************************************************************************
 *NB NB NB  NB NB NB NB NB NB NB  NB NB NB NB
 * everything set here must be written to the opacity store files
 *
 ****************************************************************************** */
#include "cddefines.h"
#define	NCSH2P	10
#include "physconst.h"
#include "dense.h"
#include "iso.h"
#include "hydrogenic.h"
#include "oxy.h"
#include "trace.h"
#include "heavy.h"
#include "kshllenr.h"
#include "rfield.h"
#include "hmi.h"
#include "path.h"
#include "atmdat.h"
#include "punch.h"
#include "grains.h"
#include "interpolate.h"
#include "hydro_bauman.h"
#include "opacity.h"
#include "helike_recom.h"

/* number of char in magic number */
#define NMAGIC 7

/* the magic number itself, the current date (YYMMDD) */
#define MAGIC "010118"

/* limit to number of opacity cells available in the opacity stack*/
static long int ndimOpacityStack = 1700000L;

/*Yan_H2_CS - cross sections for the photoionization of H2 
 * may be represented analytically in the paper 
 *>>refer	H2	photo cs	Yan, M.; Sadeghpour, H. R.; Dalgarno, A. 1998, ApJ, 496, 1044
 * This is revised version given in ERRATUM 2001, ApJ, 559, 1194 
 * return value is photo cs in cm-2 */
static double Yan_H2_CS( double energy_ryd /* photon energy in ryd */)
{

	double energy_keV; /* keV */
	double cross_section; /* barns */
	double x; /* x = E/15.4 */
	double xsqrt , x15 , x2;
	double energy = energy_ryd * EVRYD;

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

	/* energy relative to threshold */
	x = energy / 15.4;
	energy_keV = energy/1000.0;
	xsqrt = sqrt(x);
	x15 = x * xsqrt;
	x2 = x*x;

	if( energy < 15.4 )
	{
		cross_section = 0.;
	}

	else if (energy >= 15.4 && energy < 18.)
	{
		cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2);
		/* this equation is obviously negative for x = 1 */
		cross_section = MAX2( 0. , cross_section );
	}

	else if (energy >= 18. && energy <= 30.)
	{
		cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
	}

	else if (energy > 30. && energy <= 85.)
	{
		cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
	}

	/* the high-energy tail */
	else 
	{
		cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*x2) - 81.885/(x*x2))/pow(energy_keV,3.5);
	}

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

	return( cross_section * 1e-24 );
}

/*OpacityCreate1Element generate opacities for entire element by calling atmdat_phfit */
static void OpacityCreate1Element(long int nelem);

/*opacity_more_memory allocate more memory for opacity stack */
static void opacity_more_memory(void);

/*hmiopc derive total H- H minus opacity */
static double hmiopc(double freq);

/*rayleh compute Rayleigh scattering cross section for Lya */
static double rayleh(double ener);

/*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
static double Opacity_iso_photo_cs( float energy , long ipISO , long nelem , long n );

/*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */
static void OpacityCreateReilMan(long int low, 
  long int ihi, 
  float cross[], 
  long int ncross, 
  long int *ipop, 
  const char *chLabl);

static int lgRealloc = FALSE;

/*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */
static void OpacityCreatePowerLaw(
	/* lower energy limit on continuum mesh */
	long int ilo, 
	/* upper energy limit on continuum mesh */
	long int ihi, 
	/* threshold cross section */
	double cross, 
	/* power law index */
	double s, 
	/* pointer to opacity offset where this starts */
	long int *ip);

/*ofit compute cross sections for all shells of atomic oxygen */
double ofit(double e, 
	  float opart[]);

/*OpacityValenceRescale routine to rescale non-OP valence shell cross sections for atom */
static void OpacityValenceRescale(
	/* element number on C scale */
	long int nelem ,
	/* scale factor, must be >= 0. */
	double scale )
{

	long int ion , nshell , low , ihi , ipop , ip;

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

	/* return if element is not turned on 
	 * >>chng 05 oct 19, this had not been done, so low in the opacity offset below was
	 * not set, and opacity index was negative - only problem when K turned off */
	if( !dense.lgElmtOn[nelem] )
	{
#		ifdef DEBUG_FUN
		fputs( " <->OpacityValenceRescale()\n", debug_fp );
#		endif
		return;
	}

	/* scale better be >= 0. */
	ASSERT( scale >= 0. );

	ion = 0;
	/* this is valence shell on C scale */
	nshell = Heavy.nsShells[nelem][ion] - 1;

	/* set lower and upper limits to this range */
	low = opac.ipElement[nelem][ion][nshell][0];
	ihi = opac.ipElement[nelem][ion][nshell][1];
	ipop = opac.ipElement[nelem][ion][nshell][2];

	/* loop over energy range of this shell */
	for( ip=low-1; ip < ihi; ip++ )
	{
		opac.OpacStack[ip-low+ipop] *= scale;
	}

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

void OpacityCreateAll(void)
{
	long int i, 
	  ipISO ,
	  n ,
	  need ,
	  nelem;

	float opart[7];

	double crs, 
	  dx,
	  eps, 
	  thom, 
	  thres, 
	  x;

	/* this will be the file name we will try to open for the opacities */
	char chFileName[FILENAME_PATH_LENGTH_2];

	FILE *ioOPAC; /* will be pointer to opacity.opc data file */

	/* >>chng 02 may 29, change to lgOpacMalloced */
	/* remember whether opacities have ever been evaluated in this coreload
	static int lgOpEval = FALSE; */

	/* fits to cross section for photo dist of H_2^+ */
	static float csh2p[NCSH2P]={6.75f,0.24f,8.68f,2.5f,10.54f,7.1f,12.46f,
	  6.0f,14.28f,2.7f};

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

	/* H2+ h2plus h2+ */

	/* make and print dust opacities
	 * fill in dstab and dstsc, totals, zero if no dust
	 * may be different if different grains are turned on */
	InitGrains();

	/* flag lgOpacMalloced says whether opacity stack has been generated
	 * only do this one time per core load  */
	/* >>chng 02 may 29, from lgOpEval to lgOpacMalloced */
	if( lgOpacMalloced )
	{
		/* this is not the first time code called */
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " OpacityCreateAll called but NOT evaluated since already done.\n" );
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->OpacityCreateAll()\n", debug_fp );
#		endif
		return;
	}

	lgOpacMalloced = TRUE;

	/* create the space for the opacity stack */
	if( (opac.OpacStack = (double*)MALLOC((size_t)ndimOpacityStack*sizeof(double)) ) == NULL )
		BadMalloc();

	/* should we even try to bring in opacity file? */
	if( opac.lgUseFileOpac )
	{
		/* option to compile opacities into file for later use */
		/* first try local directory */
		strcpy( chFileName , "opacity.opc" );
		ioOPAC = fopen( chFileName , "rb" );
		if( ioOPAC == NULL )
		{
			/* did not get it here, lets try the path if it was set */
			if( lgDataPathSet )
			{
				/* first stuff the path name itself */
				strcpy( chFileName , chDataPath );
				strcat( chFileName , "opacity.opc" );
				/* now try to open this file */
				ioOPAC = fopen( chFileName , "rb" );
				if( ioOPAC == NULL )
				{
					/* still not here, generate new opacities */
					opac.lgOpacExist = FALSE;
				}
				else
				{
					/* found it, will use this file */
					opac.lgOpacExist = TRUE;
				}
			}
			else
			{
				opac.lgOpacExist = FALSE;
			}
		}
		else
		{
			opac.lgOpacExist = TRUE;
		}

		/*	opac.lgCompileOpac is set TRUE with compile opacities command */
		if( !opac.lgCompileOpac && opac.lgOpacExist && (ioOPAC!=NULL) )
		{
			char chMagic[NMAGIC];
			/* opacity file exists and have not been asked to compile it */

			/* first get the magic number */
			n = (long)fread( chMagic , NMAGIC , sizeof(char), ioOPAC );
			if( strcmp( chMagic , MAGIC ) != 0 )
			{
				fprintf( ioQQQ," The starting magic number in the compiled opacity file is incorrect.\n");
				fprintf( ioQQQ," Please recompile opacities, or delete opacity.opc in the data directory.\n");
				fprintf( ioQQQ," The magic number was %s\n", MAGIC );
				puts( "[Stop in OpacityCreateAll]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* next get the big stack of opacity data */
			n = (long)fread( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOPAC );
			if( ((unsigned)n-sizeof(opac.OpacStack)) != 0 )
			{
				fprintf( ioQQQ, " problem trying to read opacity.opc\n" );
				fprintf( ioQQQ, " I expected to read %li words, but fread returned only %li\n",
					(long)sizeof(opac.OpacStack),n);
				fprintf( ioQQQ, " Try recompiling the opacities with the COMPILE OPACITY command,\n" );
				fprintf( ioQQQ, " or tell the code not to use the compiled opacities with the NO FILE OPACITY command.\n" );

				/* did we hit end of file? */
				if( feof(ioOPAC ) )
				{
					fprintf( ioQQQ, " end of file hit\n" );
				}
				/* some other error condition? */
				if( ferror(ioOPAC ) )
				{
					fprintf( ioQQQ, " end of file hit\n" );
				}
				puts( "[Stop in OpacityCreateAll]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* NG following must exactly mirror contents of common opac */
			/* we will add unwritten bytes to this, then check if still zero at end */
			n = 0;
			/* abs cannot be used here since it takes an int, not a long */
			n -= (long)(fread( iso.ipOpac[ipH_LIKE],	1,sizeof(iso.ipOpac[ipH_LIKE]),
				ioOPAC ) - sizeof(iso.ipOpac[ipH_LIKE]));

			n -= (long)(fread( &opac.ipRayScat,1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.iopcom,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.ippr,		1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.ioppr,		1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.ipBrems,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.iphmra,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.iphmop,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( opac.ih2pnt,		1,sizeof(opac.ih2pnt),
				ioOPAC ) - sizeof(opac.ih2pnt));

			n -= (long)(fread( &opac.ih2pof,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( opac.iophe1,		1,sizeof(opac.iophe1),
				ioOPAC ) - sizeof(opac.iophe1));

			n -= (long)(fread( &opac.ioptri,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( opac.ipElement,	1,sizeof(opac.ipElement),
				ioOPAC ) - sizeof(opac.ipElement));

			n -= (long)(fread( opac.in1,			1,sizeof(opac.in1),
				ioOPAC ) - sizeof(opac.in1));

			n -= (long)(fread( opac.ipo3exc,	1,sizeof(opac.ipo3exc),
				ioOPAC ) - sizeof(opac.ipo3exc));

			n -= (long)(fread( opac.ipo3exc3,	1,sizeof(opac.ipo3exc3),
				ioOPAC ) - sizeof(opac.ipo3exc3));

			n -= (long)(fread( opac.ipo1exc,	1,sizeof(opac.ipo1exc),
				ioOPAC ) - sizeof(opac.ipo1exc));

			n -= (long)(fread( &opac.iopo2d,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.ipmgex,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( &opac.ipOpMgEx,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			n -= (long)(fread( opac.ica2ex,		1,sizeof(opac.ica2ex),
				ioOPAC ) - sizeof(opac.ica2ex));

			n -= (long)(fread( &opac.ica2op,	1,sizeof(long) ,
				ioOPAC ) - sizeof(long));

			/* following contains AnuOrg and will be checked against it as a sanity check*/
			n -= (long)(fread( opac.tmn,		1,sizeof(opac.tmn),
				ioOPAC ) - sizeof(opac.tmn));

			/* now check whether n is still zero, problem if not */
			if( n!= 0 )
			{
				fprintf( ioQQQ, " problem trying to read opacity.opc pointers\n" );
				fprintf( ioQQQ, " fread was short by %li\n",
					n);
				fprintf( ioQQQ, " Try recompiling the opacities with the COMPILE OPACITY command,\n" );
				fprintf( ioQQQ, " or tell the code not to use the compiled opacities with the NO FILE OPACITY command.\n" );
				puts( "[Stop in OpacityCreateAll]" );
				cdEXIT(EXIT_FAILURE);
			}

			/*lint -e777 test floats equality */
			if( opac.tmn[rfield.nupper-1] != rfield.AnuOrg[rfield.nupper-1] )
			/*lint +e777 test floats equality */
			{
				fprintf( ioQQQ, " Energy grid in opacity file opacity.opc opacity.pnt do not agree with this version of the code - please recompile.\n" );
				fprintf( ioQQQ, " Values are:%12.4e%12.4e\n", 
				  rfield.AnuOrg[rfield.nupper-1], opac.tmn[rfield.nupper-1] );
				puts( "[Stop in OpacityCreateAll]" );
				cdEXIT(EXIT_FAILURE);
			}

			if( trace.lgTrace )
			{
				fprintf( ioQQQ, " OpacityCreateAll called, values from file =%s=\n",chFileName );
			}

			/* note that opacities now exist */
			opac.lgOpacExist = TRUE;

			/* zero out opac since this array sometimes addressed before OpacityAddTotal called */
			for( i=0; i < rfield.nupper; i++ )
			{
				opac.opacity_abs[i] = 0.;
			}
			
#			ifdef DEBUG_FUN
			fputs( " <->OpacityCreateAll()\n", debug_fp );
#			endif
			return;
		}
		else
		{
			opac.lgOpacExist = FALSE;
		}
	}
	else 
	{
		/* this branch if NO FILE OPACITY command entered */
		opac.lgOpacExist = FALSE;
	}

	/* initialize the arrays of hydro photo fitting coefficients */
	atmdat_H_phot_cs();

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " OpacityCreateAll called, evaluating.\n" );
	}

	/* zero out opac since this array sometimes addressed before OpacityAddTotal called */
	for( i=0; i < rfield.nupper; i++ )
	{
		opac.opacity_abs[i] = 0.;
	}

	/* nOpacTot is number of opacity cells in OpacStack filled so far by opacity generting routines */
	opac.nOpacTot = 0;

	/* photoionization of h, he-like iso electronic sequences */
	for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.lgElmtOn[nelem] )
			{
				long int nupper;

				/* this is the opacity offset in the general purpose pointer array */
				/* indices are type, shell. ion, element, so this is the inner shell,
				 * NB - this works for H and He, but the second index should be 1 for Li */
				opac.ipElement[nelem][nelem-ipISO][0][2] = opac.nOpacTot + 1;

				/* gound state goes to high-energy limit of code, 
				 * but excited states only go up to edge for ground */
				nupper = rfield.nupper;
				for( n=0; n < iso.numLevels_max[ipISO][nelem]; n++ )
				{
					/* this is array index to the opacity offset */
					iso.ipOpac[ipISO][nelem][n] = opac.nOpacTot + 1;

					/* first make sure that first energy point is at least near the limit */
					/* >>chng 01 sep 23, increased factor form 0.98 to 0.94, needed since cells now go
					 * so far into radio, where resolution is poor */
					ASSERT( rfield.AnuOrg[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0.94f * 
						iso.xIsoLevNIonRyd[ipISO][nelem][n] );

					/* number of cells we will need to do this level */
					need = nupper - iso.ipIsoLevNIonCon[ipISO][nelem][n] + 1;
					ASSERT( need > 0 );

					if( opac.nOpacTot + need > ndimOpacityStack )
						opacity_more_memory();

					for( i=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; i < nupper; i++ )
					{
						opac.OpacStack[i-iso.ipIsoLevNIonCon[ipISO][nelem][n]+iso.ipOpac[ipISO][nelem][n]] = 
							Opacity_iso_photo_cs( rfield.AnuOrg[i] , ipISO , nelem , n );
					}

					opac.nOpacTot += need;
					/* for all excited levels high-energy limit is edge for ground state */
					nupper = iso.ipIsoLevNIonCon[ipISO][nelem][0];
				}
			}
		}
	}

	/* This check will get us through Klein-Nishina below.	*/
	/* >>chng 02 may 08, by Ryan.  Added this and other checks for allotted memory.	*/
	if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] + rfield.nupper > ndimOpacityStack )
		opacity_more_memory();

	/* Lyman alpha damping wings - Rayleigh scattering */
	opac.ipRayScat = opac.nOpacTot + 1;
	for( i=0; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
	{
		opac.OpacStack[i-1+opac.ipRayScat] = rayleh(rfield.AnuOrg[i]);
	}
	opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];

	/* ==============================================================
	 * this block of code defines the electron scattering cross section
	 * for all energies */

	/* assume Thomson scattering up to ipCKshell, 20.6 Ryd=0.3 keV */
	thom = 6.65e-25;
	opac.iopcom = opac.nOpacTot + 1;
	for( i=0; i < opac.ipCKshell; i++ )
	{
		opac.OpacStack[i-1+opac.iopcom] = thom;
		/*fprintf(ioQQQ,"%.3e\t%.3e\n", 
			rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/
	}

	/* Klein-Nishina from eqn 7.5, 
	 * >>refer	Klein-Nishina	cs	Rybicki and Lightman */
	for( i=opac.ipCKshell; i < rfield.nupper; i++ )
	{
		dx = rfield.AnuOrg[i]/3.7573e4;

		opac.OpacStack[i-1+opac.iopcom] = thom*3.e0/4.e0*((1.e0 + 
		  dx)/POW3(dx)*(2.e0*dx*(1.e0 + dx)/(1.e0 + 2.e0*dx) - log(1.e0+
		  2.e0*dx)) + 1.e0/2.e0/dx*log(1.e0+2.e0*dx) - (1.e0 + 3.e0*
		  dx)/POW3(1.e0 + 2.e0*dx));
		/*fprintf(ioQQQ,"%.3e\t%.3e\n", 
			rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-1+opac.iopcom] );*/
	}
	opac.nOpacTot += rfield.nupper - 1 + 1;

	/* ============================================================== */

	/* This check will get us through "H- hminus H minus bound-free opacity" below.	*/
	/* >>chng 02 may 08, by Ryan.  Added this and other checks for allotted memory.	*/
	if( opac.nOpacTot + 3*rfield.nupper - opac.ippr + iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 2 > ndimOpacityStack )
		opacity_more_memory();

	/* pair production */
	opac.ioppr = opac.nOpacTot + 1;
	for( i=opac.ippr-1; i < rfield.nupper; i++ )
	{
		/* pair production heating rate for unscreened H + He
		 * fit to figure 41 of Experimental Nuclear Physics,
		 * Vol1, E.Segre, ed */

		x = rfield.AnuOrg[i]/7.512e4*2.;

		opac.OpacStack[i-opac.ippr+opac.ioppr] = 5.793e-28*
		  POW2((-0.46737118 + x*(0.349255416 + x*0.002179893))/(1. + 
		  x*(0.130471301 + x*0.000524906)));
		/*fprintf(ioQQQ,"%.3e\t%.3e\n", 
			rfield.AnuOrg[i]*EVRYD , opac.OpacStack[i-opac.ippr+opac.ioppr] );*/
	}
	opac.nOpacTot += rfield.nupper - opac.ippr + 1;

	/* brems (free-free) opacity */
	opac.ipBrems = opac.nOpacTot + 1;

	for( i=0; i < rfield.nupper; i++ )
	{
		/* missing factor of 1E-20 to avoid underflow
		 * free free opacity needs g(ff)*(1-exp(hn/kT))/SQRT(T)*1E-20 */
		opac.OpacStack[i-1+opac.ipBrems] = 
			/*(float)(1.03680e-18/POW3(rfield.AnuOrg[i]));*/
			/* >>chng 00 jun 05, descale by 1e10 so that underflow at high-energy
			 * end does not occur */
			1.03680e-8/POW3(rfield.AnuOrg[i]);
	}
	opac.nOpacTot += rfield.nupper - 1 + 1;

	opac.iphmra = opac.nOpacTot + 1;
	for( i=0; i < rfield.nupper; i++ )
	{
		/* following is ratio of h minus to neut h bremss opacity */
		opac.OpacStack[i-1+opac.iphmra] = 0.1175*rfield.anusqr[i];
	}
	opac.nOpacTot += rfield.nupper - 1 + 1;

	opac.iphmop = opac.nOpacTot + 1;
	for( i=hmi.iphmin-1; i < iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]; i++ )
	{
		/* H- hminus H minus bound-free opacity */
		opac.OpacStack[i-hmi.iphmin+opac.iphmop] = 
			hmiopc(rfield.AnuOrg[i]);
	}
	opac.nOpacTot += iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - hmi.iphmin + 1;

	opac.ipH2_photo_opac_offset = opac.nOpacTot + 1;
	for( i=opac.ipH2_photo_thresh-1; i < rfield.nupper; i++ )
	{
		/* >>chng 05 nov 24, add H2 photoionization cross section */
		opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] = 
			Yan_H2_CS(rfield.AnuOrg[i]);
		/*fprintf(ioQQQ,"DEBUG H2 photo\t%.3e\t%.3e\t%.3e\n",
			rfield.AnuOrg[i],
			opac.OpacStack[i-opac.ipH2_photo_thresh + opac.ipH2_photo_opac_offset] ,
			Opacity_iso_photo_cs( rfield.AnuOrg[i] , 0 , 0 , 0 )*2. );*/
	}
	opac.nOpacTot += rfield.nupper - opac.ipH2_photo_thresh + 1;

	/* H2+ H2P h2plus photoabsorption
	 * cs from 
	 * >>refer	H2+	photodissoc	Buckingham, R.A., Reid, S., & Spence, R. 1952, MNRAS 112, 382, 0 K temp */
	OpacityCreateReilMan(opac.ih2pnt[0],opac.ih2pnt[1],csh2p,NCSH2P,&opac.ih2pof,
	  "H2+ ");

	/* This check will get us through "HeI singlets neutral helium ground" below.	*/
	/* >>chng 02 may 08, by Ryan.  Added this and other checks for allotted memory.	*/
	if( opac.nOpacTot + rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1 > ndimOpacityStack )
		opacity_more_memory();

	/* HeI singlets neutral helium ground */
	opac.iophe1[0] = opac.nOpacTot + 1;
	opac.ipElement[ipHELIUM][0][0][2] = opac.iophe1[0];
	for( i=iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1; i < rfield.nupper; i++ )
	{
		crs = atmdat_phfit(2,2,1,rfield.AnuOrg[i]*EVRYD,PhFitOn.lgPhFit);
		opac.OpacStack[i-iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]+opac.iophe1[0]] = 
			crs*1e-18;
	}
	opac.nOpacTot += rfield.nupper - iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] + 1;

	/* these are opacity offset points that would be defined in OpacityCreate1Element,
	 * but this routine will not be called for H and He
	 * generate all heavy element opacities, everything heavier than He,
	 * nelem is on the C scale, so Li is 2 */
	/*>>chng 99 jan 27, do not reevaluate hydrogenic opacity below */
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		if( dense.lgElmtOn[nelem] )
		{
			OpacityCreate1Element(nelem);
		}
	}

	/* option to rescale atoms of some elements that were not done by opacity project
	 * the valence shell - two arguments - element number on C scale, and scale factor */
	/*>>chng 05 sep 26, fudge factor to get atomic K fraction along well defined line of sight
	 * to be observed value - this is ratio of cross sections, actual value is very uncertain since
	 * differences betweeen Verner & opacity project are huge */
	OpacityValenceRescale( ipPOTASSIUM , 5. );

	/* now add on some special cases, where exicted states, etc, come in */
	/* Nitrogen
	 * >>refer	n1	photo	Henry, R., ApJ 161, 1153.
	 * photoionization of excited level of N+ */
	OpacityCreatePowerLaw(opac.in1[0],opac.in1[1],9e-18,1.75,&opac.in1[2]);

	/* atomic Oxygen
	 * only do this if 1996 Verner results are used */
	if( dense.lgElmtOn[ipOXYGEN] && !PhFitOn.lgPhFit )
	{
		/* This check will get us through this loop.	*/
		/* >>chng 02 may 08, by Ryan.  Added this and other checks for allotted memory.	*/
		if( opac.nOpacTot + opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1 > ndimOpacityStack )
			opacity_more_memory();

		/* integrate over energy range of the valence shell of atomic oxygen*/
		for( i=opac.ipElement[ipOXYGEN][0][2][0]-1; i < opac.ipElement[ipOXYGEN][0][2][1]; i++ )
		{
			/* call special routine to evaluate partial cross section for OI shells */
			eps = rfield.AnuOrg[i]*EVRYD;
			crs = ofit(eps,opart);

			/* this will be total cs of all processes leaving shell 3 */
			crs = opart[0];
			for( n=1; n < 6; n++ )
			{
				/* add up table of cross sections */
				crs += opart[n];
			}
			/* convert to cgs and overwrite cross sections set by OpacityCreate1Element */
			crs *= 1e-18;
			opac.OpacStack[i-opac.ipElement[ipOXYGEN][0][2][0]+opac.ipElement[ipOXYGEN][0][2][2]] = crs;
		}
		/* >>chng 02 may 09 - this was a significant error */
		/* >>chng 02 may 08, by Ryan.  This loop did not update total slots filled.	*/
		opac.nOpacTot += opac.ipElement[ipOXYGEN][0][2][1] - opac.ipElement[ipOXYGEN][0][2][0] + 1;
	}

	/* Henry nubmers for 1S excit state of OI, OP data very sparse */
	OpacityCreatePowerLaw(opac.ipo1exc[0],opac.ipo1exc[1],4.64e-18,0.,&opac.ipo1exc[2]);

	/* photoionization of excited level of O2+ 1D making 5007
	 * fit to TopBase Opacity Project cs */
	OpacityCreatePowerLaw(opac.ipo3exc[0],opac.ipo3exc[1],3.8e-18,0.,&opac.ipo3exc[2]);

	/* photoionization of excited level of O2+ 1S making 4363 */
	OpacityCreatePowerLaw(opac.ipo3exc3[0],opac.ipo3exc3[1],5.5e-18,0.01,
	  &opac.ipo3exc3[2]);

	/* This check will get us through the next two steps.	*/
	/* >>chng 02 may 08, by Ryan.  Added this and other checks for allotted memory.	*/
	if( opac.nOpacTot + iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1 
		+ iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1 > ndimOpacityStack )
		opacity_more_memory();

	/* photoionization to excited states of O+ */
	opac.iopo2d = opac.nOpacTot + 1;
	thres = rfield.AnuOrg[oxy.i2d-1];
	for( i=oxy.i2d-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]; i++ )
	{
		crs = 3.85e-18*(4.4*pow(rfield.AnuOrg[i]/thres,-1.5) - 3.38*
		  pow(rfield.AnuOrg[i]/thres,-2.5));

		opac.OpacStack[i-oxy.i2d+opac.iopo2d] = crs;
	}
	opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s] - oxy.i2d + 1;

	/* magnesium
	 * photoionization of excited level of Mg+
	 * fit to opacity project data Dima got */
	opac.ipOpMgEx = opac.nOpacTot + 1;
	for( i=opac.ipmgex-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
	{
		opac.OpacStack[i-opac.ipmgex+opac.ipOpMgEx] = 
			(0.2602325880970085 + 
		  445.8558249365131*exp(-rfield.AnuOrg[i]/0.1009243952792674))*
		  1e-18;
	}
	opac.nOpacTot += iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - opac.ipmgex + 1;

	ASSERT( opac.nOpacTot < ndimOpacityStack );

	/* Calcium
	 * excited states of Ca+ */
	OpacityCreatePowerLaw(opac.ica2ex[0],opac.ica2ex[1],4e-18,1.,&opac.ica2op);

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, 
			" OpacityCreateAll return OK, number of opacity cells used in OPSC= %ld and OPSV is dimensioned %ld\n", 
		  opac.nOpacTot, ndimOpacityStack );
	}

	/* option to compile opacities into file for later use 
	 * this is executed if the 'compile opacities' command is entered */
	if( opac.lgCompileOpac )
	{
		ioOPAC = fopen("opacity.opc","wb");
		if( ioOPAC == NULL )
		{
			fprintf( ioQQQ, " problem trying to open opacity.opc\n" );
			puts( "[Stop in OpacityCreateAll]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* first get the magic number */
		n = (long)fwrite( MAGIC , NMAGIC , sizeof(char), ioOPAC );

		/* write the data, remembering how many bytes were dumped */
		n = (long)fwrite( opac.OpacStack , 1, sizeof(opac.OpacStack), ioOPAC );
		if( ((size_t)n-sizeof(opac.OpacStack)) != 0 )
		{
			fprintf( ioQQQ, " problem trying to write opacity.opc\n" );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite returned only %li\n",
				(long)sizeof(opac.OpacStack),n);
			puts( "[Stop in OpacityCreateAll]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* NG following must exactly mirror contents of common opac */
		/* we will add unwritte bytes to this, then check if still zero at end */
		n = 0;
		/* abs cannot be used here since it takes an int, not a long */
		n -= (long)( fwrite( iso.ipOpac[ipH_LIKE], 1,sizeof(iso.ipOpac[ipH_LIKE]),
			ioOPAC ) - sizeof(iso.ipOpac[ipH_LIKE]));

		n -= (long)( fwrite( &opac.ipRayScat, 1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.iopcom, 1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.ippr, 1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.ioppr,		1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.ipBrems,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.iphmra,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.iphmop,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( opac.ih2pnt,		1,sizeof(opac.ih2pnt),
			ioOPAC ) - sizeof(opac.ih2pnt));

		n -= (long)( fwrite( &opac.ih2pof,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( opac.iophe1,		1,sizeof(opac.iophe1),
			ioOPAC ) - sizeof(opac.iophe1));

		n -= (long)( fwrite( &opac.ioptri,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( opac.ipElement,	1,sizeof(opac.ipElement),
			ioOPAC ) - sizeof(opac.ipElement));

		n -= (long)( fwrite( opac.in1,			1,sizeof(opac.in1),
			ioOPAC ) - sizeof(opac.in1));

		n -= (long)( fwrite( opac.ipo3exc,	1,sizeof(opac.ipo3exc),
			ioOPAC ) - sizeof(opac.ipo3exc));

		n -= (long)( fwrite( opac.ipo3exc3,	1,sizeof(opac.ipo3exc3),
			ioOPAC ) - sizeof(opac.ipo3exc3));

		n -= (long)( fwrite( opac.ipo1exc,	1,sizeof(opac.ipo1exc),
			ioOPAC ) - sizeof(opac.ipo1exc));

		n -= (long)( fwrite( &opac.iopo2d,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.ipmgex,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( &opac.ipOpMgEx,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( opac.ica2ex,		1,sizeof(opac.ica2ex),
			ioOPAC ) - sizeof(opac.ica2ex));

		n -= (long)( fwrite( &opac.ica2op,	1,sizeof(long) ,
			ioOPAC ) - sizeof(long));

		n -= (long)( fwrite( rfield.AnuOrg,		1,sizeof(rfield.AnuOrg),
			ioOPAC ) - sizeof(rfield.AnuOrg));

		/* now check whether n is still zero, problem if not */
		if( n!= 0 )
		{
			fprintf( ioQQQ, " problem trying to write opacity.opc pointers\n" );
			fprintf( ioQQQ, " fwrite was short by %li\n",
				n);
			puts( "[Stop in OpacityCreateAll]" );
			cdEXIT(EXIT_FAILURE);
		}

		fclose(ioOPAC);
		fprintf( ioQQQ, "\n\nCompile opacity completed ok  - stopping.\n" );
		fprintf( ioQQQ, "The file opacity.opc was created.\n" );
		fprintf( ioQQQ, "Make sure this file lies somewhere on the path.\n\n\n" );
		puts( "[Stop in OpacityCreateAll]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( lgRealloc )
		fprintf(ioQQQ,
		" Please consider revising ndimOpacityStack in OpacityCreateAll, a total of %li cells were needed.\n\n" , opac.nOpacTot);

#	ifdef DEBUG_FUN
	fputs( " <->OpacityCreateAll()\n", debug_fp );
#	endif
	return;
}
/*OpacityCreatePowerLaw generate array of cross sections using a simple power law fit */
static void OpacityCreatePowerLaw(
	/* lower energy limit on continuum mesh */
	long int ilo, 
	/* upper energy limit on continuum mesh */
	long int ihi, 
	/* threshold cross section */
	double cross, 
	/* power law index */
	double s, 
	/* pointer to opacity offset where this starts */
	long int *ip)
{
	long int i;
	double thres;

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

	/* non-positive cross section is unphysical */
	ASSERT( cross > 0. );

	/* place in the opacity stack where we will stuff cross sections */
	*ip = opac.nOpacTot + 1;
	ASSERT( *ip > 0 );
	ASSERT( ilo > 0 );
	thres = rfield.anu[ilo-1];

	if( opac.nOpacTot + ihi - ilo + 1 > ndimOpacityStack )
		opacity_more_memory();

	for( i=ilo-1; i < ihi; i++ )
	{
		opac.OpacStack[i-ilo+*ip] = cross*pow(rfield.anu[i]/thres,-s);
	}

	opac.nOpacTot += ihi - ilo + 1;

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

/*OpacityCreateReilMan generate photoionization cross sections from Reilman and Manson points */
static void OpacityCreateReilMan(long int low, 
  long int ihi, 
  float cross[], 
  long int ncross, 
  long int *ipop, 
  const char *chLabl)
{
	long int i, 
	  ics, 
	  j, 
	  ncr;

#define	NOP	100
	float cs[NOP], 
	  en[NOP], 
	  slope;

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

	/* this is the opacity entering routine designed for
	 * the Reilman and Manson tables.  It works with incident
	 * photon energy (entered in eV) and cross sections in megabarns
	 * */
	*ipop = opac.nOpacTot + 1;
	ASSERT( *ipop > 0 );

	ncr = ncross/2;
	if( ncr > NOP )
	{
		fprintf( ioQQQ, " Too many opacities were entered into OpacityCreateReilMan.  Increase the value of NOP.\n" );
		fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
		puts( "[Stop in OpacityCreateReilMan]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* the array CROSS has ordered pairs of elements.
	 * the first is the energy in eV (not Ryd)
	 * and the second is the cross section in megabarns */
	for( i=0; i < ncr; i++ )
	{
		en[i] = (float)(cross[i*2]/13.6);
		cs[i] = (float)(cross[(i+1)*2-1]*1e-18);
	}

	ASSERT( low>0 );
	if( en[0] > rfield.anu[low-1] )
	{
		fprintf( ioQQQ, 
			" OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (low fail).\n" );
		fprintf( ioQQQ, 
			" The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n", 
		  rfield.anu[low-1]*EVRYD, en[0]*EVRYD );

		fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
		fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" );
		fprintf( ioQQQ, " " );

		for( i=0; i < ncross; i++ )
		{
			fprintf( ioQQQ, "%11.4e", cross[i] );
		}

		fprintf( ioQQQ, "\n" );
		puts( "[Stop in OpacityCreateReilMan]" );
		cdEXIT(EXIT_FAILURE);
	}

	slope = (cs[1] - cs[0])/(en[1] - en[0]);
	ics = 1;

	if( opac.nOpacTot + ihi - low + 1 > ndimOpacityStack )
	  opacity_more_memory();

	/* now fill in the opacities using linear interpolation */
	for( i=low-1; i < ihi; i++ )
	{
		if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
		{
			opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 
			  en[ics-1]);
		}

		else
		{
			ics += 1;
			if( ics + 1 > ncr )
			{
				fprintf( ioQQQ, " OpacityCreateReilMan: The entered opacity energy bandwidth is not large enough (high fail).\n" );
				fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n", 
				  rfield.anu[i]*13.6, en[ncr-1]*13.6 );
				fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl
				   );
				fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n", 
				  en[0]*13.65 );
				fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n", 
				  rfield.anu[ihi-1]*13.6 );
				fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n", 
				  rfield.anu[low-1]*13.6 );
				puts( "[Stop in OpacityCreateReilMan]" );
				cdEXIT(EXIT_FAILURE);
			}

			slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
			if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
			{
				opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 
				  en[ics-1]);
			}
			else
			{
				ASSERT( i > 0);
				fprintf( ioQQQ, " Internal logical error in OpacityCreateReilMan.\n" );
				fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n", 
				  rfield.anu[i]*13.6, i, en[ics-1], en[ics] );

				fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n", 
				  rfield.anu[i-1]*13.6 );

				fprintf( ioQQQ, " Here comes the energy array.  ICS=%4ld\n", 
				  ics );

				for( j=0; j < ncr; j++ )
				{
					fprintf( ioQQQ, "%10.2e", en[j] );
				}
				fprintf( ioQQQ, "\n" );

				fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
				puts( "[Stop in OpacityCreateReilMan]" );
				cdEXIT(EXIT_FAILURE);
			}
		}
	}
	/* >>chng 02 may 09, this was a significant logcal error */
	/* >>chng 02 may 08, by Ryan.  This routine did not update the total slots filled.	*/
	opac.nOpacTot += ihi - low + 1;

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


/*ofit compute cross sections for all shells of atomic oxygen */
double ofit(double e, 
	  float opart[])
{
	long int i, 
	  _r;

	double otot,
		q, 
	  x;

	static double y[7][5];
	static double eth[7]={13.62,16.94,18.79,28.48,50.,110.5,538.};
	static long l[7]={1,1,1,0,1,1,0};
	static int _aini = 1;

	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		{ static double _itmp0[] = {8.915,3995.,3.242,10.44,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[0][i-1] = _itmp0[_r++];
			}
		}
		{ static double _itmp1[] = {11.31,1498.,5.27,7.319,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[1][i-1] = _itmp1[_r++];
			}
		}
		{ static double _itmp2[] = {10.5,1.059e05,1.263,13.04,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[2][i-1] = _itmp2[_r++];
			}
		}
		{ static double _itmp3[] = {19.49,48.47,8.806,5.983,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[3][i-1] = _itmp3[_r++];
			}
		}
		{ static double _itmp4[] = {50.,4.244e04,0.1913,7.012,4.454e-02};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[4][i-1] = _itmp4[_r++];
			}
		}
		{ static double _itmp5[] = {110.5,0.1588,148.3,-3.38,3.589e-02};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[5][i-1] = _itmp5[_r++];
			}
		}
		{ static double _itmp6[] = {177.4,32.37,381.2,1.083,0.0};
		for( i=1, _r = 0; i <= 5; i++ )
		{
			y[6][i-1] = _itmp6[_r++];
			}
		}
		_aini = 0;
	}

#	ifdef DEBUG_FUN
	fputs( "<+>ofit()\n", debug_fp );
#	endif
	/*compute cross sections for all shells of atomic oxygen
	 * Photoionization of OI
	 * Input parameter:   e - photon energy, eV
	 * Output parameters: otot - total photoionization cross section, Mb
	 *  opart(1) - 2p-shell photoionization, goes to 4So
	 *  opart(2) - 2p-shell photoionization, goes to 2Do
	 *  opart(3) - 2p-shell photoionization, goes to 2Po
	 *  opart(4) - 2s-shell photoionization
	 *  opart(5) - double photoionization, goes to O++
	 *  opart(6) - triple photoionization, goes to O+++
	 *  opart(7) - 1s-shell photoionization */

	otot = 0.0;

	for( i=0; i < 7; i++ )
	{
		opart[i] = 0.0;
	}

	for( i=0; i < 7; i++ )
	{
		if( e >= eth[i] )
		{
			q = 5.5 - 0.5*y[i][3] + l[i];

			x = e/y[i][0];

			opart[i] = (float)(y[i][1]*(POW2(x - 1.0) + POW2(y[i][4]))/
			  pow(x,q)/pow(1.0 + sqrt(x/y[i][2]),y[i][3]));

			otot += opart[i];
		}
	}

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

/******************************************************************************/

/*OpacityCreate1Element generate ionic subshell opacities by calling atmdat_phfit */
static void OpacityCreate1Element(
		  /* atomic number on the C scale, lowest ever called will be Li=2 */
		  long int nelem)
{
	long int ihi, 
	  ip, 
	  ipop, 
	  limit, 
	  low, 
	  need, 
	  nelec, 
	  ion, 
	  nshell;
	double cs; 
	double energy;

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

	/* confirm range of validity of atomic number, Li=2 should be the lightest */
	ASSERT( nelem >= 2 );
	ASSERT( nelem < LIMELM );

	/*>>chng 99 jan 27, no longer redo hydrogenic opacity here */
	/* for( ion=0; ion <= nelem; ion++ )*/
	for( ion=0; ion < nelem; ion++ )
	{

		/* will be used for a sanity check on number of hits in a cell*/
		for( ip=0; ip < rfield.nupper; ip++ )
		{
			opac.opacity_abs[ip] = 0.;
		}

		/* number of bound electrons */
		nelec = nelem+1 - ion;

		/* loop over all shells, from innermost K shell to valence */
		for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
		{
			/* this is array index for start of this shell within large opacity stack */
			opac.ipElement[nelem][ion][nshell][2] = opac.nOpacTot +  1;

			/* this is continuum upper limit to array index for energy range of this shell */
			limit = opac.ipElement[nelem][ion][nshell][1];

			/* this is number of cells in continuum needed to store opacity */
			need = limit - opac.ipElement[nelem][ion][nshell][0] + 1;

			/* check that opac will have enough frequeny cells */
			if( opac.nOpacTot + need > ndimOpacityStack )
				opacity_more_memory();

			/* set lower and upper limits to this range */
			low = opac.ipElement[nelem][ion][nshell][0];
			ihi = opac.ipElement[nelem][ion][nshell][1];
			ipop = opac.ipElement[nelem][ion][nshell][2];

			/* make sure indices are within correct bounds,
			 * mainly check on logic for detecting missing shells */
			ASSERT( low <= ihi || low<5 );

			/* loop over energy range of this shell */
			for( ip=low-1; ip < ihi; ip++ )
			{
				/* photo energy MAX so that we never eval below threshold */
				energy = MAX2(rfield.AnuOrg[ip]*EVRYD , 
					PH1COM.PH1[nshell][nelec-1][nelem][0]);

				/* the cross section in mega barns */
				cs = atmdat_phfit(nelem+1,nelec,nshell+1,energy,PhFitOn.lgPhFit);
				/* cannot assert that cs is positive since, at edge of shell,
				 * energy might be slightly below threshold and hence zero,
				 * due to finite size of continuum bins */
				opac.OpacStack[ip-low+ipop] = cs*1e-18;

				/* add this to total opacity, which we will confirm to be greater than zero below */
				opac.opacity_abs[ip] += cs;
			}

			opac.nOpacTot += ihi - low + 1;

			/* punch pointers option */
			if( punch.lgPunPoint )
			{
				fprintf( punch.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n", 
				  nelem, ion, nshell, rfield.anu[low-1], rfield.anu[ihi-1], 
				  opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] );
			}
		}

		/*confirm that total opacity is greater than zero  */
		for( 
			ip=opac.ipElement[nelem][ion][Heavy.nsShells[nelem][ion]-1][0]-1;
		    ip < KshllEnr.KshellLimit; ip++ )
		{
			ASSERT( opac.opacity_abs[ip] > 0. );
		}

	}

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

/*opacity_more_memory allocate more memory for opacity stack */
static void opacity_more_memory(void)
{

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

	/* double size */
	ndimOpacityStack *= 2;
	if( (opac.OpacStack = (double *)realloc(  opac.OpacStack , (size_t)ndimOpacityStack*sizeof(double) ))==NULL )
		BadMalloc();
	fprintf( ioQQQ, " NOTE OpacityCreate1Element needed more opacity cells than ndimOpacityStack,  please consider increasing it.\n" );
	fprintf( ioQQQ, " NOTE OpacityCreate1Element doubled memory allocation to %li.\n",ndimOpacityStack );
	lgRealloc = TRUE;

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

/*Opacity_iso_photo_cs returns photoionization cross section for isoelectronic sequences */
static double Opacity_iso_photo_cs( 
		/* photon energy ryd */
		float energy , 
		/* iso sequence */
		long ipISO , 
		/* charge, 0 for H */
		long nelem , 
		/* n - meaning depends on iso */
		long n )
{
	/* >>chng 01 dec 23, from float to double */
	double thres/*,ejected_electron_energy*/;
	double crs=-DBL_MAX;

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

	if( ipISO==ipH_LIKE )
	{
		double photon;
		if( n==0 )
		{
			/* this is the ground state, use Dima's routine, which works in eV
			 * and returns megabarns */
			thres = MAX2(energy*(float)EVRYD , PH1COM.PH1[0][0][nelem][0]);
			crs = atmdat_phfit(nelem+1,1,1,thres,PhFitOn.lgPhFit)* 1e-18;
			/* make sure cross section is reasonable */
			ASSERT( crs > 0. && crs < 1e-10 );
		}
		else if( n==1 )
		{
			/* 2s */
			/* Changed this to ejected electron energy because H_photo_cs, when it was 
			 * being passed a photon energy, subtracted off a threshold energy that was
			 * hydrogenic.  But the thresholds of most helike levels, having positive 
			 * quantum defects, are actually less than hydrogenic.  This created a clumsy
			 * situation in which the photon energy passed to H_photo_cs had to be bumped
			 * up by about 1.001, which still did not guarantee that H_photo_cs would find
			 * a meaningful, positive electron energy for an arbitrary level and ion.  So
			 * H_photo_cs was changed to accept ejected electron energy instead, with its
			 * functionality completely unchanged.	*/
			
			/*ejected_electron_energy = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]) 
				- iso.xIsoLevNIonRyd[ipISO][nelem][n];	*/

			/* >>chng 02 apr 24 - protect against too small a photon energy */
			/* This line doesn't do anything, change photon to energy		
			photon = MAX2( POW2( (double)(nelem+1) / 2. ) , energy );*/
			energy = MAX2( POW2( (float)(nelem+1) / 2.f ) , energy );


			/* photon will now be relative to threshold */
			photon = energy / iso.xIsoLevNIonRyd[ipISO][nelem][n];
			photon = MAX2( photon , 1. + FLT_EPSILON*2. );
		
			/*crs = H_photo_cs( ejected_electron_energy + SMALLFLOAT, 2,0,nelem+1 );*/
			/* >>chng 02 apr 24, from electron to photon energy */
			crs = H_photo_cs( photon , 2,0,nelem+1 );
			/* make sure cross section is reasonable */
			ASSERT( crs > 0. && crs < 1e-10 );
		}
		else if( n==2 )
		{
			/* 2p */
			/* same as above	*/
			
			/*ejected_electron_energy = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]) 
				- iso.xIsoLevNIonRyd[ipISO][nelem][n];	*/

			/* >>chng 02 apr 24 - protect against too small a photon energy */
		 	/* This line doesn't do anything, change photon to energy		
			photon = MAX2( POW2( (double)(nelem+1) / 2. ) , energy );	*/
			energy = MAX2( POW2( (float)(nelem+1) / 2.f ) , energy );

			/* photon will now be relative to threshold */
			photon = energy / iso.xIsoLevNIonRyd[ipISO][nelem][n];
			photon = MAX2( photon , 1. + FLT_EPSILON*2. );
		
			/* >>chng 02 apr 24, from electron to photon energy */
			/*crs = H_photo_cs( ejected_electron_energy + SMALLFLOAT, 2,1,nelem+1 );*/
			crs = H_photo_cs( photon , 2,1,nelem+1 );
			/* make sure cross section is reasonable */
			ASSERT( crs > 0. && crs < 1e-10 );
		}
		else
		{
			/* for first cell, depending on the curent resolution of the energy mesh,
			 * the center of the first cell can be below the ionization limit of the
			 * level.  do not let the energy fall below this limit */
			/* This will make sure that we don't call epsilon below threshold,
		     * the factor 1.001 was chosen so that atmdat_hpfit, which works in terms od Dima's rydberg constant,
		     * is not tripped below threshold */
			thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.001f );

			crs = atmdat_hpfit(nelem+1,n,thres*EVRYD);
			/* make sure cross section is reasonable */
			ASSERT( crs > 0. && crs < 1e-10 );
		}
	}
	else if( ipISO==ipHE_LIKE )
	{
		thres = MAX2( energy , iso.xIsoLevNIonRyd[ipISO][nelem][n]);
		/* this would be a collapsed level */
		if( n >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] )
		{
			long int nup = iso.n_HighestResolved_max[ipHE_LIKE][nelem] + n + 1 -
				(iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem]);

			/* this is a collapsed level - this is hydrogenic routine and
			 * first he-like energy may not agree exactly with threshold for H */
			crs = atmdat_hpfit(nelem,nup ,thres*EVRYD);
			/* make sure cross section is reasonable if away from threshold */
			ASSERT( 
				(energy < iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.02) ||
				(crs > 0. && crs < 1e-10) );
		}
		else
		{

			/*thres = MAX2( iso.xIsoLevNIonRyd[ipISO][nelem][n]*1.01f , energy );
			He_cross_section returns cross section (cm^-2), 
			* given EgammaRyd, the photon energy in Ryd,
			* ipLevel, the index of the level, 0 is ground, 3 within 2 3P,
			* nelem is charge, equal to 1 for Helium,
			* this is a wrapper for cross_section */
			crs = He_cross_section( thres , n , nelem );
					/* make sure cross section is reasonable */
			ASSERT( crs > 0. && crs < 1e-10 );
		}
	}
	else
		TotalInsanity();

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

/*hmiopc derive total H- H minus opacity */
#ifdef NCRS
#	undef NCRS
#endif
#define	NCRS	33

static double hmiopc(double freq)
{
	double energy, 
	  hmiopc_v, 
	  x, 
	  y;
	static double y2[NCRS];
	static double crs[NCRS]={0.,0.124,0.398,0.708,1.054,1.437,1.805,
	  2.176,2.518,2.842,3.126,3.377,3.580,3.741,3.851,3.913,3.925,
	  3.887,3.805,3.676,3.511,3.306,3.071,2.810,2.523,2.219,1.898,
	  1.567,1.233,.912,.629,.39,.19};
	static double ener[NCRS]={0.,0.001459,0.003296,0.005256,0.007351,
	  0.009595,0.01201,0.01460,0.01741,0.02044,0.02375,0.02735,0.03129,
	  0.03563,0.04043,0.04576,0.05171,0.05841,0.06601,0.07469,0.08470,
	  0.09638,0.1102,0.1268,0.1470,0.1723,0.2049,0.2483,0.3090,0.4001,
	  0.5520,0.8557,1.7669};
	static int lgFirst = TRUE;

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

	/* bound free cross section (x10**-17 cm^2) from Doughty et al
	 * 1966, MNRAS 132, 255; good agreement with Wishart MNRAS 187, 59p. */

	/* photoelectron energy, add 0.05552 to get incoming energy (Ryd) */


	if( lgFirst )
	{
		/* set up coefficients for spline */
		spline(ener,crs,NCRS,2e31,2e31,y2);
		lgFirst = FALSE;
	}

	energy = freq - 0.05552;
	if( energy < ener[0] || energy > ener[NCRS-1] )
	{
		hmiopc_v = 0.;
	}
	else
	{
		x = energy;
		splint(ener,crs,y2,NCRS,x,&y);
		hmiopc_v = y*1e-17;
	}

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

/*rayleh compute Rayleigh scattering cross section for Lya */
static double rayleh(double ener)
{
	double rayleh_v;

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

	/*TODO	2	update to astro-ph/0308073, Lee, H-W, ApJ in press */
	/* do hydrogen Rayleigh scattering cross sections;
	 * fits to 
	 *>>refer	Ly	scattering	Gavrila, M., 1967, Physical Review 163, 147
	 * and Mihalas radiative damping
	 *
	 * >>chng 96 aug 15, changed logic to do more terms for each part of
	 * rayleigh scattering
	 * if( ener.lt.0.05 ) then
	 *  rayleh = 8.41e-25 * ener**4 * DampOnFac
	 * */
	if( ener < 0.05 )
	{
		rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6))*
		  hydro.DampOnFac;
	}

	else if( ener < 0.646 )
	{
		rayleh_v = (8.41e-25*powi(ener,4) + 3.37e-24*powi(ener,6) + 
		  4.71e-22*powi(ener,14))*hydro.DampOnFac;
	}

	else if( ener >= 0.646 && ener < 1.0 )
	{
		rayleh_v = fabs(0.74959-ener);
		rayleh_v = 1.788e5/POW2(FR1RYD*MAX2(0.001,rayleh_v));
		/*  typical energy between Ly-a and Ly-beta */
		rayleh_v = MAX2(rayleh_v,1e-24)*hydro.DampOnFac;
	}

	else
	{
		rayleh_v = 0.;
	}

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


#undef	NCSH2P
#undef NCRS
#undef NMAGIC
#undef MAGIC

