/* 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 */
/*OpacityAdd1Element enter total photo cross section for all subshells into opacity array */
#include "cddefines.h"
#include "iso.h"
#include "rfield.h"
#include "dense.h"
#include "heavy.h"
#include "opacity.h"

void OpacityAdd1Element(
		/* nelem is 0 for H, 1 for He, etc */
		long int nelem)
{
	long int ipHi, 
	  ipop, 
	  limit,
	  low, 
	  n, 
	  ion, 
	  nshell;
	char chStat;
	double abundance;

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

	/* this routine drives OpacityAdd1Subshell to put in total opacities for all shells*/

	/*begin sanity check */
	ASSERT( (nelem >=0 ) && (nelem < LIMELM) );

	/* first do simple two-level systems -
	 * this is number of species that are not treated on common iso-electronic series */
	limit = nelem + 1 - NISO;
	/* this can be called with hydrogen itself, in which case nelem is 0, and limit is
	 * -1 - do not do any of the simple ions */
	limit = MAX2( 0 , limit );

	/* do not include the ion stages that have complete atoms,
	 * currently H and He like iso sequences */
	for( ion=0; ion < limit; ion++ )
	{
		if( dense.xIonDense[nelem][ion] > 0. )
		{
			/*start with static opacities, then do volatile*/
			
			chStat = 's';
			/* number of bound electrons */
			for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
			{
				/* highest shell will be volatile*/
				if( nshell== Heavy.nsShells[nelem][ion]-1 )
					chStat = 'v';
				/* set lower and upper limits to this range */
				low = opac.ipElement[nelem][ion][nshell][0];
				ipHi = opac.ipElement[nelem][ion][nshell][1];
				ipop = opac.ipElement[nelem][ion][nshell][2];
				/* OpacityAdd1Subshell will not do anything if static opacities do not need to be reset*/
				OpacityAdd1Subshell(ipop,low,ipHi,dense.xIonDense[nelem][ion] , chStat );
			}
		}
	}

	/* now loop over all species done as large multi-level systems */
	/* >>chng 02 jan 17, add loop over H and He like */
	/* ion is on the c scale, =0 for HI, =1 for HeII */
	for( ion=limit; ion<nelem+1; ++ion )
	{
		/* ipISO is 0 for H-like, 1 for He-like */
		long int ipISO = nelem-ion;

		/* do multi level systems, but only if present 
		 * test for nelem+1 in case atom present but not ion, test is whether the
		 * abundance of the recombined species is present */
		/* >>chng 02 jan 17, sec dim had been nelem+1, change to ion+1 */
		/*if( dense.xIonDense[nelem][nelem] > 0. )*/
		if( dense.xIonDense[nelem][ion] > 0. )
		{
			/* do ground first, then all excited states */
			n = 0;
			/* abundance of recombined species, which can be zero if no ion present */
			/* >>chng 02 jan 17, to arbitrary iso sequence */
			/*abundance = iso.Pop2Ion[ipH_LIKE][nelem][n]*dense.xIonDense[nelem][nelem+1];*/
			abundance = iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][ion+1];

			/* >>chng 02 may 06, add second test, had been just the chck on helium,
			 * with no option to use new soln */
			if( abundance == 0.  )
			/* >>chng 02 jan 17, logical error here, must use dense.xIonDense if abundance is 0 */
			{
				/* no ionized species, assume everything in ground */
				/* >>chng 02 jan 17, to arbitrary iso sequence */
				abundance = dense.xIonDense[nelem][ion];
			}

			/* >>chng 02 jan 17, to arbitrary iso sequence */
			/* use computed opacities and departure coef for level */
			OpacityAdd1SubshellInduc(
				iso.ipOpac[ipISO][nelem][n],
				iso.ipIsoLevNIonCon[ipISO][nelem][n],
				/* the upper limit to the integration, 
				* ground opacity goes up to the high energy limit of code*/
				rfield.nflux,
				/* the abundance of the ion */
				abundance,
				/* departure coef, volatile opac, always reevaluate */
				iso.DepartCoef[ipISO][nelem][n] , 'v' );

			/* do excited levvels,
			 * this loop only if upper levels have finite population*/
			/* >>chng 02 jan 17, to arbitrary iso sequence */
			/*if( iso.Pop2Ion[ipH_LIKE][nelem][3]*dense.xIonDense[nelem][nelem+1] > 0. )*/
			if( iso.Pop2Ion[ipISO][nelem][3]*dense.xIonDense[nelem][ion+1] > 0. )
			{
				char chType = 'v';
				/* always want to evaluate all opacities for n=3, 4, use static opacities for higher levels */
				/* >>chng 02 jan 17, to arbitrary iso sequence, use single loop, setting type of opacity */
				for( n=1; n < iso.numLevels_max[ipISO][nelem]; n++ )
				{
					/* include correction for stimulated emission */
					OpacityAdd1SubshellInduc(
						iso.ipOpac[ipISO][nelem][n],
						iso.ipIsoLevNIonCon[ipISO][nelem][n],
						/* the high energy bound of excited states is the 
						* edge of the Lyman continuum */
						iso.ipIsoLevNIonCon[ipISO][nelem][0],
						iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][ion+1],
						/* departure coef, volitile opacities */
						iso.DepartCoef[ipISO][nelem][n] , chType );
					/* above 4 is static */
					if( n==4)
						chType = 's';
				}
			}
		}
	}

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

