/*AddOpac derive total opacity for this position,
 * called by by */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "nhe3lvl.h"
#include "nhe1lvl.h"
#include "grainvar.h"
#include "opacpoint.h"
#include "mgexc.h"
#include "nitexc.h"
#include "gaunts.h"
#include "opac.h"
#include "rfield.h"
#include "o3exc.h"
#include "hemase.h"
#include "converge.h"
#include "ircoil.h"
#include "trace.h"
#include "tepowers.h"
#include "elmton.h"
#include "hmi.h"
#include "phycon.h"
#include "he.h"
#include "phe1lv.h"
#include "iphmin.h"
#include "h.h"
#include "ffsum.h"
#include "hmidep.h"
#include "heopfr.h"
#include "ophf.h"
#include "stimax.h"
#include "he1bn.h"
#include "he3bn.h"
#include "nh.h"
#include "nhe1.h"
#include "he3n.h"
#include "zonecnt.h"
#include "pca2ex.h"
#include "ionfracs.h"
#include "hydrophocs.h"
#include "opacz.h"
#include "opacin.h"
#include "opinb.h"
#include "putopacity.h"
#include "addopac.h"

void AddOpac(void)
{
	long int i, 
	  ipZ, 
	  limit, 
	  n, 
	  nd;
	double depcoe, 
	  DepartCoefInv ,
	  fac, 
	  fac1, 
	  fachmi, 
	  fachmi1, 
	  sum;
	static float oldhe3;
	double dfactor ;

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

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

	/* free electron scattering opacity, Compton recoil energy loss */
	for( i=0; i < (ircoil.ipRecoilIon - 1); i++ )
	{
		/* scattering part of total opacity */
		opac.scatop[i] += opac.OpacStack[i-1+OpacPoint.iopcom]*
		  phycon.eden;
	}

	for( i=ircoil.ipRecoilIon-1; i < rfield.nflux; i++ )
	{
		/* scattering part of total opacity
		 * at high energies Compton recoil is important, simply loose photons
		 * since code not designed for Compton thick clouds */
		fac = rfield.anu[i]/7.512e4;
		fac /= 1. + fac;
		opac.scatop[i] += (opac.OpacStack[i-1+OpacPoint.iopcom]*
		  phycon.eden*(1. - fac));
		opac.opac[i] += (opac.OpacStack[i-1+OpacPoint.iopcom]*
		  phycon.eden*fac);
	}

	for( i=ircoil.ipRecoilIon-1; i < rfield.nflux; i++ )
	{
		/* add in bound hydrogen electron scattering, treated as absorption */
		opac.opac[i] += opac.OpacStack[i-1+OpacPoint.iopcom]*
		  h.hi;
	}

	/* bound helium electron scattering, treated as absorption */
	for( i=ircoil.irche-1; i < rfield.nflux; i++ )
	{
		opac.opac[i] += (opac.OpacStack[i-1+OpacPoint.iopcom]*
		  he.hei*2.);
	}

	/* opacity due to pair production */
	sum = phycon.hden + 4.*he.ahe;
	opacin(OpacPoint.ioppr,OpacPoint.ippr,rfield.nflux,(float)sum,'s');

	/* hydrogen, helium, heavy element brems (free-free) opacity,
	 * assuming Hydrogen free-free gaunt factors
	 * logic change June 4 95 to better treat h- absorption
	 * FreeFreeOpacity is evaluated in ff heating routine, then added here
	 * OPSV is missing factor of 1E-20 to avoid underflow */
	fac = (phycon.eden/1e10)*((h.hii + he.heii + 4.*he.heiii + ffsum.EdenFFSum)/
	  1e10)/TePowers.sqrte;
	fac1 = fac*TE1RYD/phycon.te;
	fachmi = (phycon.eden/1e10)*(h.hii*hydro.hn[0][IP1S]/1e10)/TePowers.sqrte;
	fachmi1 = fachmi*TE1RYD/phycon.te;

	for( i=0; i < rfield.nflux; i++ )
	{
		if( rfield.ContBoltz[i] < 0.995 )
		{
			/* use 1-exp(hn/kT)) only if exp is small
			 * last term is scaled h minus brems free-free absorption */
			dfactor = opac.OpacStack[i-1+OpacPoint.ipBrems]*(1. - 
			  rfield.ContBoltz[i])*gaunts.gff[i];
			opac.FreeFreeOpacity[i] = dfactor*fac;
			/* this is the H- part, opac here is ratio of H- to ff */
			opac.FreeFreeOpacity[i] += dfactor*fachmi*opac.OpacStack[i-1+OpacPoint.iphmra];
		}
		else
		{
			dfactor = opac.OpacStack[i-1+OpacPoint.ipBrems]*rfield.anu[i]*
			  gaunts.gff[i];
			opac.FreeFreeOpacity[i] = dfactor*fac1;
			/* h- free free absorption */
			opac.FreeFreeOpacity[i] += dfactor*fachmi1*opac.OpacStack[i-1+OpacPoint.iphmra];
		}
		opac.opac[i] += opac.FreeFreeOpacity[i];
	}

	/* H minus absorption, with correction for stimulated emission */
	limit = MIN2(rfield.nflux,nhe1Com.nhe1[0]);
	if( hmidepCom.hmidep > SMALLFLOAT )
	{
		DepartCoefInv = 1./hmidepCom.hmidep;
	}
	else
	{
		/* the hmidep departure coef can become vastly small in totally
		 * neutral gas (no electrons) */
		DepartCoefInv = 1.;
	}

	for( i=iphminCom.iphmin-1; i < limit; i++ )
	{
		opac.opac[i] += (opac.OpacStack[i-iphminCom.iphmin+OpacPoint.iphmop]*
		  hmi.hminus*(1. - rfield.ContBoltz[i]*DepartCoefInv));
	}

	/* H2P h2plus bound free opacity */
	limit = MIN2(rfield.nflux,OpacPoint.ih2pnt[1]);
	for( i=OpacPoint.ih2pnt[0]-1; i < limit; i++ )
	{
		opac.opac[i] += hmi.h2plus*opac.OpacStack[i-OpacPoint.ih2pnt[0]+
		  OpacPoint.ih2pof];
	}

	/* get total population of hydrogen ground */
	if( xIonFracs[2][0] <= 0. )
	{
		fac = xIonFracs[1][0];
	}
	else
	{
		fac = xIonFracs[2][0]*hydro.hn[0][IP1S];
	}

	/* Ly a damp wing opac (Rayleigh scattering) */
	limit = MIN2(rfield.nflux, nh.ipHn[0][IP1S]);
	for( i=0; i < limit; i++ )
	{
		opac.scatop[i] += (fac*opac.OpacStack[i-1+OpacPoint.ipRayScat]);
	}

	 /* remember largest correction for stimulated emission */
	if( hydro.hbn[0][IP1S] > 1e-30 && !conv.lgSearch )
	{
		stimaxCom.stimax[0] = 
			(float)(MAX2(stimaxCom.stimax[0],rfield.ContBoltz[nh.ipHn[0][IP1S]-1]/
			(float)(hydro.hbn[0][IP1S])));
	}

	if( hydro.hbn[0][IP2P] > 1e-30 && !conv.lgSearch )
	{
		stimaxCom.stimax[1] = 
			(float)(MAX2(stimaxCom.stimax[1],rfield.ContBoltz[nh.ipHn[0][2]-1]/
			(float)(hydro.hbn[0][IP2P])));
	}

	/* generate current grain opacities since may be function of depth */
	if( GrainVar.lgDustOn )
	{
		for( i=0; i < rfield.nupper; i++ )
		{
			/* these are total absorption and scattering cross sections */
			GrainVar.dstab[i] = 0.;
			GrainVar.dstsc[i] = 0.;
		}

		/* grain abundance may be a function of depth */
		for( nd=0; nd < NDUST; nd++ )
		{
			if( GrainVar.lgDustOn1[nd] )
			{
				for( i=0; i < rfield.nupper; i++ )
				{
					/* these are total absorption and scattering cross sections */
					GrainVar.dstab[i] += GrainVar.dstab1[nd][i]*
					  GrainVar.dstAbund[nd];
					/* this must be positive, zero in case of uncontrolled underflow */
					assert( GrainVar.dstab[i]> 0.);
					GrainVar.dstsc[i] += GrainVar.dstsc1[nd][i]*
					  GrainVar.dstAbund[nd];
				}
			}
		}

		/* add dust grain opacity if dust present */
		for( i=0; i < rfield.nflux; i++ )
		{
			opac.scatop[i] += GrainVar.dstsc[i]*phycon.hden;
			opac.opac[i] += GrainVar.dstab[i]*phycon.hden;
		}
	}

	/* check whether hydrogen or Helium singlets mased, if not in search mode */
	if( !conv.lgSearch )
	{
		if( hydro.hbn[1][IP1S] > 0. )
		{
			if( rfield.ContBoltz[nh.ipHn[1][IP1S]-1]/hydro.hbn[1][IP1S] > 1. )
				hemase.lgHeMase = TRUE;
		}
		if( he1bnCOM.he1bn[0] > 0. )
		{
			if( rfield.ContBoltz[nhe1Com.nhe1[0]-1]/he1bnCOM.he1bn[0] > 1. )
				hemase.lgHeMase = TRUE;
		}
	}

	/* helium triplets hei */
	if( conv.lgSearch )
	{
		/* neglect stimulated emission corr if search underway
		 * >>chng 96 feb 15, was 1 */
		depcoe = 0.;
	}
	else
	{
		depcoe = he3bnCom.he3bn[0];
	}

	if( ZoneCnt.nzone <= 1 )
	{
		oldhe3 = he.heii*he3nCom.he3n[0];
	}
	oldhe3 = (float)(0.75*oldhe3 + 0.25*he.heii*he3nCom.he3n[0]);

	/* add helium triplet opacity, this version of OPACIN has stimulated emission */
	opinb(OpacPoint.ioptri,he.nhei3,nhe1Com.nhe1[0],oldhe3,depcoe, 'v' );

	/* helium singlets excited states hei, 
	 * these opacities only extend up to 1.8 ryd 
	 * excited states ARE NOT picked up in PutOpacity */
	for( n=1; n < 4; n++ )
	{
		if( conv.lgSearch )
		{
			/* neglect stimulated emission corr if search underway
			 * >>chng 96 feb 15, was 1 */
			depcoe = 0.;
		}
		else
		{
			depcoe = he1bnCOM.he1bn[n];
		}
		opinb(OpacPoint.iophe1[n],nhe1Com.nhe1[n],nhe1Com.nhe1[0],
		  phe1lv.he1n[n]*he.heii,depcoe , 'v');
	}
	for( n=4; n < NHE1LVL; n++ )
	{
		if( conv.lgSearch )
		{
			/* neglect stimulated emission corr if search underway
			 * >>chng 96 feb 15, was 1 */
			depcoe = 0.;
		}
		else
		{
			depcoe = he1bnCOM.he1bn[n];
		}
		opinb(OpacPoint.iophe1[n],nhe1Com.nhe1[n],nhe1Com.nhe1[0],
		  phe1lv.he1n[n]*he.heii,depcoe , 's');
	}

	/* following loop adds standard opacities for first 30 elements 
	 * most heavy element opacity added here */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		/* this element may be turned off */
		if( elmton.lgElmtOn[ipZ] )
		{ 
			PutOpacity(ipZ);
		}
	}

	/* following are opacities due to specific levels */

	/* nitrogen opacity
	 * excited level of N+ */
	opacin(OpacPoint.in1[2],OpacPoint.in1[0],OpacPoint.in1[1],xIonFracs[1][6]*
	  nitexc.p2nit,'v');

	/* oxygen opacity
	 * excited level of Oo */
	opacin(OpacPoint.ipo1exc[2],OpacPoint.ipo1exc[0],OpacPoint.ipo1exc[1],
	  xIonFracs[1][7]*o3exc.poiexc,'v');

	/* O2+ excited states */
	opacin(OpacPoint.ipo3exc[2],OpacPoint.ipo3exc[0],OpacPoint.ipo3exc[1],
	  xIonFracs[3][7]*o3exc.poiii2,'v');

	opacin(OpacPoint.ipo3exc3[2],OpacPoint.ipo3exc3[0],OpacPoint.ipo3exc3[1],
	  xIonFracs[3][7]*o3exc.poiii3,'v');

	/* magnesium opacity
	 * excited level of Mg+ */
	opacin(OpacPoint.ipOpMgEx,OpacPoint.ipmgex,nh.ipHn[0][IP1S],
		xIonFracs[2][11]* mgexc.popmg2,'v');

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

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

	for( i=0; i < rfield.nflux; i++ )
	{
		/* OpacStatic was zeroed in opacz */
		opac.opac[i] += opac.OpacStatic[i];
		/* make sure that opacity is positive */
		assert( opac.opac[i] > 0. );
	}

	/* total opacity now defined - now look at some ratios */
	for( n=0; n < NHE1LVL; n++ )
	{
		if( nhe1Com.nhe1[n] < rfield.nflux && 
			opac.opac[nhe1Com.nhe1[n]-1] > 0. )
		{
			if( he1bnCOM.he1bn[n] > 1e-30 && (!conv.lgSearch) )
			{
				fac = 1./he1bnCOM.he1bn[n];
			}
			else if( conv.lgSearch )
			{
				fac = 0.;
			}
			else
			{
				fac = 1.;
			}

			heopfr.ophe1f[n] = (float)(((opac.opac[nhe1Com.nhe1[n]-1] - 
			  he.heii*phe1lv.he1n[n]*opac.OpacStack[OpacPoint.iophe1[n]-
			  1]*(1. - fac*rfield.ContBoltz[nhe1Com.nhe1[n]-1]))/
			  opac.opac[nhe1Com.nhe1[n]-1] + heopfr.ophe1f[n])/
			  2.);
			heopfr.ophe1f[n] = (float)MAX2(0.,heopfr.ophe1f[n]);
		}

		else
		{
			heopfr.ophe1f[n] = 0.;
		}
	}


	/* this loop defines the variable ophf.HOpacRatio[ipZ][n],
	 * the ratio of not H to Hydrogen opaity.  for grain free environments
	 * at low densities this is nearly zero.  The correction includes 
	 * stimulated emission correction */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		/* this element may be turned off */
		if( elmton.lgElmtOn[ipZ] )
		{ 
			if( ZoneCnt.nzone < 1 )
			{
				for( n=IP1S; n <= nhlevel; n++ )
				{
					if(nh.ipHn[ipZ][n] < rfield.nflux )
					{
						if( opac.opac[nh.ipHn[ipZ][n]-1] > 1e-30 )
						{
							ophf.HOpacRatio[ipZ][n] = (float)((opac.opac[nh.ipHn[ipZ][n]-1] - 
								hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ]*
							  HydroPhoCS.STH[n]/(POW2(ipZ+1.)))/opac.opac[nh.ipHn[ipZ][n]-1]);
						}
						else
						{
							ophf.HOpacRatio[ipZ][n] = 0.;
						}
					}
				}
			}

			for( n=IP1S; n <= nhlevel; n++ )
			{
				/* ratios of other to total opacity for H continua
				 * same for Lyman, Balmer continua of hydrogen */
				if(nh.ipHn[ipZ][n] < rfield.nflux )
				{
					if( opac.opac[nh.ipHn[ipZ][n]-1] > 0. )
					{
						/* first get departure coef */
						if( hydro.hbn[ipZ][n] > 1e-30 && (!conv.lgSearch ) )
						{
							/* this is the usual case, use inverse of departure coef */
							fac = 1./hydro.hbn[ipZ][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.opac[nh.ipHn[ipZ][n]-1] > 1e-30 )
						{
							ophf.HOpacRatio[ipZ][n] = (float)(
								ophf.HOpacRatio[ipZ][n]* 0.75 + 
								0.25*(
								opac.opac[nh.ipHn[ipZ][n]-1] - 
								hydro.hn[ipZ][n]*xIonFracs[ipZ+2][ipZ]*
								HydroPhoCS.STH[n]/(POW2(ipZ+1.))*
								(1. - fac*rfield.ContBoltz[nh.ipHn[ipZ][n]-1]))/
								opac.opac[nh.ipHn[ipZ][n]-1]);
						}
						else
						{
							ophf.HOpacRatio[ipZ][n] = 0.;
						}
						ophf.HOpacRatio[ipZ][n] = (float)MAX2(0.,ophf.HOpacRatio[ipZ][n]);
					}
					else
					{
						ophf.HOpacRatio[ipZ][n] = 0.;
					}
				}
				else
				{
					ophf.HOpacRatio[ipZ][n] = 0.;
				}
			}
		}
	}
	/*fprintf(ioQQQ," HOPac=%e\n", ophf.HOpacRatio[0][0]);*/

	/* >>chng 97 mar 8, zero out very small opacities, below precision on
	 * current cpu.  done to stop ots fluxes from exploding at low densities
	 * >>chng 97 may 6, decreased from 1e-30 to 1e-34, ldl.in did not work */

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

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     AddOpac returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e HeI,II;%10.2e\n", 
		  ophf.HOpacRatio[0][IP1S], ophf.HOpacRatio[0][IP2P], ophf.HOpacRatio[0][3], 
		  ophf.HOpacRatio[0][4], heopfr.ophe1f[0] );
	}

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