/*opinb add opacity of individual species, including stimulated emission */
#include "cddefines.h"
#include "opac.h"
#include "rfield.h"
#include "fluors.h"
#include "hphoto.h"
#include "opinb.h"

void opinb(
	/* pointer to opacity within opacity stack */
	long int ipOpac, 
	/* pointer to low energy in continuum array for this opacity band */
	long int ipLowEnergy, 
	/* pointer to high energy in continuum array for this opacity */
	long int ipHiEnergy, 
	/* this abundance of this species, may be zero */
	float abundance, 
	/* the departure coef, may be infinite or zero */
	double DepartCoef ,
	/* either 'v' for volitile or 's' for static opacities */
	char chStat )
{
	long int i, 
	  iup, 
	  k;
	/* will be used to store inverse of departure coef */
	float DepartCoefInv;

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

	/* add opacity of individual species, including stimulated emission
	 * abundance is the density of the lower level (cm^-3)
	 * DepartCoef is its departure coefficient, can be zero */

	/* this is opacity offset, must be positive */
	assert( ipOpac > 0 );

	/* check that chStat is either 'v' or 's' */
	assert( chStat == 'v' || chStat == 's' );

	/* do nothing if abundance is zero, or if static opacities do not
	 * need to be redone */
	if( abundance <= 0. || (chStat=='s' && !opac.lgRedoStatic) )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->opinb()\n", debug_fp );
#		endif
		return;
	}

	k = ipOpac - ipLowEnergy;

	/* DepartCoef is dep coef, lgFluor is turned off with 'no indcued' command */
	if( (DepartCoef > 1e-35 && fluors.lgFluor) && HPhoto.lgHInducImp )
	{
		iup = MIN2(ipHiEnergy,rfield.nflux);
		/* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/
		/*iup = MAX2(ipLowEnergy,iup);*/
		DepartCoefInv = (float)(1./DepartCoef);
		if( chStat == 'v' )
		{
			/* volitile opacities, always reevaluate */
			for( i=ipLowEnergy-1; i < iup; i++ )
			{
				opac.opac[i] += opac.OpacStack[i+k]*abundance*
					MAX2(0.f , 1.f-  rfield.ContBoltz[i]*DepartCoefInv);
			}
		}
		else
		{
			/* static opacities, save in special array */
			for( i=ipLowEnergy-1; i < iup; i++ )
			{
				opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance*
					MAX2(0.f , 1.f-  rfield.ContBoltz[i]*DepartCoefInv);
			}
		}
	}

	else
	{
		/* DepartCoef is the departure coef, which can go to zero, 
		 * neglect stimulated emission in this case */
		iup = MIN2(ipHiEnergy,rfield.nflux);
		/* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/
		/*iup = MAX2(ipLowEnergy,iup);*/
		if( chStat == 'v' )
		{
			for( i=ipLowEnergy-1; i < iup; i++ )
			{
				opac.opac[i] += opac.OpacStack[i+k]*abundance;
			}
		}
		else
		{
			for( i=ipLowEnergy-1; i < iup; i++ )
			{
				opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance;
			}
		}
	}


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