/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*OpacityAddTotal derive total opacity for this position,
 * called by ConvBase */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "iso.h"
#include "grainvar.h"
#include "ca.h"
#include "rfield.h"
#include "oxy.h"
#include "co.h"
#include "dense.h"
#include "atoms.h"
#include "converge.h"
#include "ionbal.h"
#include "taulines.h"
#include "trace.h"
#include "hmi.h"
#include "phycon.h"
#include "opacity.h"
#define RJRWLOOP 1
#define RJRWMACRO 1

void OpacityAddTotal(void)
{
	long int i, 
	  ion,
	  ipHi,
	  ipISO,
	  ipop,
	  limit,
	  low,
	  nelem,
	  n,
	  nshell; 
	double DepartCoefInv ,
	  fac, 
	  sum;
	float SaveOxygen1 , 
		SaveCarbon1;

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

	/* OpacityZero will zero out scattering and absorption opacities,
	 * and set OldOpacSave to opac to save it */
	OpacityZero();

	/* free electron scattering opacity, Compton recoil energy loss */
	/* >>chng 02 mar 26, make this loop only one over electrons alone */
	/*for( i=0; i < (ionbal.ipCompRecoil[ipHYDROGEN][0] - 1); i++ )*/
	for( i=0; i < rfield.nflux; i++ )
	{
		/* scattering part of total opacity */
		opac.opacity_sct[i] += opac.OpacStack[i-1+opac.iopcom]*
		  dense.eden;
	}

	/* opacity due to compton bound recoil ionization */
	/* >>chng 01 dec 19, rewrite as loop over all elements */
	/* >>chng 02 mar 26, add in ions */
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		if( dense.lgElmtOn[nelem] )
		{ 
			for( ion=0; ion<nelem+1; ++ion )
			{
				float factor = dense.xIonDense[nelem][ion];
				if( factor > 0. )
				{
					factor *= ionbal.nCompRecoilElec[nelem-ion];
					for( i=ionbal.ipCompRecoil[nelem][ion]-1; i < rfield.nflux; i++ )
					{
						/* add in bound hydrogen electron scattering, treated as absorption */
						opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor;
					}
				}
			}
		}
	}

	/* opacity due to pair production */
	sum = dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM];
	OpacityAdd1Subshell(opac.ioppr,opac.ippr,rfield.nflux,(float)sum,'s');

	/* free-free free free brems emission for all ions */

	/* >>chng 02 jul 21, use full set of ions and gaunt factor */
	/* ipEnergyBremsThin is index to energy where gas goes optically thin to brems,
		* so this loop is over optically thick frequencies */
	if (RJRWLOOP) 
	{
		static double *TotBremsAllIons;
		double BremsThisIon, bfac = (dense.eden/1e20)/phycon.sqrte/1e10, bhfac;
		static int lgFirstTime=TRUE;
		if( lgFirstTime )
		{
			/* rfield.nupper will not change in one coreload, so just malloc this once */
			TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
			lgFirstTime = FALSE;
		}
		
		bhfac = bfac*(dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]);
		for( i=0; i < rfield.nflux; i++ )
		{
			/* do hydrogen first, before main loop since want to add on H- brems */
			nelem = ipHYDROGEN;
			ion = 1;/* >>chng 02 nov 07, add charged ions as H+ */

			/* for case of hydrogen, add H- brems - OpacStack contains the ratio
			 * of the H- to H brems cross section - multiply by this and the fraction in ground */
			BremsThisIon = bhfac * rfield.gff[ion][i]*
				(1.+opac.OpacStack[i-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
			TotBremsAllIons[i] = BremsThisIon;
			/* there is at least one standard test (grainlte) which has ZERO ionization -
			 * this assert will pass that test (the ==0) and also the usual case */
			ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons[i] > 0.) );
		}
		
		/* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
		 * using gaunt factors from rfield.gff.	*/
		for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
		{
			/* for H, IonLow is 0 and IonHigh is 1 */
			/* MAX2 occurs because we want to start at first ion (or above)
			 * and do not want atom */
			for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
			{
				/* eff. charge is ion, so first rfield.gff argument must be "ion".	*/
				BremsThisIon = (double)ion*ion*dense.xIonDense[nelem][ion]*bfac;
				for( i=0; i < rfield.nflux; i++ )
				{
					TotBremsAllIons[i] += BremsThisIon*rfield.gff[ion][i];
				}
			}
		}
		
		for( i=0; i < rfield.nflux; i++ )
		{
			if( rfield.ContBoltz[i] < 0.995 )
			{
				TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
					(1. - rfield.ContBoltz[i]);
			}
			else
			{
				TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
					rfield.anu[i]*TE1RYD/phycon.te;
			}
			opac.FreeFreeOpacity[i] = TotBremsAllIons[i];
			/* >>chng 02 jul 23, from >0 to >=0, some models have no ionization,
			 * like grainlte.in */
			ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );
			opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
		}
	}
	else
	{
		for( i=0; i < rfield.nflux; i++ )
		{
			double TotBremsAllIons = 0., BremsThisIon ;
			
			/* do hydrogen first, before main loop since want to add on H- brems */
			nelem = ipHYDROGEN;
			ion = 1;/* >>chng 02 nov 07, add charged ions as H+ */
			BremsThisIon = (dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]) *
				POW2( (double)ion )/1e10*rfield.gff[ion][i] *(dense.eden/1e20)/phycon.sqrte;
			
			/* for case of hydrogen, add H- brems - OpacStack contains the ratio
			 * of the H- to H brems cross section - multiply by this and the fraction in ground */
			BremsThisIon *= (1.+opac.OpacStack[i-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
			TotBremsAllIons = BremsThisIon;
			/* there is at least one standard test (grainlte) which has ZERO ionization -
			 * this assert will pass that test (the ==0) and also the usual case */
			ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons > 0.) );
			
			/* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
			 * using gaunt factors from rfield.gff.	*/
			for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
			{
				/* for H, IonLow is 0 and IonHigh is 1 */
				/* MAX2 occurs because we want to start at first ion (or above)
				 * and do not want atom */
				for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
				{
					/* eff. charge is ion, so first rfield.gff argument must be "ion".	*/
					BremsThisIon = POW2( (double)ion )*dense.xIonDense[nelem][ion]/1e10*rfield.gff[ion][i] *
						(dense.eden/1e20)/phycon.sqrte;
					TotBremsAllIons += BremsThisIon;
				}
			}
			
			if( rfield.ContBoltz[i] < 0.995 )
			{
				TotBremsAllIons *= opac.OpacStack[i-1+opac.ipBrems]*
					(1. - rfield.ContBoltz[i]);
			}
			else
			{
				TotBremsAllIons *= opac.OpacStack[i-1+opac.ipBrems]*
					rfield.anu[i]*TE1RYD/phycon.te;
			}
			opac.FreeFreeOpacity[i] = TotBremsAllIons;
			/* >>chng 02 jul 23, from >0 to >=0, some models have no ionization,
			 * like grainlte.in */
			ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );
			opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
		}
	}

	/* H minus absorption, with correction for stimulated emission */
	if( hmi.hmidep > SMALLFLOAT )
	{
		DepartCoefInv = 1./hmi.hmidep;
	}
	else
	{
		/* the hmidep departure coef can become vastly small in totally
		 * neutral gas (no electrons) */
		DepartCoefInv = 1.;
	}

	limit = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
	if (limit > rfield.nflux)
		limit = rfield.nflux;

	for( i=hmi.iphmin-1; i < limit; i++ )
	{
#if RJRWMACRO
		double factor;
		factor = 1. - rfield.ContBoltz[i]*DepartCoefInv;
		if (factor > 0)
#if 1
			opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
											 hmi.Hmolec[ipMHm]*factor);
#else
			opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
											 hmi.hminus*factor);
#endif
#else
		opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
		  hmi.hminus*MAX2(0.,(1. - rfield.ContBoltz[i]*DepartCoefInv)));
#endif
	}

	/* H2P h2plus bound free opacity */
#if RJRWMACRO
	limit = opac.ih2pnt[1]; 
	if (limit > rfield.nflux)
		limit = rfield.nflux;
#else
	limit = MIN2(rfield.nflux,opac.ih2pnt[1]);
#endif
	for( i=opac.ih2pnt[0]-1; i < limit; i++ )
	{
		opac.opacity_abs[i] += hmi.Hmolec[ipMH2p]*opac.OpacStack[i-opac.ih2pnt[0]+
		  opac.ih2pof];
	}

	/* get total population of hydrogen ground to do Rayleigh scattering */
	if( dense.xIonDense[ipHYDROGEN][1] <= 0. )
	{
		fac = dense.xIonDense[ipHYDROGEN][0];
	}
	else
	{
		fac = dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
	}

	/* Ly a damp wing opac (Rayleigh scattering) */
#if RJRWMACRO
	limit = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];
	if (limit > rfield.nflux)
		limit = rfield.nflux;
#else
	limit = MIN2(rfield.nflux, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]);
#endif
	for( i=0; i < limit; i++ )
	{
		opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]);
	}

	 /* remember largest correction for stimulated emission */
	if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 && !conv.lgSearch )
	{
#if RJRWMACRO
		float factor;
		factor = (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]);
		if (opac.stimax[0] < factor)
			opac.stimax[0] = factor;
#else
		opac.stimax[0] = 
			MAX2(opac.stimax[0],
			     (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]));
#endif
	}

	if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p] > 1e-30 && !conv.lgSearch )
	{
#if RJRWMACRO
		float factor;
		factor = (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]);
		if (opac.stimax[1] < factor)
			opac.stimax[1] = factor;
#else
		opac.stimax[1] = 
			MAX2(opac.stimax[1],
			     (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]));
#endif
	}

	/* check whether hydrogen or Helium singlets mased, if not in search mode */
	if( !conv.lgSearch )
	{
		if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc < 0. )
		{
			hydro.lgHLyaMased = TRUE;
		}
	}

	/* following loop adds standard opacities for first 30 elements 
	 * most heavy element opacity added here */

	/* add in high-energy opacity of CO using C and O atomic opacities by
	 * creating fake atomic abundance that includes this */
	SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
	SaveCarbon1 = dense.xIonDense[ipCARBON][0];

	/*TODO	2	update photoelectric opacity for H2 to include real cross sections and energies.
	 * this is not a higher priority because when H2 is formed there can be very little ionizing
	 * radiation. this process must be trivial compared with the solomon process */
	/* >>chng 02 jan 16, approximate inclusion of H_2 photoelectric opacity */
	/* >>refer	H2	photo cs	Yan, M., Sadeghpour, H.R., & Dalgarno, A., 1998, ApJ, 496, 1044 */
	/* include H_2 in total photoelectric opacity */
	/* set lower and upper limits to this range */
	ion = 0;
	nshell = 0;
	low = opac.ipElement[ipHYDROGEN][ion][nshell][0];
	ipHi = opac.ipElement[ipHYDROGEN][ion][nshell][1];
	ipop = opac.ipElement[ipHYDROGEN][ion][nshell][2];
	/* OpacityAdd1Subshell will not do anything if static opacities do not need to be reset*/
	OpacityAdd1Subshell(ipop,low,ipHi,hmi.H2_total*2.f , 'v' );

	/* include CO */
	dense.xIonDense[ipOXYGEN][0] += co.hevmol[ipCO];
	dense.xIonDense[ipCARBON][0] += co.hevmol[ipCO];

	/* main opacity add */
	for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
	{
		/* this element may be turned off */
		if( dense.lgElmtOn[nelem] )
		{ 
			OpacityAdd1Element(nelem);
		}
	}

	/* now reset the abundances */
	dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
	dense.xIonDense[ipCARBON][0] = SaveCarbon1;
	/*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s] = SaveHydrogen1;*/

	/* following are opacities due to specific excited levels */

	/* nitrogen opacity
	 * excited level of N+ */
	OpacityAdd1Subshell(opac.in1[2],opac.in1[0],opac.in1[1],
		dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' );

	/* oxygen opacity
	 * excited level of Oo */
	OpacityAdd1Subshell(opac.ipo1exc[2],opac.ipo1exc[0],opac.ipo1exc[1],
	  dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v');

	/* O2+ excited states */
	OpacityAdd1Subshell(opac.ipo3exc[2],opac.ipo3exc[0],opac.ipo3exc[1],
	  dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v');

	OpacityAdd1Subshell(opac.ipo3exc3[2],opac.ipo3exc3[0],opac.ipo3exc3[1],
	  dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v');

	/* magnesium opacity
	 * excited level of Mg+ */
	OpacityAdd1Subshell(opac.ipOpMgEx,opac.ipmgex,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
		dense.xIonDense[ipMAGNESIUM][1]* atoms.popmg2,'v');

	/* calcium opacity
	 * photoionization of excited levels of Ca+ */
	OpacityAdd1Subshell(opac.ica2op,opac.ica2ex[0],opac.ica2ex[1],
	  ca.popca2ex,'v');

	/*******************************************************************
	 *
	 * complete evaluation of total opacity by adding in the static part and grains
	 *
	 *******************************************************************/

	/* this loop defines the variable iso.ConOpacRatio[ipH_LIKE][nelem][n],
	 * the ratio of not H to Hydrogen opacity.  for grain free environments
	 * at low densities this is nearly zero.  The correction includes 
	 * stimulated emission correction */
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			/* this element may be turned off */
			if( dense.lgElmtOn[nelem] )
			{ 
				/* this branch is for startup only */
				if( nzone < 1 )
				{
					for( n=0; n < iso.numLevels[ipISO][nelem]; n++ )
					{
						if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
						{
							if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )
							{
#if RJRWLOOP
								/* >>chng 02 may 06, use general form of threshold cs */
								/*double t1 = atmdat.STH[n]/(POW2(nelem+1.-ipISO));*/
								long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
								double t2 = csphot(
									ip ,
									ip ,
									iso.ipOpac[ipISO][nelem][n] );

								iso.ConOpacRatio[ipISO][nelem][n] = 1.f-(float)(
									(iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/
									opac.opacity_abs[ip-1]);
#else
								/* >>chng 02 may 06, use general form of threshold cs */
								/*double t1 = atmdat.STH[n]/(POW2(nelem+1.-ipISO));*/
								double t2 = csphot(
									iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
									iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
									iso.ipOpac[ipISO][nelem][n] );

								iso.ConOpacRatio[ipISO][nelem][n] = (float)(
									(opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] - 
									iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/
									opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1]);
#endif
							}
							else
							{
								iso.ConOpacRatio[ipISO][nelem][n] = 0.;
							}
						}
					}
				}
				/* end branch for startup only, start branch for all zones including startup */
				for( n=0; n < iso.numLevels[ipISO][nelem]; n++ )
				{
					/* ratios of other to total opacity for continua of all atoms done with iso model */
					if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
					{
						if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
						{
							/* first get departure coef */
							if( iso.DepartCoef[ipISO][nelem][n] > 1e-30 && (!conv.lgSearch ) )
							{
								/* this is the usual case, use inverse of departure coef */
								fac = 1./iso.DepartCoef[ipISO][nelem][n];
							}
							else if( conv.lgSearch )
							{
								/* do not make correction for stim emission during search
								* for initial temperature solution, since trys are very non-equil */
								fac = 0.;
							}
							else
							{
								fac = 1.;
							}

							/* now get opaicty ratio with correction for stimulated emission */
							if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )
							{
								/* >>chng 02 may 06, use general form of threshold cs */
#if RJRWMACRO
								long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];

								double t2 = csphot(
									ip ,
									ip ,
									iso.ipOpac[ipISO][nelem][n] );

								double opacity_this_species = 
									iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2*
									(1. - fac*rfield.ContBoltz[ip-1]);

								double opacity_fraction =  1. - opacity_this_species / opac.opacity_abs[ip-1];
								if (opacity_fraction < 0)
									opacity_fraction = 0.;
#else
								double t2 = csphot(
									iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
									iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
									iso.ipOpac[ipISO][nelem][n] );

								int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;

								double opacity_this_species = 
									iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2*
									(1. - fac*rfield.ContBoltz[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1]);

								double opacity_fraction = 
									/* >>chng 02 nov 24, add max2 */
									MAX2( 0. , opac.opacity_abs[ip] - opacity_this_species ) / opac.opacity_abs[ip];
#endif
								/* use mean of old and new ratios */
								iso.ConOpacRatio[ipISO][nelem][n] = (float)(
									iso.ConOpacRatio[ipISO][nelem][n]* 0.75 + 0.25*opacity_fraction );

								/*if( nzone==200 && ipISO==ipHE_LIKE &&
									nelem==ipHELIUM && n==0 )
								{
									fprintf(ioQQQ,"ssss\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
										opacity_this_species,
										opacity_fraction,
										opac.opacity_abs[ip],
										iso.ConOpacRatio[ipISO][nelem][n],
										iso.Pop2Ion[ipISO][nelem][n],
										dense.xIonDense[nelem][nelem+2-ipISO-1]);
								}*/
#if RJRWMACRO
								if (iso.ConOpacRatio[ipISO][nelem][n] < 0.)
									iso.ConOpacRatio[ipISO][nelem][n] = 0.;
#endif
							}
							else
							{
								iso.ConOpacRatio[ipISO][nelem][n] = 0.;
							}
#if (!RJRWMACRO)
							iso.ConOpacRatio[ipISO][nelem][n] = (float)MAX2(0.,iso.ConOpacRatio[ipISO][nelem][n]);
#endif
						}
						else
						{
							iso.ConOpacRatio[ipISO][nelem][n] = 0.;
						}
					}
					else
					{
						iso.ConOpacRatio[ipISO][nelem][n] = 0.;
					}
				}
			}
		}
	}

	/* add dust grain opacity if dust present */
	if( gv.lgDustOn )
	{
		/* generate current grain opacities since may be function of depth */
		/* >>chng 01 may 11, removed code to update grain opacities, already done by GrainChargeTemp */
		for( i=0; i < rfield.nflux; i++ )
		{
			/* units cm-1 */
			opac.opacity_sct[i] += gv.dstsc[i]*dense.gas_phase[ipHYDROGEN];
			opac.opacity_abs[i] += gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
		}
	}

	/* check that opacity is sane */
	for( i=0; i < rfield.nflux; i++ )
	{
		/* OpacStatic was zeroed in OpacityZero, incremented in opacityadd1subshell */
		opac.opacity_abs[i] += opac.OpacStatic[i];
		/* make sure that opacity is positive */
		ASSERT( opac.opacity_abs[i] > 0. );
	}

	/* this is hack (option) to hose opacities for 2004 leiden pdr meeting */
	if( gv.lgLeidenGrnOpacHack )
	{
		/* try to do the VERY simple physics included in other pdr models,
		 * clobber our beautiful, realistic treatment with the c%*p */
		/*for( i=rfield.ipG0_TH85_lo; i < rfield.nflux; i++ )*/
		for( i=0; i < rfield.nflux; i++ )
		{
			/* >>chng 04 may 21, assume albedo really is 0.5 */
			opac.opacity_sct[i] = 1.51*6.289e-22 * dense.gas_phase[ipHYDROGEN];
			opac.opacity_abs[i] = 1.51*6.289e-22 * dense.gas_phase[ipHYDROGEN];
		}
	}

	/* compute gas albedo here */
	for( i=0; i < rfield.nflux; i++ )
	{
		opac.albedo[i] = opac.opacity_sct[i]/
			(opac.opacity_sct[i] + opac.opacity_abs[i]);
	}

	/* during search phase set old opacity array to current value */
	if( conv.lgSearch )
		OpacityZeroOld();

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     OpacityAddTotal returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e HeI,II;%10.2e\n", 
		  iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH1s], 
		  iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2p], 
		  iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][3], 
		  iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][4], 
		  iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][ipHe1s1S] );
	}
	{
		/* following should be set true to print strongest ots contributors */
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && (nzone>=378)/**/ )
		{
			if( nzone > 380 ) 
				exit(1);
			for( i=0; i<rfield.nflux; ++i )
			{
				fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
					conv.nPres2Ioniz,
					rfield.anu[i],
					opac.opacity_abs[i],
					rfield.otscon[i],
					rfield.otslin[i]);
			}
		}
	}

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

