cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
pressure_change.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 /*PressureChange called by ConvPresTempEdenIoniz
4  * evaluate the current pressure, density change needed to get it to converge,
5  * applies this correction factor to all gas constituents,
6  * sets conv.lgConvPres true if good pressure, false if pressure change capped */
7 /*lgConvPres finds amount by which density should change to move towards pressure equilibrium
8  * returns true if pressure is converged */
9 #include "cddefines.h"
10 
11 #include "pressure_change.h"
12 
13 #include "colden.h"
14 #include "conv.h"
15 #include "cosmology.h"
16 #include "dark_matter.h"
17 #include "dense.h"
18 #include "dynamics.h"
19 #include "geometry.h"
20 #include "phycon.h"
21 #include "pressure.h"
22 #include "radius.h"
23 #include "struc.h"
24 #include "thermal.h"
25 #include "trace.h"
26 #include "wind.h"
27 
28 /* zoneDensity finds density where prescribed */
29 double zoneDensity()
30 {
31  double new_density;
32 
33  DEBUG_ENTRY( "zoneDensity()" );
34 
36 
37  /* evaluate a series of possible options, and set new_density */
38  /* inside out globule */
39  if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
40  {
41  /* GLBDST is distance from globule, or glbrad-DEPTH */
42  if( radius.glbdst < 0. )
43  {
44  fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
45  fprintf( ioQQQ, " This is routine lgConvPres, GLBDST is%10.2e\n",
46  radius.glbdst );
48  }
49  new_density =
52  }
53  else if( cosmology.lgDo )
54  {
55  /* cosmological - density varies because of expansion of universe */
56  double dnew = GetDensity( cosmology.redshift_current );
57  new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
58  }
59  else if( (strcmp(dense.chDenseLaw,"WIND") == 0) ) {
60 
61  if ( wind.lgStatic() )
62  {
63  fprintf( ioQQQ, " PROBLEM WIND called with zero velocity - this is impossible.\n Sorry.\n" );
64  TotalInsanity();
65  }
66 
67  /* this is positive wind velocity the outflowing wind beyond sonic point */
68  else if( wind.lgBallistic() )
69  {
70 
71  /* this is logic for supersonic wind solution,
72  * well above sonic point. */
73  if( nzone > 1 )
74  {
75  /* Wind model */
76 
77  if( trace.lgTrace && trace.lgWind )
78  {
79  fprintf(ioQQQ," lgConvPres sets AccelGravity %.3e lgDisk?%c\n",
81  TorF(wind.lgDisk) );
82  }
83 
84  /* following is form of constant acceleration equation, v^2 = u^2 + 2 a s
85  * struc.windv[nzone-2] is velocity of previous zone
86  * this increments that velocity to form square of new wind
87  * velocity for outer edge of this zone */
88  double term = POW2(struc.windv[nzone-2]) + 2.*(wind.AccelTotalOutward - wind.AccelGravity)* radius.drad;
89 
90  /* increment velocity if it is substantially positive */
91  fixit("RHS of comparison should be ~ sound speed squared");
92  if( term <= 1e3 )
93  {
94  /* wind velocity is well below sonic point, give up,
95  * do not change velocity */
96  wind.lgVelPos = false;
97  }
98  else
99  {
100  /* wind.windv is velocity at OUTER edge of this zone */
101  term = sqrt(term);
102  if (wind.windv > 0)
103  {
104  wind.windv = (realnum) term;
105  }
106  else
107  {
108  wind.windv = -(realnum) term;
109  }
110  wind.lgVelPos = true;
111  }
112 
113  if( trace.lgTrace && trace.lgWind )
114  {
115  fprintf(ioQQQ," lgConvPres new wind V zn%li %.3e AccelTotalOutward %.3e AccelGravity %.3e\n",
117  }
118  }
119 
120  /* conservation of mass sets density here */
121  new_density = wind.emdot/(wind.windv*radius.r1r0sq) *
123  }
124  else
125  {
126  fprintf(ioQQQ,"chDenseLaw==\"WIND\" must now be ballistic or static\n");
127  TotalInsanity();
128  }
129  }
130 
131  else if( strcmp(dense.chDenseLaw,"SINE") == 0 )
132  {
133  /* rapid density fluctuation */
134  if( dense.lgDenFlucRadius )
135  {
136  new_density =
139  }
140  else
141  {
142  new_density =
145  }
146  }
147 
148  else if( strcmp(dense.chDenseLaw,"POWR") == 0 )
149  {
150  /* power law function of radius */
151  double dnew = dense.den0*pow(radius.Radius/radius.rinner,(double)dense.DensityPower);
152  new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
153  }
154 
155  else if( strcmp(dense.chDenseLaw,"POWD") == 0 )
156  {
157  /* power law function of depth */
158  double dnew = dense.den0*pow(1. + radius.depth/dense.rscale,(double)dense.DensityPower);
159  new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
160  }
161 
162  else if( strcmp(dense.chDenseLaw,"POWC") == 0 )
163  {
164  /* power law function of column density */
165  double dnew = dense.den0*pow(1.f + colden.colden[ipCOL_HTOT]/
167  new_density = dnew*(scalingDensity()/dense.gas_phase[ipHYDROGEN]);
168  }
169 
170 
171  else if( strncmp( dense.chDenseLaw ,"DLW" , 3) == 0 )
172  {
173  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
174  {
175  /* call ACF sub */
176  new_density =
179  }
180  else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
181  {
182  /* call table interpolation subroutine
183  * >>chng 96 nov 29, added dense_tabden */
184  new_density = dense.DLW.tabval(radius.Radius,radius.depth) *
186  }
187  else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
188  {
189  /* call parametrized wind subroutine */
190  new_density = dense_parametric_wind(radius.Radius) *
192  }
193  else
194  {
195  fprintf( ioQQQ, " Insanity, lgConvPres gets chCPres=%4.4s\n",
196  dense.chDenseLaw );
198  }
199 
200  if( new_density/scalingDensity() > 3. || new_density/scalingDensity()< 1./3 )
201  {
202  static bool lgWARN2BIG=false;
203  if( !lgWARN2BIG )
204  {
205  lgWARN2BIG = true;
206  fprintf(ioQQQ,"\n\n >========== Warning! The tabulated or functional change in density as a function of depth was VERY large. This is zone %li.\n",nzone);
207  fprintf(ioQQQ," >========== Warning! This will cause convergence problems. \n");
208  fprintf(ioQQQ," >========== Warning! The current radius is %.3e. \n",radius.Radius);
209  fprintf(ioQQQ," >========== Warning! The current depth is %.3e. \n",radius.depth);
210  fprintf(ioQQQ," >========== Warning! Consider using a more modest change in density vs radius. \n\n\n");
211  }
212  }
213  }
214 
215  else if( strcmp(dense.chDenseLaw,"CDEN") == 0 )
216  {
217  /* this is the default, constant density */
218  new_density = scalingDensity();
219  }
220 
221  else
222  {
223  fprintf( ioQQQ, " Unknown density law=%s= This is a major internal error.\n",
224  dense.chDenseLaw );
225  ShowMe();
227  }
228 
229  return new_density;
230 }
231 
232 enum {
241 };
242 
244 STATIC double stepDensity(const PresMode & presmode, solverState &st);
245 
247 {
248  /* >>chng 04 feb 11, add option to remember current density and pressure */
252 }
253 
254 STATIC bool lgTestPressureConvergence( double new_density)
255 {
256  double density_change_factor = new_density/scalingDensity();
257 
258  /* now see whether current pressure is within error bounds */
259  if( density_change_factor > 1. + conv.PressureErrorAllowed ||
260  density_change_factor < 1. - conv.PressureErrorAllowed )
261  {
262  return false;
263  }
264  else
265  {
266  return true;
267  }
268 }
269 
270 STATIC double limitedDensityScaling( double new_density, double dP_chng_factor )
271 {
272  double density_change_factor = new_density/scalingDensity()-1.;
273 
274  /* dP_chng_factor is initially 1, becomes smaller if sign of pressure change, changes */
276 
277  /* make sure that change is not too extreme */
278  density_change_factor *= dP_chng_factor;
279  density_change_factor = MIN2(density_change_factor,+pdelta);
280  density_change_factor = MAX2(density_change_factor,-pdelta);
281  return 1.+density_change_factor;
282 }
283 
284 /*PressureChange evaluate the current pressure, and change needed to
285  * get it to PresTotlInit,
286  * return value is true if density was changed, false if no changes were necessary */
288  /* this is change factor, 1 at first, becomes smaller as oscillations occur */
289  double dP_chng_factor, const PresMode & presmode, solverState &st,
290  bool& lgAbort, bool &lgStable )
291 {
292  DEBUG_ENTRY( "PressureChange()" );
293 
294  lgAbort = false;
295 
296  /* first evaluate total pressure for this location, and current conditions
297  * CurrentPressure is just sum of gas plus local line radiation pressure */
298  /* this sets values of pressure.PresTotlCurr, also wind velocity for dynamics */
299  PresTotCurrent();
301 
302  double old_density = scalingDensity();
303  double old_pressure = pressure.PresTotlCurr;
304 
305  double new_density = stepDensity(presmode, st);
306 
308 
309  double density_change_factor = limitedDensityScaling(new_density, dP_chng_factor);
310 
311  {
312  /*@-redef@*/
313  enum{DEBUG_LOC=false};
314  static long int nsave=-1;
315  /*@+redef@*/
316  if( DEBUG_LOC /*&& nzone > 150 && iteration > 1*/ )
317  {
318  if( nsave-nzone ) fprintf(ioQQQ,"\n");
319  nsave = nzone;
320  fprintf(ioQQQ,"nnzzone\t%li\t%.2f%%\t%.3f\t%.2e\t%.2e\t%.2e\n",
321  nzone,
322  density_change_factor,
323  /* when this is negative we need to raise the density */
324  pressure.PresTotlError*100.,
328  }
329  }
330 
331  if( trace.lgTrace )
332  {
333  fprintf( ioQQQ,
334  " PressureChange called, changing HDEN from %10.3e to %10.3e Set fill fac to %10.3e\n",
335  dense.gas_phase[ipHYDROGEN], density_change_factor*dense.gas_phase[ipHYDROGEN],
336  geometry.FillFac );
337  }
338 
339  bool lgChange = ( density_change_factor != 1. );
340 
341  lgStable = true;
342  if( lgChange )
343  {
344  ScaleAllDensities((realnum) density_change_factor);
345 
346  /* must call TempChange since ionization has changed, there are some
347  * terms that affect collision rates (H0 term in electron collision) */
348  TempChange(phycon.te , false);
349 
350  /* heating cooling balance while doing ionization,
351  * this is where the heavy lifting is done, this calls PresTotCurrent,
352  * which sets pressure.PresTotlCurr */
353  int lgTempStatus = ConvTempEdenIoniz();
354  if( lgTempStatus != 0 )
355  {
356  lgAbort = true;
357  }
358  double new_pressure = pressure.PresTotlCurr;
359  double dpdrho = (new_pressure-old_pressure)/(new_density-old_density);
360  double dpdrhoScale = (new_pressure+old_pressure)/(new_density+old_density);
361  if ( presmode.zone == SUBSONIC)
362  {
363  if ( dpdrho < -0.01*dpdrhoScale )
364  lgStable = false;
365  }
366  else if ( presmode.zone == SUPERSONIC)
367  {
368  if ( dpdrho > 0.01*dpdrhoScale )
369  lgStable = false;
370  }
371  if (0)
372  fprintf(ioQQQ,"ConvPres STABLE??? nz %ld od %.2e nd %.2e eps %.2e op %.2e np %.2e dpdn %.2e cmp %.2e %c\n",
373  nzone,old_density,new_density,new_density/old_density-1.0,
374  old_pressure,new_pressure,dpdrho,dpdrhoScale,TorF(lgStable));
375  if (0)
376  if (!lgStable)
377  {
378  conv.lgConvPres = false; // If evidence of instability found, ensure at least another iteration
379  fixit("Disabled test on stability");
380  // log zone warning when pressure instability diagnosed
381  // possibly reduce zone stepsize when close
382  }
383  }
384 
385  {
386  /*@-redef@*/
387  enum {DEBUG_LOC=false};
388  /*@+redef@*/
389  if( DEBUG_LOC && (nzone>215) )
390  {
391  fprintf( ioQQQ,
392  "%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%c\n",
393  radius.depth,
400  /* subtract continuum rad pres which has already been added on */
402  wind.windv/1e5,
403  sqrt(5.*pressure.PresGasCurr/3./dense.xMassDensity)/1e5,
404  TorF(conv.lgConvPres) );
405  }
406  }
407 }
408 
409 /* ============================================================================== */
410 /* DynaPresChngFactor, called from PressureChange to evaluate factor needed
411  * to find new density needed for
412  * current conditions and wind solution, returns ratio of new to old density */
413 
414 /* object is to get the local ram pressure
415  * RamNeeded = pressure.PresTotlInit + pressure.PresGasCurr + pressure.PresInteg;
416  * to equal the sum of the inital pressur at the illuminated face, the local gas pressure,
417  * and the integrate radiative acceleration from the incident continuum
418  *
419  * the local gas pressure is linear in density if the temperature is constant,
420  *
421  * the local ram pressure is inversely linear in density because of the relationship
422  * between velocity and density introduced by the mass flux conservation
423  */
424 
425 /* stepDensity finds a density which should move towards pressure equilibrium
426  * sets pressure.PresTotlError */
427 STATIC double stepDensity(const PresMode &presmode, solverState &st)
428 {
429  DEBUG_ENTRY( "stepDensity()" );
430 
431  double PresTotlCorrect = st.press;
432 
433  double densityCurrent = scalingDensity();
434 
435  double er = pressure.PresTotlCurr-PresTotlCorrect;
436  pressure.PresTotlError = er/PresTotlCorrect;
437 
438  if( trace.lgTrace && presmode.global != CPRE )
439  {
440  fprintf( ioQQQ,
441  " DynaPresChngFactor set PresTotlCorrect=%.3e PresTotlInit=%.3e PresInteg=%.3e DivergePresInteg=%.3e\n",
444  }
445 
446  if( dynamics.lgTracePrint && presmode.global != CPRE)
447  fprintf(ioQQQ,"Pressure: init %g +rad %g +diverge %g = %g cf %g\n",
449  dynamics.DivergePresInteg, PresTotlCorrect, pressure.PresTotlCurr);
450  /*fprintf(ioQQQ,"DEBUG\t%.2f\thden\t%.4e\tPtot\t%.4e\tPgas\t%.4e\n",
451  fnzone, scalingDensity(),pressure.PresTotlCurr,pressure.PresGasCurr );*/
452 
453  double rhohat, dpdrho;
454  if( presmode.global == CPRE || presmode.global == ORIGINAL ||
455  st.lastzone != nzone || fabs(er-st.erp) < SMALLFLOAT
456  || fabs(densityCurrent-st.dp) < SMALLFLOAT)
457  {
458  if ( presmode.zone == SUBSONIC )
459  {
460  dpdrho = pressure.PresTotlCurr/densityCurrent;
461  }
462  else
463  {
464  dpdrho = -pressure.PresTotlCurr/densityCurrent;
465  }
466  rhohat = densityCurrent - er/dpdrho;
467  }
468  else
469  {
470  /* second iteration on this zone, do linear fit to previous Pres - rho curve
471  * Linear approximation to guess root with two independent samples */
472  dpdrho = (st.erp-er)/(st.dp-densityCurrent);
473  rhohat = densityCurrent - er/dpdrho;
474 
475  /* Subsonic case: pressure ^ with density ^ => increase density further */
476  /* Super " case: pressure ^ with density v => decrease density further */
477 
478  /* Force the solution towards the required branch when it
479  * appears to have the "wrong" value of dP/drho */
480  if(presmode.zone == SUBSONIC && dpdrho <= 0)
481  {
482  if (er < 0.0)
483  rhohat = 1.03*densityCurrent;
484  else
485  rhohat = 0.97*densityCurrent;
486  }
487  else if(presmode.zone == SUPERSONIC && dpdrho >= 0)
488  {
489  if (er > 0.0)
490  rhohat = 1.03*densityCurrent;
491  else
492  rhohat = 0.97*densityCurrent;
493  }
494  }
495 
496  st.dp = densityCurrent;
497  st.erp = er;
498  st.lastzone = nzone;
499 
500  if( presmode.global != CPRE && dynamics.lgTracePrint )
501  fprintf(ioQQQ,"windv %li r %g v %g f %g\n",
503 
504  /* convergence trace at this level */
505  if( trace.nTrConvg>=1 )
506  {
507  if ( presmode.global == CPRE )
508  {
509  fprintf( ioQQQ,
510  " PresChng %s density old %.3e new %.3e current %.3e correct %.3e dpdn %.3e\n",
511  dynamics.chPresMode , densityCurrent, rhohat , pressure.PresTotlCurr,
512  PresTotlCorrect, dpdrho
513  );
514  }
515  else
516  {
517  fprintf( ioQQQ,
518  " DynaPresChngFactor mode %s new density %.3f vel %.3e\n",
519  dynamics.chPresMode , rhohat , wind.windv );
520  }
521  }
522 
523  return rhohat;
524 }
525 
527 {
528  DEBUG_ENTRY( "PresMode::set()" );
529  /* dynamics.lgSetPresMode is flag to indicate sane value of dynamics.chPresMode.
530  * If set true with SET DYNAMICS PRESSURE MODE
531  * then do not override that choice */
532  if( !dynamics.lgSetPresMode )
533  {
534  /* above set true if pressure mode was set with
535  * SET DYNAMICS PRESSURE MODE - if we got here
536  * it was not set, and must make a guess */
538  strcpy( dynamics.chPresMode , "supersonic" );
539  else
540  strcpy( dynamics.chPresMode , "subsonic" );
541  /* clear the flag - pressure mode has been set */
542  dynamics.lgSetPresMode = true;
543  }
544 
545  /* if globally looking for transonic solution, then locally sometimes
546  * supersonic, sometimes subsonic - this branch sets global flag,
547  * which can also be set with SET DYNAMICS PRESSURE MODE.
548  * Under default conditions, ChPresMode was set in previous branch
549  * to sub or super sonic depending on current velocity on first time*/
550 
551  if (strcmp(dense.chDenseLaw,"CPRE") != 0)
552  {
553  ASSERT (strcmp(dense.chDenseLaw,"DYNA") == 0);
554  }
555 
556  if (strcmp(dense.chDenseLaw,"CPRE") == 0)
557  {
558  global = CPRE;
560  }
561  else if( strcmp( dynamics.chPresMode , "original" ) == 0 )
562  {
563  global = ORIGINAL;
565  }
566  else if( strcmp( dynamics.chPresMode , "subsonic" ) == 0 )
567  {
568  global = SUBSONIC;
570  }
571  /* supersonic */
572  else if( strcmp( dynamics.chPresMode , "supersonic" ) == 0 )
573  {
574  global = SUPERSONIC;
576  }
577  /* strong d */
578  else if( strcmp( dynamics.chPresMode , "strongd" ) == 0 )
579  {
580  global = STRONGD;
582  }
583  else if( strcmp( dynamics.chPresMode , "shock" ) == 0 )
584  {
585  global = SHOCK;
587  }
588  else if( strcmp( dynamics.chPresMode , "antishock-by-mach" ) == 0 )
589  {
590  global = ANTISHOCK2;
592  }
593  else if( strcmp( dynamics.chPresMode , "antishock" ) == 0 )
594  {
595  global = ANTISHOCK;
597  }
598 
599  /* this sets pressure mode for the current zone based on global mode
600  * and local conditions */
601  if(global == CPRE)
602  {
603  zone = SUBSONIC;
604  }
605  else if(global == ORIGINAL)
606  {
607  /* in this mode use comparison between ram and gas pressure to determine
608  * whether sub or super sonic */
610  {
611  zone = SUBSONIC;
612  }
613  else
614  {
615  zone = SUPERSONIC;
616  }
617  }
618  else if(global == STRONGD)
619  {
620  //if(nzone <= 1)
621  zone = SUPERSONIC;
622  }
623  else if(global == SUBSONIC)
624  {
625  zone = SUBSONIC;
626  }
627  else if(global == SUPERSONIC)
628  {
629  zone = SUPERSONIC;
630  }
631  else if(global == SHOCK)
632  {
634  {
635  zone = SUBSONIC;
636  }
637  else
638  {
639  zone = SUPERSONIC;
640  }
641  }
642  else if(global == ANTISHOCK)
643  {
645  {
646  zone = SUPERSONIC;
647  }
648  else
649  {
650  zone = SUBSONIC;
651  }
652  }
653  else if(global == ANTISHOCK2)
654  {
655  /* WJH 19 May 2004: This version of the antishock mode will
656  insert the antishock at the point where the isothermal Mach
657  number falls below a certain value, dynamics.ShockMach */
659  {
660  zone = SUPERSONIC;
661  }
662  else
663  {
664  zone = SUBSONIC;
665  }
666  }
667  else
668  {
669  printf("Need to set global pressure mode\n");
670  cdEXIT( EXIT_FAILURE );
671  }
672 }
673 
674 double pressureZone(const PresMode &presmode)
675 {
676  DEBUG_ENTRY( "pressureZone()" );
677  double press;
678  if( presmode.global == CPRE )
679  {
680  /* constant pressure */
682  {
683  /* >>chng 01 oct 31, replace pneed with CorretPressure */
684  /* this has pressure due to incident continuum */
686  }
687  else
688  {
689  /* this does not have pressure due to incident continuum*/
690  press = pressure.PresTotlInit*
691  /* following term normally unity, power law set with option par on cmmnd*/
693  }
694 
696  {
697  fixit("check correct zone pressure for NFW");
698  // doesn't this exclude current zone gravity pressure.RhoGravity?
699  // and doesn't the above exclude current rad press, pressure.pinzon?
700  // That would be a second order correction, which would need
701  // a consistent second order treatment elsewhere to make it
702  // worthwhile.
703  press += pressure.IntegRhoGravity;
704  }
705  }
706  else
707  {
708  /* update the current desired pressure */
711  }
712  return press;
713 }
714 
realnum den0
Definition: dense.h:252
double PresTotlInit
Definition: pressure.h:52
bool lgContRadPresOn
Definition: pressure.h:65
vector< double > hist_pres_density
Definition: conv.h:291
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
STATIC void logPressureState()
realnum GetDensity(realnum z)
Definition: cosmology.cpp:28
t_colden colden
Definition: colden.cpp:5
double MaxFractionalDensityStepPerIteration
Definition: conv.h:270
double pressureZone(const PresMode &presmode)
realnum DensityPower
Definition: dense.h:250
realnum PresInteg
Definition: pressure.h:69
realnum flcPhase
Definition: dense.h:265
realnum csecnd
Definition: dense.h:264
realnum * windv
Definition: struc.h:25
bool lgSonicPointAbortOK
Definition: pressure.h:125
vector< double > hist_pres_error
Definition: conv.h:291
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_struc struc
Definition: struc.cpp:6
void PressureChange(double dP_chng_factor, const PresMode &presmode, solverState &st, bool &lgAbort, bool &lgStable)
const realnum SMALLFLOAT
Definition: cpu.h:246
bool lgDisk
Definition: wind.h:74
double PresRamCurr
Definition: pressure.h:39
char TorF(bool l)
Definition: cddefines.h:749
char chPresMode[20]
Definition: dynamics.h:120
double IntegRhoGravity
Definition: pressure.h:83
bool lgConvPres
Definition: conv.h:192
bool lgDenFlucRadius
Definition: dense.h:259
bool lgTracePrint
Definition: dynamics.h:183
double ShockMach
Definition: dynamics.h:127
t_conv conv
Definition: conv.cpp:5
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum FillFac
Definition: geometry.h:29
long int nzone
Definition: cddefines.cpp:14
STATIC bool lgTestPressureConvergence(double new_density)
t_dynamics dynamics
Definition: dynamics.cpp:42
#define MIN2(a, b)
Definition: cddefines.h:803
bool lgDo
Definition: cosmology.h:44
double PresTotlCurr
Definition: pressure.h:46
bool lgVelPos
Definition: wind.h:71
double zoneDensity()
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:32
realnum PressureErrorAllowed
Definition: conv.h:268
bool lgWind
Definition: trace.h:106
realnum glbpow
Definition: radius.h:134
int ConvTempEdenIoniz(void)
double pres_radiation_lines_curr
Definition: pressure.h:61
t_dense dense
Definition: global.cpp:15
void PresTotCurrent(void)
realnum glbdst
Definition: radius.h:134
Wind wind
Definition: wind.cpp:5
t_trace trace
Definition: trace.cpp:5
double drad
Definition: radius.h:31
bool lgBallistic(void) const
Definition: wind.h:31
double rinner
Definition: radius.h:31
realnum AccelGravity
Definition: wind.h:49
realnum pinzon
Definition: pressure.h:69
t_geometry geometry
Definition: geometry.cpp:5
void ScaleAllDensities(realnum factor)
Definition: dense.cpp:56
t_dark_matter dark
Definition: dark_matter.cpp:5
realnum redshift_current
Definition: cosmology.h:26
#define POW2
Definition: cddefines.h:979
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
realnum DynaFlux(double depth)
Definition: dynamics.cpp:1292
vector< double > hist_pres_current
Definition: conv.h:291
t_pressure pressure
Definition: pressure.cpp:9
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
realnum AccelTotalOutward
Definition: wind.h:52
realnum flong
Definition: dense.h:262
#define cdEXIT(FAIL)
Definition: cddefines.h:482
DepthTable DLW
Definition: dense.h:198
double tabval(double r0, double depth) const
Definition: depth_table.cpp:8
int gravity_symmetry
Definition: pressure.h:84
double PresGasCurr
Definition: pressure.h:46
double PresTotlError
Definition: pressure.h:46
int nTrConvg
Definition: trace.h:27
t_radius radius
Definition: radius.cpp:5
realnum PresPowerlaw
Definition: pressure.h:36
realnum gas_phase[LIMELM]
Definition: dense.h:76
realnum glbrad
Definition: radius.h:134
realnum rscale
Definition: dense.h:251
STATIC double stepDensity(const PresMode &presmode, solverState &st)
#define ASSERT(exp)
Definition: cddefines.h:613
char chDenseLaw[5]
Definition: dense.h:176
bool lgSetPresMode
Definition: dynamics.h:172
realnum emdot
Definition: wind.h:39
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
t_cosmology cosmology
Definition: cosmology.cpp:8
realnum scalingDensity(void)
Definition: dense.cpp:409
STATIC double limitedDensityScaling(double new_density, double dP_chng_factor)
realnum xMassDensity
Definition: dense.h:101
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
double dense_fabden(double radius, double depth)
double ShockDepth
Definition: dynamics.h:123
double dense_parametric_wind(double rad)
bool lgStatic(void) const
Definition: wind.h:24
#define fixit(a)
Definition: cddefines.h:417
double r1r0sq
Definition: radius.h:31
realnum glbden
Definition: radius.h:134
void ShowMe(void)
Definition: service.cpp:205
double te
Definition: phycon.h:21
realnum cfirst
Definition: dense.h:263
const int ipHYDROGEN
Definition: cddefines.h:349
realnum colden[NCOLD]
Definition: colden.h:32
realnum windv
Definition: wind.h:18
bool lgAbort
Definition: cddefines.cpp:10
realnum DivergePresInteg
Definition: dynamics.h:177