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