/*TauInc increment optical depths once per zone, called after RadInc */
#include <string.h>
#include "cddefines.h"
#define	CAPMASER	(-30.)
#include "physconst.h"
#include "taulines.h"
#include "nhe1lvl.h"
#include "nfeii.h"
#include "grainvar.h"
#include "hydrogenic.h"
#include "lyaext.h"
#include "densty.h"
#include "rfield.h"
#include "iphmin.h"
#include "showme.h"
#include "heat.h"
#include "hmidep.h"
#include "trace.h"
#include "he1tau.h"
#include "he1nxt.h"
#include "wind.h"
#include "ionrange.h"
#include "doppvel.h"
#include "he3lines.h"
#include "phe1lv.h"
#include "abrems.h"
#include "filfac.h"
#include "tcool.h"
#include "hehp.h"
#include "colden.h"
#include "hmi.h"
#include "he3tau.h"
#include "taucon.h"
#include "h.h"
#include "he.h"
#include "phycon.h"
#include "radius.h"
#include "colmas.h"
#include "jmin.h"
#include "fe2tau.h"
#include "logte.h"
#include "struc.h"
#include "zonecnt.h"
#include "ionfracs.h"
#include "he1stat.h"
#include "dtmase.h"
#include "velset.h"
#include "pop371.h"
#include "metprt.h"
#include "aver.h"
#include "linetauinc.h"
#include "molcol.h"
#include "mean.h"
#include "chlinelbl.h"
#include "tauinc.h"

void TauInc(void)
{

	char chLbl[11];

	long int i, 
	  ipHi, 
	  ipZ, 
	  j, 
	  ipLo, 
	  nd;

	double ajmass, 
	  dtau, 
	  factor, 
	  fnext, 
	  r, 
	  rjeans;

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

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " TauInc called.\n" );
	}

	/* following is general method to find means weighted by various functions
	 * called in StartIter to initialize to zero, called here to put in numbers
	 * results will be weighted by radius and volumn
	 * this is the only place things must be entered to create averages */
	aver("zone",1.,1.,"          ");
	aver("doit",phycon.te,1.,"    Te    ");
	aver("doit",phycon.te,phycon.eden,"  Te(Ne)  ");
	aver("doit",phycon.te,phycon.eden*h.hii," Te(NeNp) ");
	aver("doit",phycon.te,phycon.eden*he.heii," Te(NeHe+)");
	aver("doit",phycon.te,phycon.eden*he.heiii,"Te(NeHe2+)");
	aver("doit",phycon.te,phycon.eden*xIonFracs[2][7]," Te(NeO+) " );
	aver("doit",phycon.te,phycon.eden*xIonFracs[3][7]," Te(NeO2+)");
	aver("doit",phycon.hden,1.,"    NH    ");
	aver("doit",phycon.eden,xIonFracs[3][7],"  Ne(O2+) ");
	aver("doit",phycon.eden,h.hii,"   Ne(Np) ");

	/* save information about structure of model, now used to get t^2 */
	nd = MIN2(ZoneCnt.nzone , NZLIM)-1;
	struc.testr[nd] = phycon.te;
	struc.pdenstr[nd] = densty.pden;
	struc.heatstr[nd] = heat.htot;
	struc.coolstr[nd] = tcool.ctot;
	struc.volstr[nd] = (float)radius.dVeff;
	struc.radstr[nd] = (float)radius.dReff;
	struc.histr[nd] = h.hi;
	struc.hiistr[nd] = h.hii;
	struc.ednstr[nd] = (float)phycon.eden;
	struc.o3str[nd] = xIonFracs[3][7];
	struc.DenMass[nd] = densty.xMassDensity;

	abrems.dlnenp += phycon.eden*(double)(h.hii)*radius.dReff;

	/* counter for number of times Lya excit temp hotter than gas */
	if( LyaExT.TexcLya > phycon.te )
	{
		LyaExT.nLyaHot += 1;
		if( LyaExT.TexcLya > LyaExT.TLyaMax )
		{
			LyaExT.TLyaMax = LyaExT.TexcLya;
			LyaExT.TeLyaMax = phycon.te;
			LyaExT.nZTLaMax = ZoneCnt.nzone;
		}
	}

	coldenCom.colden[IPCOLUMNDENSITY-1] += (float)(phycon.hden*radius.dReff);
	coldenCom.colden[IPCHMIN-1] += hmi.hminus*(float)radius.dReff;
	coldenCom.colden[IPCOLH2-1] += hmi.htwo*(float)radius.dReff;
	coldenCom.colden[IPCHEHP-1] += hehpCom.hehp*(float)radius.dReff;
	coldenCom.colden[IPCH2PLS-1] += hmi.h2plus*(float)radius.dReff;
	coldenCom.colden[IPCHII-1] += h.hii*(float)radius.dReff;
	coldenCom.colden[IPCHI-1] += h.hi*(float)radius.dReff;

	/* now add total molecular column densities */
	molcol("ADD ");

	/* increment forming the mean ionizatin and temperature */
	MeanInc();

	/*-----------------------------------------------------------------------*/

	GrainVar.reftot += (float)radius.dReff;
	GrainVar.avhden += (float)radius.dReff;

	for( nd=0; nd < NDUST; nd++ )
	{
		if( GrainVar.lgDustOn1[nd] )
		{
			GrainVar.avdust[nd] += GrainVar.tedust[nd]*(float)radius.dReff;
			GrainVar.avdft[nd] += GrainVar.DustDftVel[nd]*(float)radius.dReff;
			GrainVar.avdpot[nd] += (float)(GrainVar.dstpot[nd]*EVRYD*radius.dReff);
		}
	}

	GrainVar.thtot += (float)(phycon.te*radius.dReff*phycon.hden);
	colmas.TotMassColl += densty.xMassDensity*(float)radius.dReff;
	colmas.tmas += phycon.te*densty.xMassDensity*(float)radius.dReff;
	colmas.wmas += densty.wmole*densty.xMassDensity*(float)radius.dReff;

	/* now find minimum Jeans length and mass; length in cm */
	rjeans = 7.79637 + (logte.alogte - log10(densty.wmole) - log10(densty.xMassDensity*
	  filfac.FillFac))/2.;

	/* minimum Jeans mass in gm */
	ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*densty.xMassDensity*
	  filfac.FillFac);

	/* now remember smallest */
	jmin.rjnmin = (float)MIN2(jmin.rjnmin,rjeans);
	jmin.ajmmin = (float)MIN2(jmin.ajmmin,ajmass);

	taucon.telec += (float)(radius.dReff*phycon.eden*6.65e-25);
	taucon.thmin += 
		(float)(radius.dReff*hmi.hminus*3.9e-17*(1. - rfield.ContBoltz[iphminCom.iphmin-1]/
	  hmidepCom.hmidep));

	/* find line widths for thermal and turbulent motion
	 * CLOUDY uses line center optical depths, =0.015 F / DELTA NU */
	velset();

	/* following used for wind solution in TAUCHN */
	if( wind.windv > 0. )
		wind.wndeff = 
		(float)(MIN2(radius.depth,radius.Radius)/MAX2(DoppVel.doppler[0],
		  wind.windv)*filfac.FillFac);

	/* prevent maser runaway */
	dtMase.dTauMase = 0.;

	/* if following is still zero when maser used, then maser was in
	 * heavy element line */
	dtMase.ijMase = 0;

	/* hydrogenic species */
	for( ipZ=0; ipZ < LIMELM; ipZ++ )
	{
		if( (IonRange.IonHigh[ipZ] == ipZ + 2)  )
		{
			/* must convert to physical scale for opacities,
			 * get ionic abundances first */
			if( xIonFracs[ipZ+2][ipZ] > 1e-30 )
			{
				factor = xIonFracs[ipZ+2][ipZ];
			}
			else
			{
				/* case where almost no parent ion - this will make
				 * very large line opacity, so background dest small */
				factor = 1.;
			}

			for( ipLo=0; ipLo <= (nhlevel - 1); ipLo++ )
			{
				for( ipHi=ipLo + 1; ipHi <= nhlevel; ipHi++ )
				{
					/* must temporarily make ipLnPopOpc physical */
					HydroLines[ipZ][ipHi][ipLo].PopOpc *= 
						(float)factor;
					/* actually do the work */
					LineTauInc(&HydroLines[ipZ][ipHi][ipLo]);
					/* go back to original units so that final correction ok */
					HydroLines[ipZ][ipHi][ipLo].PopOpc /= 
						(float)factor;
				}
			}
		}
	}

	/* increment optical depths for all heavy element lines
	 * same routine does wind and static */
	for( i=1; i <= nLevel1; i++ )
	{
		LineTauInc(&TauLines[i]);
		if( TauLines[i].TauIn < CAPMASER )
		{
			fprintf( ioQQQ, " TauInc, Tau line maser tau too small,%10.2e\n", 
			  TauLines[i].TauIn );
			strcpy( chLbl, chLineLbl(&TauLines[i] ) );
			fprintf( ioQQQ, " line=%3ld zone=%4ld %10.10s\n", 
			  i, ZoneCnt.nzone, chLbl );
			ShowMe();
			puts( "[Stop in tauinc]" );
			exit(1);
		}
	}

	/* all lines in cooling with g-bar */
	for( i=0; i < nWindLine; i++ )
	{
		LineTauInc(&TauLine2[i]);
		if( TauLine2[i].TauIn < CAPMASER )
		{
			fprintf( ioQQQ, " TauInc, Wind line maser tau too small,%10.2e\n", 
			  TauLine2[i].TauIn );
			fprintf( ioQQQ, " line=%3ld zone=%4ld\n", i, ZoneCnt.nzone );
			ShowMe();
			puts( "[Stop in tauinc]" );
			exit(1);
		}
	}

	/* do large FeII atom if this is enabled */
	if( FeII.lgFeIION )
	{
		FeIITauInc();
	}

	/* following is for static atmosphere */
	if( wind.windv == 0. )
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " TauInc updates for static atmos\n" );
		}

		/* He I lines
		 * Helium lines, optical depths are total, not He alone
		 * hydrogen and helium were added in with correction for stimulated emission */
		factor = he.heii*radius.dReff/DoppVel.doppler[1];
		fnext = he.heii*radius.dRNeff/DoppVel.doppler[1]*0.6;
		for( ipLo=0; ipLo < (NHE1LVL - 1); ipLo++ )
		{
			for( ipHi=ipLo + 1; ipHi < NHE1LVL; ipHi++ )
			{
				/* line absorption, corrected for stimulated emission */
				if( ipHi+1 < NHE1LVL - 2 )
				{
					dtau = (phe1lv.he1n[ipLo] - phe1lv.he1n[ipHi]*
					  he1statCOM.he1stat[ipLo]/he1statCOM.he1stat[ipHi])*
					  he1tauCOM.he1opc[ipLo][ipHi];
				}
				else
				{
					/* do not allow stim emission for psuedo-states */
					dtau = phe1lv.he1n[ipLo]*he1tauCOM.he1opc[ipLo][ipHi];
				}
				he1tauCOM.he1tau[ipLo][ipHi] += (float)(dtau*factor);
				he1nxtCOM.he1nxt[ipLo][ipHi] = (float)(he1tauCOM.he1tau[ipLo][ipHi] + 
				  dtau*fnext);
				if( he1nxtCOM.he1nxt[ipLo][ipHi] < CAPMASER )
				{
					fprintf( ioQQQ, " TauInc, He1 line maser tau too small,%10.2e\n", 
					  he1nxtCOM.he1nxt[ipLo][ipHi] );
					fprintf( ioQQQ, " upper, lower=%3ld%3ld zone=%4ld\n", 
					  ipHi+1, ipLo+1, ZoneCnt.nzone );
					ShowMe();
					puts( "[Stop in tauinc]" );
					exit(1);
				}
			}
		}

		t206.TauIn += (float)((phe1lv.he12s - phe1lv.he12p*
		  0.333)*1.162e-6*factor);
		if( t206.TauIn < CAPMASER )
		{
			fprintf( ioQQQ, " TauInc, He1 2.06 mic line maser tau too small,%10.2e\n", 
			  t206.TauIn );
			fprintf( ioQQQ, " zone=%4ld\n", ZoneCnt.nzone );
			ShowMe();
			puts( "[Stop in tauinc]" );
			exit(1);
		}

		factor = radius.dReff/DoppVel.doppler[1];
		he3tau[IPT10830-1].TauIn += (float)(8.80e-7*(he3lines.p2s - 
		  he3lines.p2p/3.)*factor);
		he3tau[IPT3889-1].TauIn += (float)(3.75e-8*he3lines.p2s*
		  factor);
		he3tau[IPT5876-1].TauIn += (float)(5.35e-7*he3lines.p2p*
		  factor);
		he3tau[IPT7065-1].TauIn += (float)(2.95e-7*he3lines.p2p*
		  factor);

		/* iron fe feii fe2 
		 * first are overlapping feii lines */
		for( i=0; i < NFEII; i++ )
		{
			fe2tau.Fe2TauLte[i] += fe2tau.feopc[i]*(float)radius.dReff;
		}

	}

	/* this is windy solution*/
	else
	{
		if( trace.lgTrace )
		{
			fprintf( ioQQQ, " TauInc updates for wind\n" );
		}

		/* expanding atmosphere
		 * min of depth or radius is to stop (nearly) infinite line
		 * line optical depths when ionization parameter set so
		 * radius itself is very large */
		r = MIN2(radius.depth,radius.Radius)/MAX2(DoppVel.doppler[0],
		  wind.windv)*filfac.FillFac;

		/* He I lines
		 * Helium lines and continua, optical depths are total, not He alone
		 * hydrogen and helium were added in with correction for stimulated emission */
		factor = he.heii*r;
		for( i=0; i < (NHE1LVL - 1); i++ )
		{
			for( j=i + 1; j < NHE1LVL; j++ )
			{
				/* line absorption, corrected for stimulated emission */
				if( j+1 <= NHE1LVL - 2 )
				{
					dtau = (phe1lv.he1n[i] - phe1lv.he1n[j]*he1statCOM.he1stat[i]/
					  he1statCOM.he1stat[j])*he1tauCOM.he1opc[i][j];
				}
				else
				{
					/* do not include stim emission for psuedo-states */
					dtau = phe1lv.he1n[i]*he1tauCOM.he1opc[i][j];
				}
				he1tauCOM.he1tau[i][j] = (float)(dtau*factor);
				he1nxtCOM.he1nxt[i][j] = he1tauCOM.he1tau[i][j];
			}
		}
		t206.TauIn = (float)(phe1lv.he12s*he.heii*1.162e-6*
		  r);
		t206.TauTot = t206.TauIn;

		he3tau[IPT10830-1].TauIn = (float)(r*(he3lines.p2s - 
		  he3lines.p2p/3.)*8.80e-7);
		he3tau[IPT10830-1].TauTot = he3tau[IPT10830-1].TauIn;
		he3tau[IPT3889-1].TauIn = (float)(3.75e-8*he3lines.p2s* r);
		he3tau[IPT3889-1].TauTot = he3tau[IPT3889-1].TauIn;
		he3tau[IPT5876-1].TauIn = (float)(5.35e-7*he3lines.p2p* r);
		he3tau[IPT5876-1].TauTot = he3tau[IPT5876-1].TauIn;
		he3tau[IPT7065-1].TauIn = (float)(2.95e-7*he3lines.p2p* r);
		he3tau[IPT7065-1].TauTot = he3tau[IPT7065-1].TauIn;

	}

	if( trace.lgTrace && trace.lgOptcBug )
	{
		metprt();
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " TauInc returns.\n" );
	}

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