/* 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 */
/*PresTotCurrent determine the gas and line radiation pressures for current conditions,
 * this sets values of pressure.PresTotlCurr, also calls tfidle */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "iso.h"
#include "wind.h"
#include "magnetic.h"
#include "rfield.h"
#include "phycon.h"
#include "hmi.h"
#include "h2.h"
#include "dense.h"
#include "atomfeii.h"
#include "co.h"
#include "tfidle.h"
#include "dynamics.h"
#include "trace.h"
#include "rt.h"
#include "atmdat.h"
#include "heavy.h"
#include "lines_service.h"
#include "pressure.h"
#include "radius.h"

/* this sets values of pressure.PresTotlCurr, also calls tfidle */
void PresTotCurrent(void)
{
	long int i, 
		ion,
		/* will be used to point to line causing radiation pressure */
	  ip1=-LONG_MAX,
	  ip2=-LONG_MAX,
	  ip3=-LONG_MAX,
	  ip4=-LONG_MAX,
	  ipISO ,
	  ipHi, 
	  ipLo, 
	  ipSet=-1 , /* used in debug print statement to see which line adds most pressure */
	  nelem;

	/* will be used to check on size of opacity, was capped at this value */
	float smallfloat=SMALLFLOAT*10.f;

	double 
		Energy12,
		Energy13,
		Piso_seq[NISO],
		PLevel1=0.,
		PLevel2=0.,
		PHfs = 0.,
		PCO=0.,
		P_H2=0.,
		P_FeII=0.,
	  RadPres1, 
	  TotalPressure_v, 
	  pmx;
	float den_hmole ,
		den_comole ,
		den_ions;

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

	/* PresTotCurrent - evaluate total pressure
	 * this function evaluates the following
	 * PresTotCurrent - gas and radiation pressure, dynes/cm^2
	 * GasPres gas pressure, same units
	 *
	 * HDEN is defined as total hydrogen nuclei phycon; all forms
	 *
	 * dense.xNucleiTotal is total number of atoms, ions, and molecules
	 * nelem,0 is abundance of nelem */

	/* update density - temperature variables which may have changed */
	/* >>chng 02 jun 10, add this here, remove from calling routine */
	tfidle(FALSE);

	/* find total number of nuclei, and total ionization energy of gas */
	den_ions = 0.;
	for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
	{
		for( ion=0; ion<=nelem+1; ++ion )
		{
			den_ions += dense.xIonDense[nelem][ion];
		}
	}

	/* evaluate the total ionization energy on a scale where neutral atoms
	 * correspond to an energy of zero, so the ions have a positive energy */
	phycon.EnergyIonization = 0.;
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
		{
			int kadvec = dynamics.lgMETALS;
			/* this gives the iso sequence for this ion - should it be included in the
			 * advected energy equation*/
			ipISO = nelem-ion;
			if( ipISO >= 0 && ipISO<NISO )
				kadvec = dynamics.lgISO[ipISO] ;
			for( i=1; i<=ion; ++i )
			{
				long int nelec = nelem+1 - i ;
				/* this is the sum of all the energy needed to bring the atom up
				 * to the ion+1 level of ionization - at this point a positive number */
				phycon.EnergyIonization += dense.xIonDense[nelem][ion] * 
					PH1COM.PH1[Heavy.nsShells[nelem][i-1]-1][nelec][nelem][0]/EVRYD* 0.9998787*kadvec;
			}
		}
	}
	/* convert to ergs/cm^3 */
	phycon.EnergyIonization *= EN1RYD;

	/*TODO	2	this is the molecular binding energy (?) not set to anything */
	phycon.EnergyBinding = 0.;

	/* add on molecules of the heavy elements, the atoms and ions of the elements within the
	 * network at beyond NUM_HEAVY_MOLEC and so are not counted.  NUM_COMOLE_CALC is the total
	 * number calculated, including the atoms and ions */
	den_comole = 0.;
	for( i=0; i < NUM_HEAVY_MOLEC; i++ )
	{
		den_comole += co.hevmol[i];
	}

	/* add on hydrogen molecules */
	/* Including ipMHeHp induces instability in fastpn.in 
     * >> Derivative of net cooling negative and so possibly thermally unstable in   1 zones.
	 * but probably no correction for HeH+ has been made to the He abundance, so leaving it out
	 * is correct.
	 * xNucleiTotal is a confusing name once there are molecules! */
	/* if N_H_MOLEC != 8, we will need to change this line... */
	ASSERT(N_H_MOLEC == 8);
	/*>>chng 04 sep 22, make this a loop over all H molecular species */
	den_hmole = 0.;
	for( i=0; i<N_H_MOLEC; ++i )
	{
		if( i!=ipMH && i!=ipMHp && i!=ipMHm )
			den_hmole += hmi.Hmolec[i];
		/*dense.xNucleiTotal += hmi.Hmolec[ipMH2g] + hmi.Hmolec[ipMH2p] + hmi.Hmolec[ipMHm] + hmi.Hmolec[ipMH3p]
		+ hmi.Hmolec[ipMH2s];*/
	}

	/* pden is the total number of particles per unit vol */
	dense.xNucleiTotal = den_hmole + den_comole + den_ions;
	/*fprintf(ioQQQ,"DEBUG xNucleiTotal %.2f %.3e %.3e %.3e %.3e \n",
		fnzone,dense.xNucleiTotal , den_ions , den_hmole , den_comole);*/

	/* this is the particle density that enters into the pressure */
	dense.pden = (float)(dense.eden + dense.xNucleiTotal);

	/* the current gas pressure */
	pressure.PresGasCurr = (float)(dense.pden*phycon.te*BOLTZMANN);
	/*fprintf(ioQQQ,"DEBUG gassss %.2f %.4e %.4e %.4e \n", 
		fnzone, pressure.PresGasCurr , dense.pden , pressure.PresInteg );*/

	/* dense.wmole is molecular weight, AtomicMassUnit per particle */
	dense.wmole = 0.;
	for( i=0; i < LIMELM; i++ )
	{
		dense.wmole += dense.gas_phase[i]*(float)dense.AtomicWeight[i];
	}

	/* dense.wmole is now molecular weight, AtomicMassUnit per particle */
	dense.wmole /= dense.pden;
	ASSERT( dense.wmole > 0. && dense.pden > 0. );

	/* xMassDensity is density in grams per cc */
	/*TODO	2	- should this include mass in grains? */
	dense.xMassDensity = (float)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden);
	
	/*>>chng 04 may 25, introduce folling test on xMassDensity0 
	 * WJH 21 may 04, this is the mass density that corresponds to the hden 
	 * specified in the init file. It is used to calculate the mass flux in the 
	 * dynamics models. It may not necessarily be the same as struc.DenMass[0],
	 * since the pressure solver may have jumped to a different density at the
	 * illuminated face from that specified.*/   
	if ( dense.xMassDensity0 < 0.0 )
		 dense.xMassDensity0 = dense.xMassDensity;
	
	/* >>chng 01 nov 02, move to here from dynamics routines */
	/* >>chng 02 jun 18 (rjrw), fix to use local values */
	/* WJH 21 may 04, now recalculate wind v for the first zone too.
	 * This is necessary when we are forcing the sub or supersonic branch */
	if (wind.windv < 0)
		 wind.windv = DynaFlux(radius.depth)/(dense.xMassDensity);

	/* this is the current ram pressure, at this location */
	pressure.PresRamCurr = dense.xMassDensity*POW2( wind.windv );

	/* RadPres is pressure due to lines, lgRadPresON turns off or on */
	pressure.PresRadCurr = 0.;
	pmx = 0.;

#	if 0
	/* >>chng 04 jun 01, do not treat as special case, include as part of loop below */
	/* this is the total Lya line width, in velocity units, cm/s
	 * >>chng 00 dec 04, this had been inside if( pressure.lgLineRadPresOn but
	 * must always be evaluated since needed for pumping of FeII by Lya */
	/* if statement stops evaluation before initial quantities set */
	if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc > 0. )
	{
		hydro.HLineWidth = (float)(
			RT_LyaWidth(
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot,
			4.72e-2/phycon.sqrte,
			DoppVel.doppler[ipHYDROGEN] ) );
		ASSERT( hydro.HLineWidth >= 0.);
	}
	else
	{
		hydro.HLineWidth = 0.;
	}
#	endif

	/* must always evaluate H La width since used elsewhere */
	if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc > 0. )
	{
		hydro.HLineWidth = (float)RT_LineWidth(&EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s]);
#		if 0
		fprintf(ioQQQ,"DEBUG widLya\t%li\t%.2f\t%.3e",
			iteration,
			fnzone,
			hydro.HLineWidth);
		hydro.HLineWidth = (float)(
			RT_LyaWidth(
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
			EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot,
			4.72e-2/phycon.sqrte,
			DoppVel.doppler[ipHYDROGEN] ) );
		fprintf(ioQQQ,"\t%.3e\n",
			hydro.HLineWidth);
#		endif
	}
	else
		hydro.HLineWidth = 0.;


	/* find max rad pressure even if capped
	 * lgLineRadPresOn is turned off by NO RADIATION PRESSURE command */
	if( pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
	{
		/* >>chng 04 jun 01, only use one loop for all line rad pres, even lya
		 * had been special case done above */
		pressure.PresRadCurr = 0.;
		ipSet = -1;
		for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
		{
			Piso_seq[ipISO] = 0.;
			/* >>chng 01 sep 23, rewrote for iso sequences */
			for( nelem=ipISO; nelem < LIMELM; nelem++ )
			{
				/* does this ion stage exist? */
				if( dense.IonHigh[nelem] >= nelem + 1 - ipISO  )
				{
					/* use 2s - 2p to get pops in 2p */
					double sum_pop_split_term = BIGDOUBLE;
					if( ipISO==ipH_LIKE )
					{
						/* in H-like, all Balmer lines are split into lower 2s and 2p */
						sum_pop_split_term = EmisLines[ipISO][nelem][ipH2s][ipH1s].PopHi + 
							EmisLines[ipISO][nelem][ipH2p][ipH1s].PopHi;
					}
					else if( ipISO==ipHE_LIKE )
					{
						/* in He-like, 2 ^3P is resolved into J=0, 1, 2 */
						sum_pop_split_term = EmisLines[ipISO][nelem][ipHe2p3P0][ipHe2s3S].PopHi + 
							EmisLines[ipISO][nelem][ipHe2p3P1][ipHe2s3S].PopHi +
							EmisLines[ipISO][nelem][ipHe2p3P2][ipHe2s3S].PopHi;
					}
					else
						TotalInsanity();
					/* >>chng 01 jan 08, do not include highest level since maser can occur,
					* and pops are set to small number in this case */
					for( ipLo=0; ipLo < (iso.numLevels[ipISO][nelem] - 2); ipLo++ ) 
					{
						double poplo=BIGDOUBLE;
						int lgPopSet;
						lgPopSet = FALSE;
						if( (ipISO==ipH_LIKE && ipLo >=ipH2s && ipLo <=ipH2p ) )
						{
							if( ipLo ==ipH2s )
							{
								/* for lowest of two, do the rt */
								poplo = sum_pop_split_term;
								lgPopSet = TRUE;
							}
							else if( ipLo==ipH2p )
							{
								poplo = 0;
								continue;
							}
						}
						else if ( ipISO==ipHE_LIKE )
						{
							if( ipLo ==ipHe2p3P0 )
							{
								poplo = sum_pop_split_term;
								lgPopSet = TRUE;
							}
							else if( ipLo > ipHe2p3P0 && ipLo <=ipHe2p3P2 )
							{
								poplo = 0;
								continue;
							}
						}

						for( ipHi=ipLo+1; ipHi < iso.numLevels[ipISO][nelem]-1; ipHi++ )
						{
							/* get lower population - two special cases done above */
							if( !lgPopSet )
								poplo = EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc;

							/* needs to be larger than a small float since this is capped whe
							 * defined*/
							/* >>chng 02 jun 19, add test for small float - spends a lot of time
							 * in this routine, and many he-like transitions are bogus */
							/* >>chng 04 may 13, add Aul > 10 to prevent forbidden lines with
							 * tiny optical depths - taumin doubled - occur, but still opt thin 
							 * sympton was very large radiation pressure from Paschen beta */
							/* >>chng 04 aug 28, add test on ipCont being positive to skip bogus lines*/
							/*NB - this code must be kep parrallel with code in prt_linepres */
							if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont > 0 &&
								EmisLines[ipISO][nelem][ipHi][ipLo].Aul > 10. &&
							    EmisLines[ipISO][nelem][ipHi][ipLo].PopHi > smallfloat &&
								/* >>chng 04 may 21, change PopOpc to poplo as defined above,
								 * to not overcount 23P of He */
								/*EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc > smallfloat &&*/
								poplo > smallfloat &&
								/* >>chng 01 nov 01, add test that have not overrun optical depth scale */
								( (EmisLines[ipISO][nelem][ipHi][ipLo].TauTot - EmisLines[ipISO][nelem][ipHi][ipLo].TauIn) > smallfloat ) &&
								EmisLines[ipISO][nelem][ipHi][ipLo].ipCont>0 &&
								EmisLines[ipISO][nelem][ipHi][ipLo].Pesc > FLT_EPSILON*100. )
							{
								double width = RT_LineWidth(&EmisLines[ipISO][nelem][ipHi][ipLo]);
								RadPres1 = 5.551e-2*
									(powi(EmisLines[ipISO][nelem][ipHi][ipLo].EnergyWN/1.e6,4))*
									(EmisLines[ipISO][nelem][ipHi][ipLo].PopHi/
									iso.stat[ipISO][nelem][ipHi])/
									/* >>chng 04 may 21, change PopOpc to poplo as defined above,
									 * to not overcount 23P of He */
									/*(EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc/ */
									(poplo/
									iso.stat[ipISO][nelem][ipLo])*
									/* RT_LineWidth gets line width in terms of RT effects */
									width;
								if( RadPres1 > pmx )
								{
									ipSet = 2;
									pmx = RadPres1;
									ip4 = ipISO;
									ip3 = nelem;
									ip2 = ipHi;
									ip1 = ipLo;
								}
								Piso_seq[ipISO] += RadPres1;
								pressure.PresRadCurr += (float)RadPres1;
								{
									/* option to print particulars of some line when called */
									/*@-redef@*/
									enum {DEBUG_LOC=FALSE};
									/*@+redef@*/
									if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
									{
										fprintf(ioQQQ,
											"DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
											fnzone,
											RadPres1,
											EmisLines[ipISO][nelem][ipHi][ipLo].PopOpc,
											EmisLines[ipISO][nelem][ipHi][ipLo].PopLo,
											EmisLines[ipISO][nelem][ipHi][ipLo].PopHi,
											EmisLines[ipISO][nelem][ipHi][ipLo].Pesc,
											width);
									}
								}							
							}
						}
					}
				}
			}
			ASSERT( Piso_seq[ipISO] >= 0. );
		}
		{
			/* option to print particulars of some line when called */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
#			if 0
			if( DEBUG_LOC /*&& iteration > 1*/ && nzone > 260 )
			{
				fprintf(ioQQQ,
					"DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
					nzone,
					pmx,
					EmisLines[ipISO][ip3][ip1][ip2].PopOpc,
					EmisLines[ipISO][ip3][7][0].PopLo,
					EmisLines[ipISO][ip3][7][2].PopLo,
					EmisLines[ipISO][ip3][7][3].PopLo,
					EmisLines[ipISO][ip3][7][4].PopLo,
					EmisLines[ipISO][ip3][7][5].PopLo,
					EmisLines[ipISO][ip3][7][6].PopLo);
			}
#			endif
			if( DEBUG_LOC /*&& iteration > 1 && nzone > 150 */)
			{
				fprintf(ioQQQ,
					"DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
					pmx,
					ip4,
					ip3,
					ip2,
					ip1,
					EmisLines[ip4][ip3][ip2][ip1].PopOpc,
					EmisLines[ip4][ip3][ip2][ip1].PopHi,
					1.-EmisLines[ip4][ip3][ip2][ip1].Pesc );
			}
		}

		PLevel1 = 0.;
		/* 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)*
				  RT_LineWidth(&TauLines[i]);
				if( RadPres1 > pmx )
				{
					ipSet = 3;
					pmx = RadPres1;
					ip1 = i;
				}
				PLevel1 += (float)RadPres1;
				pressure.PresRadCurr += (float)RadPres1;
			}
		}
		ASSERT( PLevel1 >= 0. );

		/* level 2 lines */
		PLevel2 = 0.;
		for( i=0; i < nWindLine; i++ )
		{
			if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
			{
				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)*
					RT_LineWidth(&TauLine2[i]);
					pressure.PresRadCurr += (float)RadPres1;
					PLevel2 += (float)RadPres1;
					if( RadPres1 > pmx )
					{
						ipSet = 4;
						pmx = RadPres1;
						ip1 = i;
					}
				}
			}
		}
		ASSERT( PLevel2 >= 0. );

		/* fine structure lines */
		PHfs = 0.;
		for( i=0; i < nHFLines; i++ )
		{
			if( HFLines[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(HFLines[i].EnergyWN/
				  1.e6,4))*(HFLines[i].PopHi/HFLines[i].gHi)/
				  (HFLines[i].PopLo/HFLines[i].gLo)*
				  RT_LineWidth(&HFLines[i]);
				pressure.PresRadCurr += (float)RadPres1;
				PHfs += (float)RadPres1;
				if( RadPres1 > pmx )
				{
					ipSet = 8;
					pmx = RadPres1;
					ip1 = i;
				}
			}
		}
		ASSERT( PHfs >= 0. );

		/* radiation pressure due to H2 lines */
		P_H2 = H2_RadPress();
		/*fprintf(ioQQQ,"DEBUG ph2\t%li\t%.2f\t%.3e\n",
			iteration, fnzone, P_H2 );*/
		/* flag to remember H2 radiation pressure */
		if( P_H2 > pmx )
		{
			pmx = P_H2;
			ipSet = 9;
		}
		pressure.PresRadCurr += (float)P_H2;

		/* radiation pressure due to FeII lines and large atom */
		P_FeII = FeIIRadPress();
		if( P_FeII > pmx )
		{
			pmx = P_FeII;
			ipSet = 7;
		}
		pressure.PresRadCurr += (float)P_FeII;

		PCO = 0.;
		/* co carbon monoxide lines */
		for( i=0; i < nCORotate; i++ )
		{
			if( C12O16Rotate[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(C12O16Rotate[i].EnergyWN/
				  1.e6,4))*(C12O16Rotate[i].PopHi/C12O16Rotate[i].gHi)/
				  (C12O16Rotate[i].PopLo/C12O16Rotate[i].gLo)*
				  RT_LineWidth(&C12O16Rotate[i]);
				pressure.PresRadCurr += (float)RadPres1;
				PCO += (float)RadPres1;
				if( RadPres1 > pmx )
				{
					ipSet = 5;
					pmx = RadPres1;
					ip1 = i;
				}
			}
			if( C13O16Rotate[i].PopHi > 1e-30 )
			{
				RadPres1 = 5.551e-2*(powi(C13O16Rotate[i].EnergyWN/
				  1.e6,4))*(C13O16Rotate[i].PopHi/C13O16Rotate[i].gHi)/
				  (C13O16Rotate[i].PopLo/C13O16Rotate[i].gLo)*
				  RT_LineWidth(&C13O16Rotate[i]);
				pressure.PresRadCurr += (float)RadPres1;
				PCO += (float)RadPres1;
				if( RadPres1 > pmx )
				{
					ipSet = 6;
					pmx = RadPres1;
					ip1 = i;
				}
			}
		}
		ASSERT( PCO >= 0. );
	}
	else
	{
		/* case where radiation pressure is not turned on */
		ip1 = 0;
		ipSet = 0;
	}

	/* PBETA is ratio of radiation to total (gas + continuum + etc) pressure */
	/* >>chng 04 may 15, define pressure beta as radiation to total pressure rather than just gas */
	pressure.pbeta = (float)(pressure.PresRadCurr/SDIV(pressure.PresTotlCurr));
#	if 0
	if( pressure.lgContRadPresOn )
	{
		/* this includes integrated incident continuum */
		pressure.pbeta = (float)(pressure.PresRadCurr/(pressure.PresGasCurr + pressure.PresInteg));
	}
	else
	{
		/* this does not */
		pressure.pbeta = (float)(pressure.PresRadCurr/pressure.PresGasCurr);
	}
#	endif

	/*lint -e644 Piso_seq not init */
	if( trace.lgTrace && (pressure.PresTotlCurr > SMALLFLOAT) )
	{
		fprintf(ioQQQ,
			"     PresTotCurrent %.2f Ptot:%.3e Pgas:%.3e Prads:%.2e , constituents follow H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e\n",
		fnzone,
		pressure.PresTotlCurr, 
		pressure.PresGasCurr/pressure.PresTotlCurr,
		pressure.PresRadCurr/pressure.PresTotlCurr,
		Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
		Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
		PLevel1/pressure.PresTotlCurr,
		PLevel2/pressure.PresTotlCurr,
		PCO/pressure.PresTotlCurr,
		PHfs/pressure.PresTotlCurr,
		P_H2/pressure.PresTotlCurr);
		/*lint +e644 Piso_seq not initialized */
	}

	/* this would be a major logical error .. .. . */
	if( pressure.PresRadCurr < 0. )
	{
		fprintf(ioQQQ,
			"PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
		Piso_seq[ipH_LIKE],
		Piso_seq[ipHE_LIKE],
		PLevel1,
		PLevel2,
		PCO);
		ShowMe();
		puts( "[Stop in PresTotCurrent]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	/* this is the total internal energy of the gas, the amount of energy needed
	 * to produce the current level populations, relative to ground,
	 * only used for energy terms in advection */
	phycon.EnergyExcitation = 0.;
	ipLo = 0;
	for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
	{
		for( nelem=ipISO; nelem < LIMELM; nelem++ )
		{
			if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
			{
				for( ipHi=1; ipHi < iso.numLevels[ipISO][nelem]; ipHi++ ) 
				{
					phycon.EnergyExcitation += 
						EmisLines[ipISO][nelem][ipHi][ipLo].PopHi * 
						EmisLines[ipISO][nelem][ipHi][ipLo].EnergyErg*
						/* last term is option to turn off energy term, no advection hlike, he-like */
						dense.xIonDense[nelem][nelem+1-ipISO]*dynamics.lgISO[ipISO];
				}
			}
		}
	}

	if( dynamics.lgISO[ipH_LIKE] )
		/* internal energy of H2 */
		phycon.EnergyExcitation += H2_InterEnergy();

	/* this is option to turn off energy effects of advection, no advection metals */
	if( dynamics.lgMETALS )
	{
		for( i=1; i <= nLevel1; i++ )
		{
			phycon.EnergyExcitation += 
				TauLines[i].PopHi* TauLines[i].EnergyErg;
		}
		for( i=0; i < nWindLine; i++ )
		{
			if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
			{
				phycon.EnergyExcitation += 
					TauLine2[i].PopHi* TauLine2[i].EnergyErg;
			}
		}
		for( i=0; i < nHFLines; i++ )
		{
			phycon.EnergyExcitation += 
				HFLines[i].PopHi* HFLines[i].EnergyErg;
		}
		Energy12 = 0.;
		Energy13 = 0.;
		for( i=0; i < nCORotate; i++ )
		{
			Energy12 += C12O16Rotate[i].EnergyErg;
			phycon.EnergyExcitation += 
				C12O16Rotate[i].PopHi* Energy12;

			Energy13 += C13O16Rotate[i].EnergyErg;
			phycon.EnergyExcitation += 
				C13O16Rotate[i].PopHi* Energy13;
		}

		/* internal energy of large FeII atom */
		phycon.EnergyExcitation += FeII_InterEnergy();
	}

	/* ==================================================================
	 * end internal energy of atoms and molecules */

	/* evaluate some parameters to do with magnetic field */
	Magnetic_evaluate();

	/* Conserved quantity in steady-state energy equation for dynamics:
	 * thermal energy, since recombination is treated as cooling
	 * (i.e. is loss of electron k.e. to emitted photon), so don't
	 * include
	 * ...phycon.EnergyExcitation + phycon.EnergyIonization + phycon.EnergyBinding
	 * */

	phycon.EnthalpyDensity =  
		0.5*POW2(wind.windv)*dense.xMassDensity +	/* KE */
		5./2.*pressure.PresGasCurr +				/* thermal plus PdV work */
		magnetic.EnthalpyDensity;						/* magnetic terms */

	/* 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 jul 23, only cap radiation pressure on first iteration */
	if( iteration <= 1 && pressure.PresRadCurr >= pressure.PresGasCurr )
	{
		pressure.RadBetaMax = (float)MAX2(pressure.RadBetaMax,pressure.pbeta);
		/* stop RadPres from exceeding the gas pressure on first iteration */
		pressure.PresRadCurr = (float)MIN2(pressure.PresRadCurr,pressure.PresGasCurr);
		pressure.lgPradCap = TRUE;
	}

	/* remember the globally most important line, in the entire model */
	if( pressure.pbeta > pressure.RadBetaMax )
	{
		pressure.RadBetaMax = pressure.pbeta;
		pressure.ipPradMax_line = ip1;
		pressure.ipPradMax_nzone = nzone;
		if( ipSet == 2 )
		{
			/* hydrogenic */
			strcpy( pressure.chLineRadPres, "ISO   ");
			ASSERT( ip4 < NISO && ip3<LIMELM );
			ASSERT( ip1>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
			strcat( pressure.chLineRadPres, chLineLbl(&EmisLines[ip4][ip3][ip2][ip1]) );
			{
				/* option to print particulars of some line when called */
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				/*lint -e644 Piso_seq not initialized */
				/* trace verious constituents in radiation pressure */
				if( DEBUG_LOC  )
				{
					fprintf(ioQQQ,"DEBUG iso prad\t%li\t%li\tISO,nelem=\t%li\t%li\tnu, no=\t%li\t%li\t%.4e\t%.4e\t%e\t%e\t%e\n",
					iteration, nzone, 
					ip4,ip3,ip2,ip1,
					EmisLines[ip4][ip3][ip2][ip1].TauIn,
					EmisLines[ip4][ip3][ip2][ip1].TauTot,
					EmisLines[ip4][ip3][ip2][ip1].Pesc,
					EmisLines[ip4][ip3][ip2][ip1].Pelec_esc,
					EmisLines[ip4][ip3][ip2][ip1].Pdest);
					if( ip2==5 && ip1==4 )
					{
						double width;
						fprintf(ioQQQ,"hit it\n");
						width = RT_LineWidth(&EmisLines[ip4][ip3][ip2][ip1]);
						fprintf(ioQQQ,"DEBUG width %e\n", width);
					}
				}
			}
		}
		else if( ipSet == 3 )
		{
			/* level 1 lines */
			ASSERT( ip1>=0  );
			strcpy( pressure.chLineRadPres, "Level1 ");
			strcat( pressure.chLineRadPres, chLineLbl(&TauLines[ip1]) );
		}
		else if( ipSet == 4 )
		{
			/* level 2 lines */
			ASSERT( ip1>=0  );
			strcpy( pressure.chLineRadPres, "Level2 ");
			strcat( pressure.chLineRadPres, chLineLbl(&TauLine2[ip1]) );
		}
		else if( ipSet == 5 )
		{
			/* c12o16 carbon monoxide lines */
			ASSERT( ip1>=0  );
			strcpy( pressure.chLineRadPres, "12CO ");
			strcat( pressure.chLineRadPres, chLineLbl(&C12O16Rotate[ip1]) );
		}
		else if( ipSet == 6 )
		{
			/* c13o16 carbon monoxide lines */
			ASSERT( ip1>=0  );
			strcpy( pressure.chLineRadPres, "13CO ");
			strcat( pressure.chLineRadPres, chLineLbl(&C13O16Rotate[ip1]) );
		}
		else if( ipSet == 7 )
		{
			/* helium lines */
			strcpy( pressure.chLineRadPres, "Fe II ");
		}
		else if( ipSet == 8 )
		{
			/* helium lines */
			strcpy( pressure.chLineRadPres, "hyperf ");
			ASSERT( ip1>=0  );
			strcat( pressure.chLineRadPres, chLineLbl(&HFLines[ip1]) );
		}
		else if( ipSet == 9 )
		{
			/* helium lines */
			strcpy( pressure.chLineRadPres, "large H2 ");
		}
		else
		{
			fprintf(ioQQQ," PresTotCurrent ipSet set to %li, this is impossible.\n", ipSet);
			ShowMe();
			puts( "[Stop in PresTotCurrent]" );
			cdEXIT(EXIT_FAILURE);
		}
	}

	/* radiation pressure only included in total eqn of state when this flag is
	 * true, set with constant pressure command */
	/* >>chng 01 nov 02, add ram pressure here */

	/* >>chng 02 jun 27 - add in magnetic pressure */
	/* start to bring total pressure together */
	TotalPressure_v = pressure.PresGasCurr + pressure.PresRamCurr + magnetic.pressure;
	
	/* option to add in internal line radiation pressure */
	if( pressure.lgRadPresON )
	{
		TotalPressure_v += pressure.PresRadCurr;
	}

	{
		/*@-redef@*/
		enum{DEBUG_LOC=FALSE};
		/*@+redef@*/
		if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT /*&& iteration > 1*/ )
		{
			fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n", 
				nzone,
				pressure.PresTotlCorrect,
			    pressure.PresTotlCurr, 
				TotalPressure_v ,
				pressure.PresGasCurr/pressure.PresTotlCurr,
				pressure.PresRadCurr/pressure.PresTotlCurr ,
				pressure.PresRamCurr/pressure.PresTotlCurr
				);
		}
	}

	ASSERT( TotalPressure_v > 0. );

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

	/* this is what we came for - set the current pressure */
	pressure.PresTotlCurr = TotalPressure_v;

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

	return;
}

