cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pressure_total.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*PresTotCurrent determine the gas and line radiation pressures for current conditions,
4  * this sets values of pressure.PresTotlCurr, also calls tfidle */
5 #include "cddefines.h"
6 #include "taulines.h"
7 #include "opacity.h"
8 #include "hydrogenic.h"
9 #include "conv.h"
10 #include "iso.h"
11 #include "wind.h"
12 #include "magnetic.h"
13 #include "phycon.h"
14 #include "thermal.h"
15 #include "h2.h"
16 #include "dense.h"
17 #include "dynamics.h"
18 #include "trace.h"
19 #include "rt.h"
20 #include "pressure.h"
21 #include "radius.h"
22 #include "rfield.h"
23 #include "doppvel.h"
24 #include "rt_escprob.h"
25 
26 /* this sets values of pressure.PresTotlCurr, also calls tfidle */
27 void PresTotCurrent(void)
28 {
29  static long int
30  /* used in debug print statement to see which line adds most pressure */
31  ipLineTypePradMax=-1 ,
32  /* points to line causing radiation pressure */
33  ipLinePradMax=-LONG_MAX,
34  ip2=-LONG_MAX,
35  ip3=-LONG_MAX,
36  ip4=-LONG_MAX;
37 
38  /* the line radiation pressure variables that must be preserved since
39  * a particular line radiation pressure may not be evaluated if it is
40  * not important */
41  static double Piso_seq[NISO]={0.,0.},
42  PLevel2=0.,
43  PHfs=0.,
44  PCO=0.,
45  P_H2=0.,
46  P_dBase=0.;
47 
48  double
49  RadPres1,
50  TotalPressure_v,
51  pmx;
52 
53  DEBUG_ENTRY( "PresTotCurrent()" );
54 
55  if( !conv.nTotalIoniz )
56  {
57  /* resetting ipLinePradMax, necessary for
58  * multiple cloudy calls in single coreload. */
59  ipLinePradMax=-LONG_MAX;
60  //pressure.PresTotlCurr = 0.;
61  }
62 
63  /* PresTotCurrent - evaluate total pressure, dyne cm^-2
64  * and radiative acceleration */
65 
66  /* several loops over the emission lines for radiation pressure and
67  * radiative acceleration are expensive due to the number of lines and
68  * the memory they occupy. Many equations of state do not include
69  * radiation pressure or radiative acceleration. Only update these
70  * when included in EOS - when not included only evaluate once per zone.
71  * do it once per zone since we will still report these quantities.
72  * this flag says whether we must update all terms */
73  bool lgMustEvaluate = false;
74 
75  /* this says we already have an evaluation in this zone so do not
76  * evaluate terms known to be trivial, even when reevaluating major terms */
77  bool lgMustEvaluateTrivial = false;
78  /* if an individual agent is larger than this fraction of the total current
79  * radiation pressure then it is not trivial
80  * conv.PressureErrorAllowed is relative error allowed in pressure */
81  double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. *
83 
84  /* reevaluate during search phase when conditions are changing dramatically */
85  if( conv.lgSearch )
86  {
87  lgMustEvaluate = true;
88  lgMustEvaluateTrivial = true;
89  }
90 
91  /* reevaluate if zone or iteration has changed */
92  static long int nzoneEvaluated=0, iterationEvaluated=0;
93  if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated )
94  {
95  lgMustEvaluate = true;
96  lgMustEvaluateTrivial = true;
97  /* this flag says which source of radiation pressure was strongest */
98  ipLineTypePradMax = -1;
99  }
100 
101  /* constant pressure or dynamical sim - we must reevaluate terms
102  * but do not need to reevaluate trivial contributors */
103  if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) ||
104  (strcmp(dense.chDenseLaw,"DYNA") == 0 ) ||
105  (strcmp(dense.chDenseLaw,"CPRE") == 0 ) )
106  lgMustEvaluate = true;
107 
108  if( lgMustEvaluate )
109  {
110  /* we are reevaluating quantities in this zone and iteration,
111  * remember that we did it */
112  nzoneEvaluated = nzone;
113  iterationEvaluated = iteration;
114  }
115 
116  /* update density - temperature variables which may have changed */
117  /* must call TempChange since ionization has changed, there are some
118  * terms that affect collision rates (H0 term in electron collision) */
119  TempChange(phycon.te , false);
120 
122 
123  /* evaluate the total ionization energy on a scale where neutral atoms
124  * correspond to an energy of zero, so the ions have a positive energy */
126 #if 0
127  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
128  {
129  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
130  {
131  /* lgMETALS is true, set false with "no advection metals" command */
132  int kadvec = dynamics.lgMETALS;
133  /* this gives the iso sequence for this ion - should it be included in the
134  * advected energy equation? lgISO is true, set false with
135  * "no advection H-like" and "he-like" - for testing*/
136  ipISO = nelem-ion;
137  fixit(); /* should this be kadvec = kadvec && dynamics.lgISO[ipISO]; ? */
138  if( ipISO >= 0 && ipISO<NISO )
139  kadvec = dynamics.lgISO[ipISO];
140  for( long i=1; i<=ion; ++i )
141  {
142  long int nelec = nelem+1 - i;
143  /* this is the sum of all the energy needed to bring the atom up
144  * to the ion+1 level of ionization - at this point a positive number */
145  phycon.EnergyIonization += dense.xIonDense[nelem][ion] *
146  t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD*kadvec;
147  }
148  }
149  }
150  /* convert to ergs/cm^3 */
151  phycon.EnergyIonization *= EN1RYD;
152 #endif
153 
157  phycon.EnergyBinding = 0.;
158 
159  /* find total number of atoms and ions */
160  SumDensities();
161 
162  /* the current gas pressure */
163  pressure.PresGasCurr = dense.pden*phycon.te*BOLTZMANN;
164  /*fprintf(ioQQQ,"DEBUG gassss %.2f %.4e %.4e %.4e \n",
165  fnzone, pressure.PresGasCurr , dense.pden , pressure.PresInteg );*/
166 
167  /* >>chng 01 nov 02, move to here from dynamics routines */
168  /* >>chng 02 jun 18 (rjrw), fix to use local values */
169  /* WJH 21 may 04, now recalculate wind v for the first zone too.
170  * This is necessary when we are forcing the sub or supersonic branch */
171  if(!(wind.lgBallistic() || wind.lgStatic()))
172  {
174  }
175 
176  /* this is the current ram pressure, at this location */
178 
179  /* this is the current turbulent pressure, not now added to equation of state
180  * >>chng 06 mar 29, add Heiles_Troland_F / 6. term*/
183 
186  /* radiative acceleration, lines and continuum */
187  if( lgMustEvaluate )
188  {
189  /* continuous radiative acceleration */
190  double rforce = 0.;
191  double relec = 0.;
192  for( long i=0; i < rfield.nflux; i++ )
193  {
194  rforce += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
195  rfield.anu(i)*(opac.opacity_abs[i] + opac.opacity_sct[i]);
196 
197  /* electron scattering acceleration - used only to derive force multiplier */
198  relec +=
199  (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+
200  rfield.ConInterOut[i]) *
202  }
203  /* total continuum radiative acceleration */
204  wind.AccelCont = (realnum)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
205 
206  wind.AccelElectron = (realnum)(relec*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
207 
208  /* line acceleration; xMassDensity is gm per cc */
210 
211  /* total acceleration cm s^-2 */
213 
214  /* find wind.fmul - the force multiplier; wind.AccelElectron can be zero for very low
215  * densities - fmul is of order unity - wind.AccelLine and wind.AccelCont
216  * are both floats to will underflow long before wind.AccelElectron will - fmul is only used
217  * in output, not in any physics */
220  else
221  wind.fmul = 0.;
222 
223  /* inward acceleration of gravity cm s^-2 */
224  if (wind.comass == 0.)
225  {
226  wind.AccelGravity = 0.;
227  }
228  else
229  {
230  double reff = radius.Radius-radius.drad/2.;
232  /* wind.comass - mass of star in solar masses */
233  GRAV_CONST*wind.comass*SOLAR_MASS/POW2(reff)*
234  /* wind.DiskRadius normally zero, set with disk option on wind command */
235  (1.-wind.DiskRadius/reff) );
236  }
237 
238 # if 0
239  if( fudge(-1) )
240  fprintf(ioQQQ,"DEBUG pressure_total updates AccelTotalOutward to %.2e grav %.2e\n",
242 # endif
243  }
244 
245  /* must always evaluate H La width since used elsewhere */
246  if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().PopOpc() > 0. )
247  {
249  }
250  else
251  hydro.HLineWidth = 0.;
252 
253 
254  /* find max rad pressure even if capped
255  * lgLineRadPresOn is turned off by NO RADIATION PRESSURE command */
256  if( trace.lgTrace )
257  {
258  fprintf(ioQQQ,
259  " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
260  "lgMustEvaluateTrivial:%c "
261  "pressure.lgLineRadPresOn:%c "
262  "rfield.lgDoLineTrans:%c \n",
263  nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial),
265  }
266 
267  if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
268  {
269  /* RadPres is pressure due to lines, lgPres_radiation_ON turns off or on */
271  /* used to remember largest radiation pressure source */
272  pmx = 0.;
273  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
274  {
275  if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
276  {
277  Piso_seq[ipISO] = 0.;
278  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
279  {
280  /* does this ion stage exist? */
281  if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
282  {
283  /* do not include highest levels since maser can occur due to topoff,
284  * and pops are set to small number in this case */
285  for( long ipHi=1; ipHi <iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
286  {
287  for( long ipLo=0; ipLo < ipHi; ipLo++ )
288  {
289  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
290  continue;
291 
292  ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > iso_ctrl.SmallA );
293 
294  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() > SMALLFLOAT &&
295  /* test that have not overrun optical depth scale */
296  ( (iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
297  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()) > SMALLFLOAT ) &&
298  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() > FLT_EPSILON*100. )
299  {
300  RadPres1 = PressureRadiationLine( iso_sp[ipISO][nelem].trans(ipHi,ipLo), GetDopplerWidth(dense.AtomicWeight[nelem]) );
301 
302  if( RadPres1 > pmx )
303  {
304  ipLineTypePradMax = 2;
305  pmx = RadPres1;
306  ip4 = ipISO;
307  ip3 = nelem;
308  ip2 = ipHi;
309  ipLinePradMax = ipLo;
310  }
311  Piso_seq[ipISO] += RadPres1;
312  {
313  /* option to print particulars of some line when called */
314  enum {DEBUG_LOC=false};
315  if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
316  {
317  fprintf(ioQQQ,
318  "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
319  fnzone,
320  RadPres1,
321  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc(),
322  iso_sp[ipISO][nelem].st[ipLo].Pop(),
323  iso_sp[ipISO][nelem].st[ipHi].Pop(),
324  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc());
325  }
326  }
327  }
328  }
329  }
330  }
331  }
332  ASSERT( Piso_seq[ipISO] >= 0. );
333  }
334  pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
335  }
336  {
337  /* option to print particulars of some line when called */
338  enum {DEBUG_LOC=false};
339 # if 0
340  if( DEBUG_LOC /*&& iteration > 1*/ && nzone > 260 )
341  {
342  fprintf(ioQQQ,
343  "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
344  nzone,
345  pmx,
346  iso_sp[ipISO][ip3].trans(ipLinePradMax,ip2).Emis().PopOpc(),
347  iso_sp[ipISO][ip3].st[0].Pop(),
348  iso_sp[ipISO][ip3].st[2].Pop(),
349  iso_sp[ipISO][ip3].st[3].Pop(),
350  iso_sp[ipISO][ip3].st[4].Pop(),
351  iso_sp[ipISO][ip3].st[5].Pop(),
352  iso_sp[ipISO][ip3].st[6].Pop());
353  }
354 # endif
355  if( DEBUG_LOC /*&& iteration > 1 && nzone > 150 */)
356  {
357  fprintf(ioQQQ,
358  "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
359  pmx,
360  ip4,
361  ip3,
362  ip2,
363  ipLinePradMax,
364  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().PopOpc(),
365  iso_sp[ip4][ip3].st[ip2].Pop(),
366  1.-iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc() );
367  }
368  }
369 
370  if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
371  {
372  /* level 2 lines */
373  PLevel2 = 0.;
374  for( long i=0; i < nWindLine; i++ )
375  {
376  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
377  {
378  if( (*TauLine2[i].Hi()).Pop() > 1e-30 )
379  {
380  RadPres1 = PressureRadiationLine( TauLine2[i], GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
381 
382  PLevel2 += RadPres1;
383  if( RadPres1 > pmx )
384  {
385  ipLineTypePradMax = 4;
386  pmx = RadPres1;
387  ipLinePradMax = i;
388  }
389  }
390  }
391  }
392  ASSERT( PLevel2 >= 0. );
393  }
395 
396  /* fine structure lines */
397  if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
398  {
399  PHfs = 0.;
400  for( size_t i=0; i < HFLines.size(); i++ )
401  {
402  if( (*HFLines[i].Hi()).Pop() > 1e-30 )
403  {
404  RadPres1 = PressureRadiationLine( HFLines[i], GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
405 
406  PHfs += RadPres1;
407  if( RadPres1 > pmx )
408  {
409  ipLineTypePradMax = 8;
410  pmx = RadPres1;
411  ipLinePradMax = i;
412  }
413  }
414  }
415  ASSERT( PHfs >= 0. );
416  }
418 
419  /* radiation pressure due to H2 lines */
420  if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
421  {
422  P_H2 = 0.;
423  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
424  {
425  P_H2 += (*diatom)->H2_RadPress();
426  /* flag to remember H2 radiation pressure */
427  if( P_H2 > pmx )
428  {
429  pmx = P_H2;
430  ipLineTypePradMax = 9;
431  }
432  ASSERT( P_H2 >= 0. );
433  }
434  }
436 
437  /* do lines from third-party databases */
438  if( lgMustEvaluateTrivial || P_dBase > TrivialLineRadiationPressure )
439  {
440  P_dBase = 0.;
441  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
442  {
443  if( dBaseSpecies[ipSpecies].lgActive )
444  {
445  realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
446  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
447  tr != dBaseTrans[ipSpecies].end(); ++tr)
448  {
449  int ipHi = (*tr).ipHi();
450  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
451  continue;
452  int ipLo = (*tr).ipLo();
453  if( (*tr).ipCont() > 0 && (*(*tr).Hi()).Pop() > 1e-30 )
454  {
455  RadPres1 = PressureRadiationLine( *tr, DopplerWidth );
456 
457  if( RadPres1 > pmx )
458  {
459  ipLineTypePradMax = 10;
460  pmx = RadPres1;
461  ip3 = ipSpecies;
462  ip2 = ipHi;
463  ipLinePradMax = ipLo;
464  }
465  P_dBase += RadPres1;
466  }
467  }
468  }
469  }
470  ASSERT( P_dBase >= 0. );
471  }
473 
474  }
476  {
477  /* case where radiation pressure is not turned on */
478  ipLinePradMax = -1;
479  ipLineTypePradMax = 0;
480  }
481 
482  fixit("all other line stacks need to be included here");
483  // can we just sweep over line stack? Is that ready yet?
484 
485  /* the ratio of radiation to total (gas + continuum + etc) pressure */
487 
488  /* this would be a major logical error */
490  {
491  fprintf(ioQQQ,
492  "PresTotCurrent: negative pressure, constituents follow %e %e %e %e \n",
493  Piso_seq[ipH_LIKE],
494  Piso_seq[ipHE_LIKE],
495  PLevel2,
496  PCO);
497  ShowMe();
499  }
500 
501  /* following test will never succeed, here to trick lint, ipLineTypePradMax is only used
502  * when needed for debug */
503  if( trace.lgTrace && ipLineTypePradMax <0 )
504  {
505  fprintf(ioQQQ,
506  " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
507  pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
508  }
509 
510  /* this is the total internal energy of the gas, the amount of energy needed
511  * to produce the current level populations, relative to ground,
512  * only used for energy terms in advection */
514 #if 0
515  fixit(); /* the quantities phycon.EnergyExcitation, phycon.EnergyIonization, phycon.EnergyBinding
516  * are not used anywhere, except in print statements... */
517  broken(); /* the code below contains serious bugs! It is supposed to implement loops
518  * over quantum states in order to evaluate the potential energy stored in
519  * excited states of all atoms, ions, and molecules. The code below however
520  * often implements loops over all combinations of upper and lower levels!
521  * This code needs to be rewritten once quantumstates are fully implemented. */
522  ipLo = 0;
523  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
524  {
525  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
526  {
527  if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
528  {
529  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
530  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
531  {
533  iso_sp[ipISO][nelem].st[ipHi].Pop() *
534  iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() *
535  /* last term is option to turn off energy term, no advection hlike, he-like */
536  dynamics.lgISO[ipISO];
537  }
538  }
539  }
540  }
541 
542  if( dynamics.lgISO[ipH_LIKE] )
543  /* internal energy of H2 */
544  phycon.EnergyExcitation += H2_InterEnergy();
545 
546  /* this is option to turn off energy effects of advection, no advection metals */
547  if( dynamics.lgMETALS )
548  {
549  for( long i=0; i < nWindLine; i++ )
550  {
551  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
552  {
554  (*TauLine2[i].Hi()).Pop()* TauLine2[i].EnergyErg;
555  }
556  }
557  for( long i=0; i < nHFLines; i++ )
558  {
560  (*HFLines[i].Hi()).Pop()* HFLines[i].EnergyErg;
561  }
562 
563  /* internal energy of large FeII atom */
564  }
565 #endif
566 
567  /* ==================================================================
568  * end internal energy of atoms and molecules */
569 
570  /* evaluate some parameters to do with magnetic field */
572 
573  /*lint -e644 Piso_seq not init */
575  {
576  fprintf(ioQQQ,
577  " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
578  "gas parts; H:%.2e He:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
579  fnzone,
585  Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
586  Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
587  PLevel2/pressure.PresTotlCurr,
589  PHfs/pressure.PresTotlCurr,
590  P_H2/pressure.PresTotlCurr,
592  /*lint +e644 Piso_seq not initialized */
593  }
594 
595  /* Conserved quantity in steady-state energy equation for dynamics:
596  * thermal energy, since recombination is treated as cooling
597  * (i.e. is loss of electron k.e. to emitted photon), so don't
598  * include
599  * ...phycon.EnergyExcitation + phycon.EnergyIonization + phycon.EnergyBinding
600  * */
601 
602  /* constant density case is also hypersonic case */
604  {
605  /* this branch is time dependent single-zone */
606  /* \todo 1 this has 3/2 on the PresGasCurr while true dynamics case below
607  * has 5/2 - so this is not really the enthalpy density - better
608  * would be to always use this term and add the extra PresGasCurr
609  * when the enthalpy is actually needed */
611  0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
612  3./2.*pressure.PresGasCurr + /* thermal plus PdV work */
613  magnetic.EnthalpyDensity; /* magnetic terms */
614  //pressure.RhoGravity * (radius.Radius/(radius.Radius+radius.drad)); /* gravity */
615  }
616  else
617  {
618  /* this branch either advective or constant pressure */
619  /*fprintf(ioQQQ,"DEBUG enthalpy HIT2\n");*/
620  /* this is usual dynamics case, or time-varying case where pressure
621  * is kept constant and PdV work is added to the cell */
623  0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
624  5./2.*pressure.PresGasCurr + /* thermal plus PdV work */
625  magnetic.EnthalpyDensity; /* magnetic terms */
626  //pressure.RhoGravity * (radius.Radius/(radius.Radius+radius.drad)); /* gravity */
627  }
628 
629  /* stop model from running away on first iteration, when line optical
630  * depths are not defined correctly anyway.
631  * if( iter.le.itermx .and. RadPres.ge.GasPres ) then
632  * >>chng 97 jul 23, only cap radiation pressure on first iteration */
633  if( iteration <= 1 && pressure.pres_radiation_lines_curr >= pressure.PresGasCurr )
634  {
635  /* stop RadPres from exceeding the gas pressure on first iteration */
638  pressure.lgPradCap = true;
639  }
640 
641  /* remember the globally most important line, in the entire model
642  * test on nzone so we only do test after solution is stable */
644  {
646  pressure.ipPradMax_line = ipLinePradMax;
648  if( ipLineTypePradMax == 2 )
649  {
650  /* hydrogenic */
651  pressure.chLineRadPres = "ISO ";
652  ASSERT( ip4 < NISO && ip3<LIMELM );
653  ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
654  pressure.chLineRadPres += chLineLbl(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax));
655  {
656  /* option to print particulars of some line when called */
657  enum {DEBUG_LOC=false};
658  /*lint -e644 Piso_seq not initialized */
659  /* trace serious constituents in radiation pressure */
660  if( DEBUG_LOC )
661  {
662  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",
663  iteration, nzone,
664  ip4,ip3,ip2,ipLinePradMax,
665  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauIn(),
666  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauTot(),
667  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc(),
668  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pelec_esc(),
669  iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pdest());
670  if( ip2==5 && ipLinePradMax==4 )
671  {
672  double width;
673  fprintf(ioQQQ,"hit it\n");
674  width = RT_LineWidth(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax),GetDopplerWidth(dense.AtomicWeight[ip3]));
675  fprintf(ioQQQ,"DEBUG width %e\n", width);
676  }
677  }
678  }
679  }
680  else if( ipLineTypePradMax == 4 )
681  {
682  /* level 2 lines */
683  ASSERT( ipLinePradMax>=0 );
684  pressure.chLineRadPres = "Level2 ";
685  pressure.chLineRadPres += chLineLbl(TauLine2[ipLinePradMax]);
686  }
687  else if( ipLineTypePradMax == 5 )
688  {
689  cdEXIT( EXIT_FAILURE );
690  }
691  else if( ipLineTypePradMax == 6 )
692  {
693  cdEXIT( EXIT_FAILURE );
694  }
695  else if( ipLineTypePradMax == 7 )
696  {
697  /* FeII lines */
698  pressure.chLineRadPres = "Fe II ";
699  }
700  else if( ipLineTypePradMax == 8 )
701  {
702  /* hyperfine struct lines */
703  pressure.chLineRadPres = "hyperf ";
704  ASSERT( ipLinePradMax>=0 );
705  pressure.chLineRadPres += chLineLbl(HFLines[ipLinePradMax]);
706  }
707  else if( ipLineTypePradMax == 9 )
708  {
709  /* large H2 lines */
710  pressure.chLineRadPres = "large H2 ";
711  }
712  else if( ipLineTypePradMax == 10 )
713  {
714  /* database lines */
715  pressure.chLineRadPres = "dBaseLin " ;
716  pressure.chLineRadPres += dBaseSpecies[ip3].chLabel ;
717  }
718  else
719  {
720  fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
721  ShowMe();
723  }
724  }
725 
726  if( trace.lgTrace && pressure.pbeta > 0.5 )
727  {
728  fprintf(ioQQQ,
729  " PresTotCurrent Pbeta:%.2f due to %s\n",
730  pressure.pbeta ,
731  pressure.chLineRadPres.c_str()
732  );
733  }
734 
735  /* >>chng 02 jun 27 - add in magnetic pressure */
736  /* start to bring total pressure together */
737  TotalPressure_v = pressure.PresGasCurr;
738 
739  /* these flags are set false by default since constant density is default,
740  * set true for constant pressure or dynamics */
741  TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON;
742 
743  /* magnetic pressure, evaluated in magnetic.c - this can be negative for an ordered field!
744  * option is on by default, turned off with constant density, or constant gas pressure, cases */
747  TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON;
748 
749  /* turbulent pressure
750  * >>chng 06 mar 24, include this by default */
751  TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure;
752 
753  /* radiation pressure only included in total eqn of state when this flag is
754  * true, set with constant pressure command */
755  /* option to add in internal line radiation pressure */
757 
758  {
759  enum{DEBUG_LOC=false};
760  if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT /*&& iteration > 1*/ )
761  {
762  fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
763  nzone,
766  TotalPressure_v ,
770  );
771  }
772  }
773 
774  if( TotalPressure_v< 0. )
775  {
776  ASSERT( magnetic.pressure < 0. );
777 
778  /* negative pressure due to ordered field overwhelms total pressure - punt */
779  fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total outward pressure - please reconsider the geometry & field.\n");
781  }
782 
783  ASSERT( TotalPressure_v > 0. );
784 
785  /* remember highest pressure anywhere */
786  pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v);
787 
788  /* this is what we came for - set the current pressure */
789  pressure.PresTotlCurr = TotalPressure_v;
790 
791  return;
792 }
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
realnum AccelLine
Definition: wind.h:61
size_t size(void) const
Definition: transition.h:331
realnum EnergyErg() const
Definition: transition.h:90
double * OpacStack
Definition: opacity.h:164
double * opacity_abs
Definition: opacity.h:104
qList st
Definition: iso.h:482
const int ipHE_LIKE
Definition: iso.h:65
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
t_opac opac
Definition: opacity.cpp:5
realnum Heiles_Troland_F
Definition: doppvel.h:37
realnum ** flux
Definition: rfield.h:68
t_Heavy Heavy
Definition: heavy.cpp:5
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:311
long int ipPradMax_nzone
Definition: pressure.h:106
double EnergyIonization
Definition: phycon.h:41
realnum * outlin_noplot
Definition: rfield.h:189
long int IonHigh[LIMELM+1]
Definition: dense.h:130
double PresRamCurr
Definition: pressure.h:39
realnum HLineWidth
Definition: hydrogenic.h:100
char TorF(bool l)
Definition: cddefines.h:749
t_magnetic magnetic
Definition: magnetic.cpp:17
bool lgMETALS
Definition: dynamics.h:92
bool lgTimeDependentStatic
Definition: dynamics.h:102
t_conv conv
Definition: conv.cpp:5
realnum PresMax
Definition: pressure.h:96
TransitionList HFLines("HFLines",&AnonStates)
t_phycon phycon
Definition: phycon.cpp:6
double EnthalpyDensity
Definition: phycon.h:50
realnum TurbVel
Definition: doppvel.h:21
double fudge(long int ipnt)
Definition: service.cpp:555
realnum AccelElectron
Definition: wind.h:58
double pressure
Definition: magnetic.h:41
realnum ** outlin
Definition: rfield.h:189
FILE * ioQQQ
Definition: cddefines.cpp:7
double * opacity_sct
Definition: opacity.h:107
long int nzone
Definition: cddefines.cpp:14
t_DoppVel DoppVel
Definition: doppvel.cpp:5
t_dynamics dynamics
Definition: dynamics.cpp:42
TransitionList TauLine2("TauLine2",&AnonStates)
#define MIN2(a, b)
Definition: cddefines.h:803
double PresTotlCurr
Definition: pressure.h:46
double anu(size_t i) const
Definition: mesh.h:120
long int nSpecies
Definition: taulines.cpp:22
bool lgISO[NISO]
Definition: dynamics.h:89
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:32
realnum PressureErrorAllowed
Definition: conv.h:268
double pres_radiation_lines_curr
Definition: pressure.h:61
t_dense dense
Definition: global.cpp:15
long int ipPradMax_line
Definition: pressure.h:103
static t_ADfA & Inst()
Definition: cddefines.h:209
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void PresTotCurrent(void)
realnum SmallA
Definition: iso.h:391
double EnthalpyDensity
Definition: magnetic.h:38
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int iteration
Definition: cddefines.cpp:16
bool lgSearch
Definition: conv.h:168
t_trace trace
Definition: trace.cpp:5
double drad
Definition: radius.h:31
bool lgBallistic(void) const
Definition: wind.h:31
void broken(void)
Definition: service.cpp:1067
realnum AccelGravity
Definition: wind.h:49
realnum AccelCont
Definition: wind.h:55
long int IonLow[LIMELM+1]
Definition: dense.h:129
realnum pden
Definition: dense.h:108
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
#define POW2
Definition: cddefines.h:979
const int ipH1s
Definition: iso.h:29
bool lgTrace
Definition: trace.h:12
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition: pressure.cpp:50
double RT_line_driving(void)
realnum DynaFlux(double depth)
Definition: dynamics.cpp:1292
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
realnum * ConInterOut
Definition: rfield.h:154
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
realnum AccelTotalOutward
Definition: wind.h:52
bool lgPres_magnetic_ON
Definition: pressure.h:91
long int iopcom
Definition: opacity.h:223
vector< diatomics * > diatoms
Definition: h2.cpp:8
t_hydro hydro
Definition: hydrogenic.cpp:5
#define cdEXIT(FAIL)
Definition: cddefines.h:482
realnum RadBetaMax
Definition: pressure.h:96
string chLineRadPres
Definition: pressure.h:109
realnum GetDopplerWidth(realnum massAMU)
bool lgPradCap
Definition: pressure.h:113
double PresGasCurr
Definition: pressure.h:46
long nWindLine
Definition: cdinit.cpp:19
double PresTotlError
Definition: pressure.h:46
t_radius radius
Definition: radius.cpp:5
realnum pbeta
Definition: pressure.h:96
long int nTotalIoniz
Definition: conv.h:159
bool lgDoLineTrans
Definition: rfield.h:99
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
const int ipH2p
Definition: iso.h:31
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
#define ASSERT(exp)
Definition: cddefines.h:613
char chDenseLaw[5]
Definition: dense.h:176
double EnergyBinding
Definition: phycon.h:54
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
bool lgLineRadPresOn
Definition: pressure.h:117
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
realnum xMassDensity
Definition: dense.h:101
vector< species > dBaseSpecies
Definition: taulines.cpp:15
double eden
Definition: dense.h:201
double EnergyExcitation
Definition: phycon.h:47
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
bool lgTurb_pressure
Definition: doppvel.h:42
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
bool lgPres_ram_ON
Definition: pressure.h:92
bool lgStatic(void) const
Definition: wind.h:24
void SumDensities(void)
Definition: dense.cpp:235
double DiskRadius
Definition: wind.h:78
double PresTurbCurr
Definition: pressure.h:42
#define fixit(a)
Definition: cddefines.h:417
bool lgElemsConserved(void)
Definition: dense.cpp:119
double fnzone
Definition: cddefines.cpp:15
void ShowMe(void)
Definition: service.cpp:205
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:349
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
Definition: rt_escprob.cpp:946
long int nflux
Definition: rfield.h:46
bool lgPres_radiation_ON
Definition: pressure.h:90
long int numLevels_local
Definition: iso.h:529
long int nCollapsed_local
Definition: iso.h:519
realnum windv
Definition: wind.h:18
realnum comass
Definition: wind.h:14
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
realnum fmul
Definition: wind.h:65
void Magnetic_evaluate(void)
Definition: magnetic.cpp:39