/*PutOpacity enter total photo cross section for all subshells into opacity array */
#include "cddefines.h"
#include "opacpoint.h"
#include "hydrogenic.h"
#include "nh.h"
#include "rfield.h"
#include "nsshells.h"
#include "ionfracs.h"
#include "opacin.h"
#include "opinb.h"
#include "opinpl.h"
#include "putopacity.h"

void PutOpacity(
		/* ipZ is 0 for H, 1 for He, etc */
		long int ipZ)
{
	long int ipHi, 
	  ipop, 
	  low, 
	  n, 
	  nion, 
	  nshell;
	char chStat;
	double abundance;

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

	/* this routine drives opacin to put in total opacities */

	/*begin sanity check */
	assert (ipZ >=0 );
	assert( ipZ < LIMELM );

	/* do not include the hydrogenic level */
	for( nion=0; nion < ipZ; nion++ )
	{
		if( xIonFracs[nion+1][ipZ] > 0. )
		{
			/*start with static opacities, then go to volatile*/
			chStat = 's';
			/* number of bound electrons */
			for( nshell=0; nshell < nsShellsCom.nsShells[nion][ipZ]; nshell++ )
			{
				/* highest shell will be volatile*/
				if( nshell== nsShellsCom.nsShells[nion][ipZ]-1 )
					chStat = 'v';
				/* set lower and upper limits to this range */
				low = OpacPoint.ipElement[0][nshell][nion][ipZ];
				ipHi = OpacPoint.ipElement[1][nshell][nion][ipZ];
				ipop = OpacPoint.ipElement[2][nshell][nion][ipZ];
				/* opacin will not do anything if static opacities do not need to be reset*/
				opacin(ipop,low,ipHi,xIonFracs[nion+1][ipZ] , chStat );
			}
		}
	}

	/* now do hydrogenic, but only if present 
	 * test for ipZ+1 in case atom present but not ion, test is whether the
	 * abundance of the recombined species is present */
	if( xIonFracs[ipZ+1][ipZ] > 0. )
	{
		/* do ground first, then all excited states */
		n = IP1S;
		/* abundance of recombined species, which can be zero if no ion present */
		abundance = hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ];
		if( abundance == 0. )
		{
			abundance = hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ];
		}
		else 
		{
			/* no ionized species, assume everything in ground */
			abundance = xIonFracs[ipZ+1][ipZ];
		}

		/* use computed opacities and departure coef for level */
		opinb(
			OpacPoint.ipHOpac[ipZ][n],
			nh.ipHn[ipZ][n],
			/* the upper limit to the integration, 
			 * ground opacity goes up to the high energy limit of code*/
			rfield.nflux,
			hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ],
			/* departure coef, volatile opac, always reevaluate */
			hydro.hbn[ipZ][n] , 'v' );

		/* do excited levvels,
		 * this loop only if upper levels have finite population*/
		if( hydro.hn[ipZ][3]*xIonFracs[ipZ+2][ipZ] > 0. )
		{
			/* always want to evaluate all opacities for n=3, 4 */
			for( n=IP2S; n < 5; n++ )
			{
				/* this is for lower levels of H, uses computed opacities */
				/* this is departure coef for level */
				opinb(
					OpacPoint.ipHOpac[ipZ][n],
					nh.ipHn[ipZ][n],
					/* the high energy bound of excited states is the 
					 * edge of the Lyman continuum */
					nh.ipHn[ipZ][IP1S],
					hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ],
					/* departure coef, volitile opacities */
					hydro.hbn[ipZ][n] , 'v' );
			}

			/* NHPLPHOT lowest level where hydrogenic nu^-3 photo cross sect used*/
			for( n=5; n < NHPLPHOT; n++ )
			{
				/* this is for lower levels of H, uses computed opacities */
				/* this is departure coef for level */
				opinb(
					OpacPoint.ipHOpac[ipZ][n],
					nh.ipHn[ipZ][n],
					/* the high energy bound of excited states is the 
					 * edge of the Lyman continuum */
					nh.ipHn[ipZ][IP1S],
					hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ],
					/* departure coef, static opacities */
					hydro.hbn[ipZ][n] , 's' );
			}

			/* now do upper levels that are only powerlaws */
			for( n=NHPLPHOT; n <= nhlevel; n++ )
			{
				/* this is for highly excited states, where simple
				 * power laws are used */
				opinPL(
					n,
					ipZ,
					hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ],
					hydro.hbn[ipZ][n]);
			}
		}
	}


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