/*TotalPressure determine the gas and line radiation pressures for current conditions */
/*RamPressure determine the ram pressure for current conditions */
#include <string.h>
#include <limits.h>
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "showme.h"
#include "hydrogenic.h"
#include "itercnt.h"
#include "powi.h"
#include "wind.h"
#include "phycon.h"
#include "hmi.h"
#include "densty.h"
#include "hevmolec.h"
#include "ionfracs.h"
#include "trace.h"
#include "atomcwgt.h"
#include "hwidth.h"
#include "tepowers.h"
#include "doppvel.h"
#include "widla.h"
#include "widthline.h"
#include "texcline.h"
#include "lyaext.h"
#include "ionrange.h"
#include "chlinelbl.h"
#include "pressure.h"

double TotalPressure(void)
{
	long int i, 
		/* will be used to point to line causing radiation pressure */
	  ip1=LONG_MAX,ip2=LONG_MAX,ip3=LONG_MAX, 
	  ipHi, 
	  ipLo, 
	  ipSet , /* used in debug print statement to see which line adds most pressure */
	  ipZ;
	double HHeLines, 
	  RadPres1, 
	  TotalPressure_v, 
	  p2, 
	  pmx;

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

	/* function TotalPressure - evaluate total pressure
	 * this function has a single dummy argument
	 * this function evaluates the following
	 * TotalPressure - gas and radiation pressure, dynes/cm^2
	 * GasPres gas pressure, same units
	 *
	 * HDEN is defined as total hydrogen nuclei densty; all forms
	 * present abarb with carbon monoxide
	 *
	 * densty.TotalNuclei is total number of nuclei and molecules
	 * nelem,0 is abundance of nelem */
	densty.TotalNuclei = 0.;
	for( i=0; i < LIMELM; i++ )
	{
		densty.TotalNuclei += xIonFracs[0][i];
	}

	/* add on molecules */
	for( i=0; i < NHEVML; i++ )
	{
		densty.TotalNuclei += hevmolec.hevmol[i];
	}

	/* do not double count C+, which is in this network as ipCP */
	densty.TotalNuclei -= hevmolec.hevmol[IPCP-1];

	/* add on hydrogen molecules */
	densty.TotalNuclei += hmi.htwo + hmi.h2plus + hmi.hminus + hmi.h3plus;

	/* phycon.ElecFrac is electron fraction, used for secondary ionization efficiency */
	phycon.ElecFrac = (float)(phycon.eden/densty.TotalNuclei);

	/* pden is the total number of particles per unit vol */
	densty.pden = (float)(phycon.eden + densty.TotalNuclei);

	/* this is total neutral particle densty for collisional ionization */
	phycon.CollidDensity = 0.;
	for( i=0; i < LIMELM; i++ )
	{
		phycon.CollidDensity += xIonFracs[1][i];
	}

	/* now add all the heavy molecules */
	for( i=0; i < NHEVML; i++ )
	{
		phycon.CollidDensity += hevmolec.hevmol[i];
	}

	/* do not double count C+, which is in this network as ipCP */
	phycon.CollidDensity -= hevmolec.hevmol[IPCP-1];

	phycon.CollidDensity += (float)(hmi.h2plus + hmi.hminus + hmi.h3plus + 
	  2.*hmi.htwo);

	/* gas pressure */
	pressure.GasPres = (float)(densty.pden*phycon.te*BOLTZMANN);

	/* WMOLE is molecular weight, AtomicMassUnit per particle */
	densty.wmole = 0.;
	for( i=0; i < LIMELM; i++ )
	{
		densty.wmole += xIonFracs[0][i]*AtomcWgt.AtomicWeight[i];
	}

	densty.wmole += (float)(hevmolec.hevmol[IPCH-1]*13. + hevmolec.hevmol[IPCHP-1]*
	  13. + hevmolec.hevmol[IPOH-1]*17. + hevmolec.hevmol[IPOHP-1]*
	  17. + hevmolec.hevmol[IPOTWO-1]*32. + hevmolec.hevmol[IPCTWO-1]*
	  24. + hevmolec.hevmol[IPCO-1]*28. + hevmolec.hevmol[IPCOP-1]*
	  28. + hevmolec.hevmol[IPH2O-1]*18. + hevmolec.hevmol[IPH2OP-1]*
	  18. + hevmolec.hevmol[IPO2P-1]*32. + hevmolec.hevmol[IPC2P-1]*
	  24. + hevmolec.hevmol[IPH3P-1]*3. + hevmolec.hevmol[IPH3OP-1]*
	  19. + hevmolec.hevmol[IPCH2P-1]*14. + hevmolec.hevmol[IPCH2-1]*
	  14. + hevmolec.hevmol[IPCH3-1]*15.);

	densty.wmole /= densty.pden;

	/* this is Lya excitation temperature, has always been evaluated here */
	LyaExT.TexcLya = (float)TexcLine( &HydroLines[0][IP2P][IP1S] );

	/* xMassDensity is grams per cc */
	densty.xMassDensity = (float)(densty.wmole*ATOMIC_MASS_UNIT*densty.pden);

	/* RadPres is pressure due to lines, lgRadPresON turns off or on */
	pressure.RadPres = 0.;
	pmx = 0.;
	/* find max rad pressure even if capped
	 * lgRadPrsON is turned off by NO RADIATION PRESSURE command */
	if( pressure.lgRadPrsON )
	{
		if( hydro.hn[0][IP1S] > 0. )
		{
			/* this is the 2p/1s level population, will be used for excitation
			 * temperature and radiation pressure */
			p2 = hydro.hn[0][IP2P]/hydro.hn[0][IP1S]/3.;

			/* this is the total Lya line width, in velocity units, cm/s */
			hwidth.HLineWidth[0] = (float)(
				widla(
				HydroLines[0][IP2P][IP1S].TauIn,
			   HydroLines[0][IP2P][IP1S].TauTot,
			   4.72e-2/TePowers.sqrte,
				DoppVel.doppler[0]));

			/* coefficient assumes that HLineWidth is half width, not total */
			pressure.plines[0] = (float)(p2*5.08e-6*hwidth.HLineWidth[0]);
		}
		else
		{
			p2 = 0.;
			hwidth.HLineWidth[0] = 0.;
			pressure.plines[0] = 0.;
		}
		/* plines=l-alph above
		 * plines[1] = hign n lyman ;   plines[2]=excit h lines in hlevel
		 * plines[3] = he i 10830
		 * plines[4]= HeII 304 plines[10]=HeI La (in helium) 
		 * plines[7]= HeI 584
		 *
		 * compute radiation pressure due to special H and He lines */
		HHeLines = 0.;
		ip1 = 0;
		ipSet = -1;
		for( i=0; i < 20; i++ )
		{
			if( pressure.plines[i] > pmx )
			{
				ipSet = 1;
				pmx = pressure.plines[i];
				ip1 = i;
			}
			HHeLines += pressure.plines[i];
		}

		pressure.RadPres = (float)HHeLines;

		/* line radiation pressure for hydrogen lines */
		for( ipZ=0; ipZ < LIMELM; ipZ++ )
		{
			if( IonRange.IonHigh[ipZ] == ipZ + 2 )
			{
				for( ipLo=IP1S; ipLo <= (nhlevel - 1); ipLo++ ) 
				{
					for( ipHi=MAX2( 3 , ipLo+1 ); ipHi <= nhlevel; ipHi++ )
					{
						if( HydroLines[ipZ][ipHi][ipLo].PopHi > SMALLFLOAT &&
							HydroLines[ipZ][ipHi][ipLo].PopOpc > SMALLFLOAT )
						{
							RadPres1 = 5.551e-2*
							  (powi(HydroLines[ipZ][ipHi][ipLo].EnergyWN/1.e6,4))*
							  (HydroLines[ipZ][ipHi][ipLo].PopHi/
							  HydroLines[ipZ][ipHi][ipLo].gHi)/
							  (HydroLines[ipZ][ipHi][ipLo].PopOpc/
							  HydroLines[ipZ][ipHi][ipLo].gLo)*
							  /* WidthLine gets linewidth in terms of RT effects */
							  WidthLine(&HydroLines[ipZ][ipHi][ipLo])* 2.;
							if( RadPres1 > pmx )
							{
								ipSet = 2;
								pmx = RadPres1;
								ip1 = ipHi;
								ip2 = ipLo;
								ip3 = ipZ;
							}
							pressure.RadPres += (float)RadPres1;
						}
					}
				}
			}
		}

		/* treat hydrogen's lya special, since it is special,
		 * but only do this if hydrogen is ionized */
		if( IonRange.IonHigh[0] == 2 )
		{
			/* hydrogen Lya by itself */
			if( HydroLines[0][IP2P][IP1S].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(HydroLines[0][IP2P][IP1S].EnergyWN/
				  1.e6,4))*(HydroLines[0][IP2P][IP1S].PopHi/
				  HydroLines[0][IP2P][IP1S].gHi)/
				  (HydroLines[0][IP2P][IP1S].PopLo/
				  HydroLines[0][IP2P][IP1S].gLo)*
				  WidthLine(&HydroLines[0][IP2P][IP1S])* 2.;
			}
			else
			{
				RadPres1 = 0.;
			}

			/* correct for counting Lya radiation pressure twice */
			pressure.RadPres = (float)MAX2(0.,pressure.RadPres-RadPres1);
		}

		/* line radiation pressure from large set of level 1 lines */
		for( i=1; i <= nLevel1; i++ )
		{
			if( TauLines[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(TauLines[i].EnergyWN/
				  1.e6,4))*(TauLines[i].PopHi/TauLines[i].gHi)/
				  (TauLines[i].PopLo/TauLines[i].gLo)*
				  WidthLine(&TauLines[i])*2.;
				if( RadPres1 > pmx )
				{
					ipSet = 3;
					pmx = RadPres1;
					ip1 = i;
				}
				pressure.RadPres += (float)RadPres1;
			}
		}

		for( i=0; i < nWindLine; i++ )
		{
			if( TauLine2[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(TauLine2[i].EnergyWN/
				  1.e6,4))*(TauLine2[i].PopHi/TauLine2[i].gHi)/
				  (TauLine2[i].PopLo/TauLine2[i].gLo)*
				  WidthLine(&TauLine2[i])*2.;
				pressure.RadPres += (float)RadPres1;
				if( RadPres1 > pmx )
				{
					ipSet = 4;
					pmx = RadPres1;
					ip1 = i;
				}
			}
		}

	}
	else
	{
		/* case where radiation pressure is not turned on */
		ip1 = 0;
		ipSet = 0;
	}

	/* PBETA is ratio of radiation to gas + continuum pressure */
	if( pressure.lgConPres )
	{
		/* this includes integrated incident continuum */
		pressure.pbeta = pressure.RadPres/(pressure.GasPres + pressure.PresInteg);
	}
	else
	{
		/* this does not */
		pressure.pbeta = pressure.RadPres/pressure.GasPres;
	}

	/* following test will never suceed, here to trick lint, ipSet is only used
	 * when needed for debug */
	if( trace.lgTrace && ipSet <0 )
	{
		fprintf(ioQQQ,
			"      TotalPressure, pressure.pbeta = %e, ipSet%li ip1=%li \n", 
			pressure.pbeta,ipSet,ip1 );
	}

	/* stop model from running away on first iteration, when line optical
	 * depths are not defined correctly anyway.
	 * if( iter.le.itermx .and. RadPres.ge.GasPres ) then
	 * >>chng 97 july 23, only cap radiation pressure on first iteration */
	if( IterCnt.iter <= 1 && pressure.RadPres >= pressure.GasPres )
	{
		pressure.RadBetaMax = (float)MAX2(pressure.RadBetaMax,pressure.pbeta);
		/* stop RadPres from exceeding the gas pressure on first iteration */
		pressure.RadPres = (float)MIN2(pressure.RadPres,pressure.GasPres);
		pressure.lgPradCap = TRUE;
	}

	/* radiation pressure only included in total eqn of state when this flag is
	 * true, set with constant pressure command */
	if( pressure.lgRadPresON )
	{
		TotalPressure_v = pressure.GasPres + pressure.RadPres;
	}
	else
	{
		TotalPressure_v = pressure.GasPres;
	}

	/* remember the globally most important line, in the entire model */
	if( pressure.pbeta > pressure.RadBetaMax )
	{
		pressure.RadBetaMax = pressure.pbeta;
		pressure.ipPradMax = ip1;
		if( ipSet == 1 )
		{
			/* special first hydrogenic set */
			strcpy( pressure.chLineRadPres, "First 20" );
		}
		else if( ipSet == 2 )
		{
			/* hydrogenic */
			strcpy( pressure.chLineRadPres, chLineLbl(&HydroLines[ip3][ip1][ip2]) );
		}
		else if( ipSet == 3 )
		{
			/* level 1 lines */
			strcpy( pressure.chLineRadPres, chLineLbl(&TauLines[ip1]) );
		}
		else if( ipSet == 4 )
		{
			/* level 2 lines */
			strcpy( pressure.chLineRadPres, chLineLbl(&TauLine2[ip1]) );
		}
		else
		{
			fprintf(ioQQQ," TotalPressure ipSet insanity.  this is impossible.\n");
			ShowMe();
			puts( "[Stop in TotalPressure]" );
			exit(1);
		}
	}

	/* remember highest pressure anywhere */
	pressure.PresMax = (float)MAX2(pressure.PresMax,TotalPressure_v);


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

double RamPressure(void)
{
	double RamPressure_v;

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

	/* this is ram pressure - xMassDensity is density in gm/cm3,
	 * wind.windv is current wind velocity in cm/s */
	RamPressure_v = densty.xMassDensity * POW2( wind.windv );

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