/*opinPL add opacity of individual species, with stimulated emission, for power law opacity */
#include "cddefines.h"
#include "hydrophocs.h"
#include "opac.h"
#include "nh.h"
#include "rfield.h"
#include "fluors.h"
#include "hphoto.h"
#include "opinpl.h"

void opinPL(
	long int n,	/* quantum number of level on C not physical scale, 
					 * must be less than NHYDRO_MAX_LEVEL, dim of STH */
  long int ipZ,/* charge, ipZ=0 for h, 1 for he, etc */
  float abundance, 
  double b)
{
	long int i, 
	  ipHi;
	float binv, 
	  csThresh;

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

	/*add opacity of individual species, with stimulated emission, for power law opacity
	 *abundance is the density of the lower level (cm^-3)
	 *b is its departure coefficient, can be zero, is double precision */

	/* do nothing if abundance is zero, 
	 * or if no need to redo static opacities */
	if( abundance <= 0. || !opac.lgRedoStatic )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->opinPL()\n", debug_fp );
#		endif
		return;
	}

	/* following is numer of levels known - in HydroPhoCS.STH, where threshold
	 * photoionization cross sections are stored */
	assert( n < NHYDRO_MAX_LEVEL );

	csThresh = HydroPhoCS.STH[n]*rfield.anu3[nh.ipHn[ipZ][n]-1]/
	  POW2(ipZ+1.f);

	csThresh *= abundance;

	/* pick highest continuum to fill in */
	if( n == IP1S )
	{
		/* for Lyman continuum, go up to highest defined boltzmann factor */
		ipHi = rfield.nflux;
	}
	else
	{
		/* for excited states go either up to Lyman continuum or highest
		 * defined Boltzmann factor */
		ipHi = MIN2(rfield.nflux,nh.ipHn[ipZ][IP1S]);
	}

	/* b is dep coef, lgFluor is turned off with 'no indcued' command */
	if( (b > 1e-35 && fluors.lgFluor) && HPhoto.lgHInducImp )
	{
		binv = (float)(1./b);
		for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ )
		{
			opac.OpacStatic[i] += csThresh/rfield.anu3[i]*
				MAX2(0.f,1.f- rfield.ContBoltz[i]*binv);
		}
	}
	else
	{
		/* B is the departure coef, which can go to zero, neglect stim emin */
		for( i=nh.ipHn[ipZ][n]-1; i < ipHi; i++ )
		{
			opac.OpacStatic[i] += csThresh/rfield.anu3[i];
		}
	}

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