cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
radius_next.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*radius_next use adaptive logic to find next zone thickness */
4 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
5 /*GrainRateDr called by radius_next to find grain heating rate dr */
14 #include "cddefines.h"
15 #include "iso.h"
16 #include "geometry.h"
17 #include "h2.h"
18 #include "hyperfine.h"
19 #include "opacity.h"
20 #include "dense.h"
21 #include "heavy.h"
22 #include "grainvar.h"
23 #include "elementnames.h"
24 #include "dynamics.h"
25 #include "thermal.h"
26 #include "hmi.h"
27 #include "coolheavy.h"
28 #include "timesc.h"
29 #include "stopcalc.h"
30 #include "colden.h"
31 #include "phycon.h"
32 #include "rt.h"
33 #include "trace.h"
34 #include "wind.h"
35 #include "save.h"
36 #include "pressure.h"
37 #include "iterations.h"
38 #include "struc.h"
39 #include "radius.h"
40 #include "dark_matter.h"
41 #include "mole.h"
42 #include "rfield.h"
43 #include "freebound.h"
44 #include "lines_service.h"
45 #include "cosmology.h"
46 
47 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
48 STATIC void ContRate(double *freqm,
49  double *opacm);
50 
51 /*GrainRateDr called by radius_next to find grain heating rate dr */
52 STATIC void GrainRateDr(double *grfreqm,
53  double *gropacm);
54 
56 {
57  string m_str;
58  bool *m_flag;
59 public:
60  drChoiceItem(const string& str, bool *flag) : m_str(str), m_flag(flag) {}
61  void select() const
62  {
63  if (m_flag)
64  *m_flag = true;
65  }
66  const char* c_str() const
67  {
68  return m_str.c_str();
69  }
70 };
71 
72 class drList
73 {
74  multimap<double,drChoiceItem> m_map;
75 public:
76  typedef multimap<double,drChoiceItem>::const_iterator const_iterator;
77  void insert(double val, const string& str, bool* flag = NULL)
78  {
79  if (flag != NULL)
80  *flag = false;
81  m_map.insert(pair<const double,drChoiceItem>(val,drChoiceItem(str,flag)));
82  }
83  void clear()
84  {
85  m_map.clear();
86  }
88  {
89  return m_map.begin();
90  }
92  {
93  return m_map.end();
94  }
95  size_t size() const
96  {
97  return m_map.size();
98  }
99 };
100 
101 
102 /*radius_next use adaptive logic to find next zone thickness
103  * return 0 if ok, 1 for abort */
104 int radius_next(void)
105 {
106  static double OHIIoHI,
107  OldHeat = 0.,
108  OldTe = 0.,
109  OlddTdStep = 0.,
110  OldWindVelocity = 0.,
111  Old_H2_heat_cool;
112 
113  const double BigRadius = 1e30;
114 
115  DEBUG_ENTRY( "radius_next()" );
116 
117  /*--------------------------------------------------------------
118  *
119  * this sub determines the thickness of the next zone
120  * if is called one time for each zone
121  *
122  *-------------------------------------------------------------- */
123 
124  if( trace.lgTrace )
125  fprintf( ioQQQ, " radius_next called\n" );
126 
127  /* >>chng 03 sep 21 - decide whether this is the first call */
128  bool lgFirstCall = ( nzone == 0 );
129 
130  drList drChoice;
131 
132  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
133  * check on change in hydrogen ionization */
134  if( !lgFirstCall && OHIIoHI > 1e-15 &&
135  (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) )
136  {
137  double atomic_frac = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
138  double drHion;
139  /* ratio of current HII/HI to old value - < 1 when becoming more neutral */
140  /* this is now change in atomic fraction */
141  double x = fabs(1. - atomic_frac/OHIIoHI);
142  if( atomic_frac > 0.05 && atomic_frac < 0.9 )
143  {
144  /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
145  * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
146  /* >>chng 02 jun 05 from 0.2 to 0.05 poorly resolved i-front, also added two-branch logic*/
147  drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
148  }
149  else
150  {
151  /* >>chng 96 feb 23 from 0.7 to 0.3 because punched through H i-front
152  * >>chng 97 may 5, from 0.3 to 0.2 since punched through H i-front */
153  drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
154  }
155  double SaveOHIIoHI = OHIIoHI;
157  ostringstream oss;
158  oss << "change in H ionization old=" << scientific << setprecision(3);
159  oss << SaveOHIIoHI << " new=" << OHIIoHI;
160  drChoice.insert( drHion, oss.str() );
161  }
162  else
163  {
164  if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
166  else
167  OHIIoHI = 0.;
168  }
169 
170  if( rt.dTauMase < -0.01 )
171  {
172  /* maser so powerful, must limit inc in tay
173  * >>chng 96 aug 08, from 0.3 to 0.1 due to large maser crash */
174  double drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
175  ostringstream oss;
176  //
177  // NB NB -- DON'T ALTER THIS STRING, setting rt.lgMaserSetDR below keys from this!
178  // No longer true:
179  // >>chng 13 oct 01 rjrw take logical pointer argument for
180  // >>optional flag to be set if selected
181  //
182  oss << "H maser dTauMase=" << scientific << setprecision(2) << rt.dTauMase << " ";
183  oss << rt.mas_species << " " << rt.mas_ion << " " << rt.mas_hi << " " << rt.mas_lo;
184  drChoice.insert( drHMase, oss.str(), &rt.lgMaserSetDR );
185  }
186 
187  /* check change in relative ionization - this is to insure a gradual approach to
188  * ionization fronts - were dr allowed to remain constant we would crash through
189  * a front in a single zone */
190  double change_heavy_frac_big = -1.;
191  double frac_change_heavy_frac_big = -1.;
192  const realnum CHANGE_ION_HEAV = 0.2f;
193  const realnum CHANGE_ION_HHE = 0.15f;
194  if( nzone > 4 )
195  {
196  long ichange_heavy_nelem = -1L;
197  long ichange_heavy_ion = -1L;
198  double dr_change_heavy = BigRadius;
199  for( int nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
200  {
201  if( dense.lgElmtOn[nelem] )
202  {
203  realnum change;
204  /* will only track ions with fractional abundances greater than this
205  * He and C are special cases because of special roles in PDRs */
206  realnum frac_limit;
207  if( nelem<=ipHELIUM || nelem==ipCARBON )
208  {
209  frac_limit = 1e-4f;
210  change = CHANGE_ION_HHE;
211  }
212  else
213  {
214  /* struc.dr_ionfrac_limit set in init_coreload, is limit to fractional
215  * abundances of ions that are used in prt_comment to indicate oscillations
216  * or rapid change in ionization */
217  frac_limit = struc.dr_ionfrac_limit/2.f;
218  change = CHANGE_ION_HEAV;
219  }
220 
221  for( int ion=0; ion<=nelem+1; ++ion )
222  {
223  realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
224  realnum abundold = struc.xIonDense[nzone-3][nelem][ion]/SDIV( struc.gas_phase[nzone-3][nelem]);
225  if( abundnew > frac_limit && abundold > frac_limit )
226  {
227  /* this test is not done when [nzone-n] <0 */
228  realnum abundolder = struc.xIonDense[nzone-4][nelem][ion]/SDIV( struc.gas_phase[nzone-4][nelem]);
229  realnum abundoldest = struc.xIonDense[nzone-5][nelem][ion]/SDIV( struc.gas_phase[nzone-5][nelem]);
230  /* this is fractional change, use min of old and new abundances, to pick up
231  * rapid changing Ca+ sooner */
232  double change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
233  /* want fractional change to be less than this factor */
234  if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
235  /* >>chng 03 dec 07, add test that abund is not oscillating */
236  /* also test that abundance is increasing - we are headed into a front */
237  ((abundnew-abundold)>0.) &&
238  ((abundold-abundolder)>0.) &&
239  ((abundolder-abundoldest)>0.) )
240  {
241  ichange_heavy_nelem = nelem;
242  ichange_heavy_ion = ion;
243  change_heavy_frac_big = change_heavy_frac;
244  frac_change_heavy_frac_big = abundnew;
245  /* factor in max has been adjusted to prevent code from running
246  * through various ionization fronts.
247  * >>chng 03 dec 06, from min of 0.5 to min of 0.25, crash into He i-front
248  * in hizqso.in */
249  /* >>chng 04 mar 03, min had become 0.1, forced oscillations in nova.in
250  * in silicon, zoning changed greatly, causing change in diffuse lin
251  * pumping. put back to 0.25 */
252  dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
253  }
254  }
255  }
256  }
257  }
258  /* set to -1 before loop over ions, if no significant changes in ionization occurred
259  * then may still be -1
260  */
261  ostringstream oss;
262  if(ichange_heavy_nelem>=0)
263  {
264  oss << "change in ionization, element " << elementnames.chElementName[ichange_heavy_nelem];
265  oss << " ion (C scale) " << ichange_heavy_ion << " rel change " << scientific << setprecision(2);
266  oss << change_heavy_frac_big << " ion frac " << frac_change_heavy_frac_big;
267  }
268  else
269  {
270  oss << "none";
271  dr_change_heavy = BigRadius;
272  }
273  drChoice.insert( dr_change_heavy, oss.str() );
274  }
275 
276  /* if Leiden hacks are on then use increase in dust extinction as control
277  * >>chng 05 aug 13, add this */
279  {
280  /* Draine field is only defined over narrow range in FUV - must not let change
281  * in extinction become too large -
282  * prefactor is change in optical depth */
283  double drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point ) /
285  drChoice.insert( drLeiden_hack, "Leiden hack" );
286  }
287 
288  // limit to dust optical depth per zone
289  static double extin_mag_V_point_old = -1.;
290  if( nzone>1 )
291  {
292  const double extin_mag_V_limit = 20.+0.01*extin_mag_V_point_old;
293  double dExtin = (rfield.extin_mag_V_point - extin_mag_V_point_old)/extin_mag_V_limit;
294  if( dExtin > SMALLFLOAT )
295  {
296  double drExtintion = radius.drad / dExtin;
297  ostringstream oss;
298  oss << "change in V mag extinction; current=" << scientific << setprecision(3) <<
300  oss << fixed << " delta=" << dExtin;
301  drChoice.insert( drExtintion, oss.str() );
302  }
303 
304  }
305  extin_mag_V_point_old = rfield.extin_mag_V_point;
306 
307  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
308  /* check how heating is changing */
309  if( nzone > 1 && !thermal.lgTemperatureConstant &&
310  /* if density fluctuations in place then override change in heat for dr set */
311  !( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn ) )
312  {
313  double dHdStep = safe_div(fabs(thermal.htot-OldHeat),thermal.htot,0.0);
314  if( dHdStep > 0. )
315  {
316  double drdHdStep;
317  if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
318  {
319  drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
320  }
321  else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
322  {
323  drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
324  }
325  else
326  {
327  /* changed from .15 to .12 for outer edge of coolhii too steep dT
328  * changed to .10 for symmetry, big change in some rates, 95nov14
329  * changed from .10 to .125 since parispn seemed to waste zones
330  * >>chng 96 may 21, from .125 to .15 since pn's still waste zones */
331  drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
332  }
333  ostringstream oss;
334  oss << "change in heating; current=" << scientific << setprecision(3) << thermal.htot;
335  oss << fixed << " delta=" << dHdStep;
336  drChoice.insert( drdHdStep, oss.str() );
337  }
338  }
339  OldHeat = thermal.htot;
340 
341  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
342  /* pressure due to incident continuum if in eos */
343  if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
344  {
345  if( nzone > 2 && pressure.pinzon > 0. )
346  {
347  /* pinzon is pressrue from acceleration onto previos zone
348  * want this less than some fraction of total pressure */
349  /* >>chng 06 feb 01, change from init pres to current total pressure
350  * in const press high U ulirgs current pressure may be quite larger
351  * than init pressure due to continuum absorption */
352  double drConPres = 0.05*pressure.PresTotlCurr/
354  drChoice.insert( drConPres, "change in con accel" );
355  }
356  }
357 
358  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
359  /* gravitational pressure */
360  if( strcmp(dense.chDenseLaw,"CPRE") == 0 && (dark.lgNFW_Set || pressure.gravity_symmetry>=0) )
361  {
362  if( fabs( pressure.RhoGravity ) > 0. )
363  {
364  double drGravPres = 0.05 * pressure.PresTotlCurr * radius.drad / fabs( pressure.RhoGravity );
365  ASSERT( drGravPres > 0. );
366  drChoice.insert( drGravPres, "change in gravitational pressure" );
367  }
368  }
369 
370  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
371  /* check how temperature is changing */
372  double dTdStep = FLT_MAX;
373  if( nzone > 1 )
374  {
375  /* change in temperature; current= */
376  dTdStep = (phycon.te-OldTe)/phycon.te;
377  /* >>chng 02 dec 08, desired change in temperature per zone must not
378  * be smaller than allower error in temperature. For now use relative error
379  * in heating - - cooling balance. Better would be to also use c-h deriv wrt t
380  * to get slope */
381  /* >>chng 02 dec 11 rjrw change back to 0.03 -- improve dynamics.dRad criterion instead */
382  double x = 0.03;
383  /* >>chng 02 dec 07, do not do this if there is mild te jitter,
384  * so that dT is changing sign - this happens
385  * in ism.ini, very stable temperature with slight noise up and down */
386  if( dTdStep*OlddTdStep > 0. )
387  {
388  /* >>chng 96 may 30, variable depending on temp
389  * >>chng 96 may 31, allow 0.7 smaller, was 0.8
390  * >>chng 97 may 05, from 0.7 to 0.5 stop from punching through thermal front
391  * >>chng 04 mar 30, from 0.7 to 0.5 stop from punching through thermal front,
392  * for some reason factor was 0.7, not 0.5 as per previous change */
393  double absdt = fabs(dTdStep);
394  double drdTdStep = radius.drad*MAX2(0.5,x/absdt);
395  ostringstream oss;
396  oss << "change in temperature; current=" << scientific << setprecision(3) << phycon.te;
397  oss << fixed << " dT/T=" << dTdStep;
398  drChoice.insert( drdTdStep, oss.str() );
399  }
400  }
401  OlddTdStep = dTdStep;
402  OldTe = phycon.te;
403 
404  /* >>chng 02 oct 06, only check on opacity - interaction if not
405  * constant temperature - there were constant temperature models that
406  * extended to infinity but hung with last few photons and this test.
407  * better to ignore this test which is really for thermal models */
408  /* >>chng 06 mar 20, do not call if recombination logic in place */
410  {
411  double freqm = 0., opacm = 0.;
412  /* find freq where opacity largest and interaction rate is fastest */
413  ContRate(&freqm,&opacm);
414 
415  /* use optical depth at max interaction energy
416  * >>chng 96 jun 06 was drChange=0.15 changed to 0.3 for high Z models
417  * taking too many zones
418  * drInter = drChange / MAX(1e-30,opacm*FillFac) */
419 
420  double drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
421  ostringstream oss;
422  oss << "cont inter nu=" << scientific << setprecision(2) << freqm << " opac=" << opacm;
423  drChoice.insert( drInter, oss.str() );
424  }
425 
426  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
427  /* check whether change in wind velocity constrains DRAD
428  * WJH 22 May 2004: disable when we are near the sonic point since
429  * the velocity may be jumping all over the place but we just want
430  * to push through it as quickly as we can */
432  {
433  double v = fabs(wind.windv);
434  /* this is relative change in velocity */
435  double dVelRelative = fabs(wind.windv-OldWindVelocity)/
437 
438  const double WIND_CHNG_VELOCITY_RELATIVE = 0.01;
439  /*fprintf(ioQQQ,"DEBUG rad %.3f vel %.2e dVelRelative %.2e",
440  log10(radius.Radius) , wind.windv , dVelRelative );*/
441 
442  /* do not use this on first zone since do not have old velocity */
443  double winddr;
444  if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE && nzone>1 )
445  {
446  /* factor less than one if change larger than WIND_CHNG_VELOCITY_RELATIVE */
447  double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
448  /*fprintf(ioQQQ," DEBUG factor %.2e", factor );*/
449  factor = MAX2(0.8 , factor );
450  winddr = radius.drad * factor;
451  }
452  else
453  {
454  winddr = BigRadius;
455  }
456 
457  /* set dr from advective term */
459  {
460  winddr = MIN2( winddr , dynamics.dRad );
461  /*>>chng 04 oct 06, set dVelRelative to dynamics.dRad since dVelRelative is printed as part
462  * of reason for choosing this criteria, want to reflect valid reason. */
463  dVelRelative = dynamics.dRad;
464  }
465  ostringstream oss;
466  oss << "Wind, dVelRelative=" << scientific << setprecision(3);
467  oss << dVelRelative << " windv=" << wind.windv;
468  drChoice.insert( winddr, oss.str() );
469  }
470  /* remember old velocity */
471  OldWindVelocity = wind.windv;
472 
473  const double DNGLOB = 0.10;
474 
475  /* inside out globule */
476  if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
477  {
478  if( radius.glbdst < 0. )
479  {
480  fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
481  fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n",
482  radius.glbdst );
484  }
487  /* DNGLOB is relative change in density allowed this zone, 0.10 */
488  double GlobDr = ( fac2/factor > 1. + DNGLOB ) ? radius.drad*DNGLOB/(fac2/factor - 1.) : BigRadius;
489  GlobDr = MIN2(GlobDr,radius.glbdst/20.);
490  ostringstream oss;
491  oss << "GLOB law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
492  drChoice.insert( GlobDr, oss.str() );
493  }
494 
495  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
496  double hdnew = 0.;
497  if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
498  {
499  /* one of the special density laws, first get density at possible next dr */
500  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
501  {
503  radius.drad);
504  }
505  else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
506  {
508  radius.drad);
509  }
510  else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
511  {
513  }
514  else
515  {
516  fprintf( ioQQQ, " dlw insanity in radius_next\n" );
518  }
519  double drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
520  drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
521  ostringstream oss;
522  oss << "spec den law, new old den " << scientific << setprecision(2);
523  oss << hdnew << " " << dense.gas_phase[ipHYDROGEN];
524  drChoice.insert( drTab, oss.str() );
525  }
526 
527  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
528  /* special density law */
529  if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
530  {
531  double dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
533  /* DNGLOB is relative change in density allowed this zone, 0.10 */
534  double SpecDr;
535  if( dnew == 0. )
536  {
537  SpecDr = radius.drad*3.;
538  }
539  else if( dnew/DNGLOB > 1.0 )
540  {
541  SpecDr = radius.drad/(dnew/DNGLOB);
542  }
543  else
544  {
545  SpecDr = BigRadius;
546  }
547  ostringstream oss;
548  oss << "special law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
549  drChoice.insert( SpecDr, oss.str() );
550  }
551 
552  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
553  /* check grain line heating dominates
554  * this is important in PDR and HII region calculations
555  * >>chng 97 jul 03, added following check */
556  if( safe_div(thermal.heating(0,13),thermal.htot,0.0) > 0.2 )
557  {
558  double grfreqm = 0., gropacm = 0.;
559  /* >>chng 01 jan 03, following returns 0 when NO light at all,
560  * had failed with botched monitor */
561  GrainRateDr(&grfreqm,&gropacm);
562  double DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
563  ostringstream oss;
564  oss << "grain heating nu=" << scientific << setprecision(2) << grfreqm << " opac=" << gropacm;
565  drChoice.insert( DrGrainHeat, oss.str() );
566  }
567 
568  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
569  /* check if line heating dominates
570  * this is important in high metallicity models */
571  if( safe_div(thermal.heating(0,22),thermal.htot,0.0) > 0.2 )
572  {
573  long level = -1;
574  TransitionProxy t = FndLineHt(&level);
575 
576  if( t.Coll().heat()/thermal.htot > 0.1 )
577  {
578  ASSERT( t.associated() );
579 
580  double TauInwd = t.Emis().TauIn();
581  double DopplerWidth = t.Emis().dampXvel()/t.Emis().damp();
582  double TauDTau = t.Emis().PopOpc()*t.Emis().opacity()/DopplerWidth;
583 
584  fixit("all other line stacks need to be included here.");
585  // can we just sweep over line stack? Is that ready yet?
586 
587  /* in following logic cannot use usual inverse opacity,
588  * since line heating competes with escape probability,
589  * so is effective at surprising optical depths */
590  double drLineHeat = ( TauDTau > 0. ) ? MAX2(1.,TauInwd)*0.4/TauDTau : BigRadius;
591  ostringstream oss;
592  oss << "level " << level << " line heating, " << chLineLbl(t) << " TauIn " << scientific;
593  oss << setprecision(2) << TauInwd << " " << t.Emis().pump() << " " << t.Emis().Pesc();
594  drChoice.insert( drLineHeat, oss.str() );
595  }
596  }
597 
598  /* do not let change in cooling/heating due to H2 become too large */
599  if( lgFirstCall )
600  Old_H2_heat_cool = 0.;
601  else if( !thermal.lgTemperatureConstant )
602  {
603  /* this is case where temperature has not been set */
604  /* compare total heating - cooling due to h2 with total due to everything */
605  double H2_heat_cool = safe_div((fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)), thermal.htot,0.0);
606  // Leiden hack PDR sims can have H2 heating set by sims equation, does not change
607  double dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
608  if( H2_heat_cool > 0.1 && Old_H2_heat_cool>0. && dH2_heat_cool>SMALLFLOAT )
609  {
610  /* >>chng 05 oct 27, had been taking 20% of original radius - this caused zoning
611  * to become very fine and may not have been needed - change from 0.2 to 0.5 */
612  double drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
613  ostringstream oss;
614  oss << "change in H2 heating/cooling, d(c,h)/H " << scientific << setprecision(2);
615  oss << dH2_heat_cool;
616  drChoice.insert( drH2_heat_cool, oss.str() );
617  }
618  }
619  Old_H2_heat_cool = safe_div((fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) , thermal.htot, 0.0);
620 
621  /* >>chng 03 mar 04, add this logic */
622  /* do not let change in H2 and heavy element molecular abundances become too large
623  * "change in heav ele mole abundance, d(mole)/elem" */
624  if( nzone >= 4 )
625  {
626  /* first do H2 abundance */
627  /* >>chng 04 dec 15, do not use special logic when large h2 is turned on */
628  double Old_H2_abund = struc.H2_abund[nzone-3];
629  /* >>chng 03 jun 16, limit from 0.01 to 0.001, some models fell over H2 front due to
630  * large zoning, when large H2 was just this caused oscillations in solomon process */
631  /* >>chng 03 nov 22, from > 0.001 to > 3e-4, models that start almost in H2 have
632  * rapid increase in H2 at shallow depths, start sensing this sooner */
633  /* >>chng 03 dec 10, from 3e-4 to 1e-4, to get smaller chagned in leiden1.in test */
634  /* radius_next keys from change in H2 abundance, d(H2)/H */
635  /* >>chng 04 mar 11, start sensing H2 earlier since otherwise step size
636  * needs to become way too small tooo quickly. limit changed from 1e-4 to 1e-6 */
637  /* >>chng 04 jun 29, fromo > 1e-6 to >1e-8, some pdr's had too large a change in H2 */
638 
639  if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
640  {
641  double fac = 0.2;
642  /* this is percentage change in H2 density - "change in H2 abundance" */
643  double dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
644  /* in testing th85ism the dH2_abund did come out zero */
645  /* >>chng 03 jun 16, change d(H2) from 0.05 to 0.1, fine resolution of H2/H exposed
646  * small numerical oscillations in Solomon process */
647  /* >>chng 03 nov 22, smallest possible ratio of dr(next)/dr changed from
648  * 0.2 to 0.05, models that started almost in H2 front need to be able to sense it */
649  /* >>chng 04 mar 15, with such small possible changes in dr, only 0.05, a thermal front
650  * can easily cause large changes in H2 and T that are not due to zoning, but to the
651  * discontinuity. make smallest change larger to prevent hald due to dr */
652  /* >>chng 05 aug 23, thermal front allowed dr to become much too small
653  * chng from 0.02 to .6 */
654  dH2_abund = SDIV(dH2_abund);
655  double drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
656  ostringstream oss;
657  oss << "change in H2 abundance, d(H2)/H " << scientific << setprecision(2) << dH2_abund;
658  drChoice.insert( drH2_abund, oss.str() );
659  }
660  }
661 
662  int mole_dr_change = -1;
663  if( nzone >= 4 )
664  {
665  /* check how molecular fractions of all heavy elements are chaning relative
666  * to total gas phase abundance */
667  double max_change = 0.0;
668  /* >>chng 04 jun 02, upper limit had been all species but now limit to real
669  * molecules since do not want this logic to work with the ions */
670  for( int i=0; i < mole_global.num_calc; ++i )
671  {
672  realnum abund,
673  abund1,
674  abund_limit;
675  if( !mole_global.list[i]->isMonatomic() && mole_global.list[i]->isIsotopicTotalSpecies() )
676  {
677  /* >>chng 03 sep 21, add CO logic */
678  /* >>chng 04 mar 30, generalize to any molecule at all */
679  /* >>chng 04 mar 31 lower limit to abund had been 0.01, change
680  * to 0.001 to pick up approach to molecular fronts */
681  abund = 0.;
682  for( molecule::nNucsMap::iterator atom=mole_global.list[i]->nNuclide.begin();
683  atom != mole_global.list[i]->nNuclide.end(); ++atom)
684  {
685  long int nelem = atom->first->el->Z-1;
686  if (dense.lgElmtOn[nelem])
687  {
688  abund1 = (realnum)mole.species[i].den*atom->second/SDIV(dense.gas_phase[nelem]);
689  if (abund1 > abund)
690  {
691  abund = abund1;
692  }
693  }
694  }
695  /* is this an ice? need special logic for ice since density increases
696  * exponentially when grain temp falls below sublimation temperature
697  * >>chng 05 dec 06 - detect changes for smaller abundances for ices
698  * due to large changes in ice abundances */
699  if( mole_global.list[i]->lgGas_Phase )
700  {
701  abund_limit = 1e-3f;
702  }
703  else
704  {
705  /* this is an ice - track its abundance at lower values so that
706  * we resolve the sublimation transition region */
707  abund_limit = 1e-5f;
708  }
709 
710  if( abund > abund_limit )
711  {
712  /* the relative change in the abundance */
713  double relative_change =
714  2.0 * fabs( mole.species[i].den - struc.molecules[nzone-3][i] ) /
715  SDIV ( mole.species[i].den + struc.molecules[nzone-3][i] ) ;
716  if (relative_change > max_change)
717  {
718  max_change = relative_change;
719  mole_dr_change = i;
720  }
721  }
722  }
723  }
724  if( mole_dr_change >= 0 )
725  {
726  double dr_mole_abund = radius.drad * MAX2(0.6, 0.05/SDIV(max_change));
727  ostringstream oss;
728  oss << "change in abundance species=" << mole_global.list[mole_dr_change]->label << ", ";
729  oss << scientific << setprecision(2);
730  oss << "old=" << struc.molecules[nzone-3][mole_dr_change] << " new=" << mole.species[mole_dr_change].den;
731  drChoice.insert( dr_mole_abund, oss.str() );
732  }
733  }
734 
735  /* some consideration due to big H2 molecule */
736  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
737  drChoice.insert( (*diatom)->H2_DR(), "change in big H2 Solomon rate line opt depth" );
738 
739  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
740  /* can't make drmax large deep in model since causes feedback
741  * oscillations with changes in heating or destruction rates
742  * >>chng 96 oct 15, change from 5 to 10 */
743  /* >>chng 96 nov 11, had been 4 * drad up to 11, change to following
744  * to be similar to c90.01, max of 1.3 between 5 and 11
745  * >>chng 04 oct 29 geometry.DirectionalCosin was ioncorrect applied to this factor */
746  double drmax = ( nzone < 5 ) ? 4.*radius.drad : 1.3*radius.drad;
747  drChoice.insert( drmax, "DRMAX" );
748 
749  /* time dependent recombination - conditions become very homogeneous and
750  * crash into I fronts - use last iteration's radius as guess of current case */
751  double recom_dr_last_iter = BigRadius;
753  {
754  /* first time through nzone == 1 */
755  static long int nzone_recom = -1;
756  if( nzone <= 1 )
757  {
758  /* initialization case */
759  nzone_recom = 0;
760  }
762  nzone_recom < struc.nzonePreviousIteration )
763  {
764  ASSERT( nzone_recom>=0 && nzone_recom<struc.nzonePreviousIteration );
765  /* case where we are within previous computed structure
766  * first possibly increase nzone_recom so that it points deeper
767  * than this zone */
768  while( struc.depth_last[nzone_recom] < radius.depth &&
769  nzone_recom < struc.nzonePreviousIteration-1 )
770  ++nzone_recom;
771  recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
772  drChoice.insert( recom_dr_last_iter,
773  "previous iteration recom logic" );
774  }
775  }
776 
777  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
778  /* change in electron density - radius_next keys from change in elec den,
779  * remember old electron density during first call */
780  /* this is low ionization solution */
781  if( nzone > 2 )
782  {
783  /* next is-2 since nzone is on physics not c scale, and we want previous zone */
784  double Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[nzone-3][ipHYDROGEN];
785  double Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
786  double dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
787  double drEfrac;
788  if( dEfrac > SMALLFLOAT )
789  {
790  double fac = 0.04;
791  /* >>chng 03 dec 25, use smaller rel change in elec frac when most elec in ipMole or grains */
792  /* >>chng 04 sep 14, change to from from metals but comment out */
793  /* >>chng 04 sep 17, change to from from metals - uncomment */
794  if( dense.eden_from_metals > 0.5 )
795  {
796  fac = 0.04;
797  }
798  /* >>chng 04 feb 23, add test on hydrogen being predomintly
799  * recombined due to three-body recom, which is very sensitive
800  * to the electron density - but only do this in partially ionized medium */
801  else if( iso_sp[ipH_LIKE][ipHYDROGEN].RecomCollisFrac > 0.8 &&
804 
805  {
806  fac = 0.02;
807  }
808  /* >>chng 04 sep 17, change to 0.1 from 0.2 */
809  drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
810  }
811  else
812  {
813  drEfrac = BigRadius;
814  }
815  ostringstream oss;
816  oss << "change in elec den, rel chng: " << scientific << setprecision(3) << dEfrac;
817  oss << " cur=" << Efrac_new << " old=" << Efrac_old;
818  drChoice.insert( drEfrac, oss.str() );
819  }
820 
821  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
822  /* do not let thickness get too large */
823  if( nzone > 20 )
824  {
825  /* >>chng 02 nov 05, change from 1/20 to 1/10 wasted zones early on */
826  drChoice.insert( radius.depth/10., "relative depth" );
827  }
828 
829  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
830  /* case where stopping thickness or edge specified, need to approach slowly */
831  double thickness_total = BigRadius;
832  double DepthToGo = BigRadius;
833  if( StopCalc.HColStop < 5e29 )
834  {
835  double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
836  DepthToGo = MIN2(DepthToGo ,
837  (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
838  /* >>chng 03 oct 28, forgot to div col den above by eff density */
839  thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
840  }
841 
842  if( StopCalc.colpls < 5e29 )
843  {
844  double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
845  DepthToGo = MIN2(DepthToGo ,
846  (StopCalc.colpls-findspecieslocal("H+")->column) / coleff );
847  thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
848  }
849 
850  if( StopCalc.col_h2 < 5e29 )
851  {
852  /* >>chng 03 apr 15, add this molecular hydrogen */
853  double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
854  DepthToGo = MIN2(DepthToGo ,
855  (StopCalc.col_h2-findspecieslocal("H2")->column-findspecieslocal("H2*")->column) / coleff );
856  thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
857  }
858 
859  if( StopCalc.col_h2_nut < 5e29 )
860  {
861  /* neutral column density of H, 2 H_2 + H^0 */
862  double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
863  DepthToGo = MIN2(DepthToGo ,
865  thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
866  }
867 
868  if( StopCalc.col_H0_ov_Tspin < 5e29 )
869  {
870  /* >>chng 05 jan 09, add n(H0)/Tspin */
871  double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
872  DepthToGo = MIN2(DepthToGo ,
874  thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
875  }
876 
877  if( StopCalc.col_monoxco < 5e29 )
878  {
879  /* >>chng 03 apr 15, add this, CO */
880  double coleff = (double)SDIV( (findspecieslocal("CO")->den)*geometry.FillFac);
881  DepthToGo = MIN2(DepthToGo ,
882  (StopCalc.col_monoxco-findspecieslocal("CO")->column) / coleff );
883  thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
884  }
885 
886  if( StopCalc.col_species < 5e29 )
887  {
888  /* arbitrary species */
889  static molezone *SpeciesCurrent;
890  static bool lgMustFind=true;
891  if( lgMustFind )
892  {
893  lgMustFind = false;
894  SpeciesCurrent= findspecieslocal_validate(StopCalc.chSpeciesColumn.c_str());
895  }
896  double coleff = (double)SDIV( (SpeciesCurrent->den)*geometry.FillFac);
897  DepthToGo = MIN2(DepthToGo ,
898  (StopCalc.col_species-SpeciesCurrent->column) / coleff );
899  thickness_total = MIN2(thickness_total , StopCalc.col_species / coleff );
900  }
901 
902  if( StopCalc.colnut < 5e29 )
903  {
904  double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
905  DepthToGo = MIN2(DepthToGo ,
906  (StopCalc.colnut - findspecieslocal("H")->column) / coleff );
907  thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
908  }
909 
910  /* this is case where outer radius is set */
911  if( iterations.StopThickness[iteration-1] < 5e29 )
912  {
913  thickness_total = MIN2(thickness_total , iterations.StopThickness[iteration-1] );
914  DepthToGo = MIN2(DepthToGo ,
916  }
917 
918  /* this is case where stopping optical depth was specified */
920  {
921  /* end optical depth has been specified */
923  DepthToGo = MIN2(DepthToGo ,
925  }
926 
927  if( DepthToGo <= 0. )
928  TotalInsanity();
929 
930  /* set dr if one of above tests have triggered */
931  if( DepthToGo < BigRadius )
932  {
933  double depth_min = MIN2( DepthToGo , thickness_total/100. );
934  DepthToGo = MAX2( DepthToGo / 10. , depth_min );
935  /* want to
936 
937  approach the outer edge slowly,
938  * the need for this logic is most evident in brems.in -
939  * HI fraction varies across coronal model */
940  double drThickness = MIN2( thickness_total/10. , DepthToGo );
941  drChoice.insert( drThickness, "depth to go" );
942  }
943 
944  /* stop AV - usually this is dust, but we consider all opacity sources,
945  * so always include this */
946  /* compute some average grain properties */
947  double AV_to_go = BigRadius;
949  {
950  double SAFETY = 1. + 10.*DBL_EPSILON;
951  /* by default stop av is very large, and opacity can be very small, so ratio
952  * goes to inf - work with logs to see how big the number is */
953  /* apply safety margin to avoid updates to the total extinction getting lost
954  * in the numerical precision */
955  double ave = log10(SAFETY*StopCalc.AV_extended - rfield.extin_mag_V_extended) -
957  double avp = log10(SAFETY*StopCalc.AV_point - rfield.extin_mag_V_point) -
958  log10(rfield.opac_mag_V_point);
959  AV_to_go = MIN2( ave , avp );
960  if( AV_to_go < 37. )
961  {
962  AV_to_go = exp10( AV_to_go)/geometry.FillFac;
963  /* this is to make sure that we go slightly deeper than AV so that
964  * we trigger this stop */
965  AV_to_go *= 1.0001;
966  }
967  else
968  AV_to_go = BigRadius;
969  /*fprintf(ioQQQ,"DEBUG next dr %.3e %.3e %.3e\n", AV_to_go , ave, avp );*/
970  }
971  drChoice.insert( AV_to_go, "A_V to go" );
972 
973  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
974  /* spherical models, do not want delta R/R big */
975  double drSphere = radius.Radius*0.04;
976  drChoice.insert( drSphere, "sphericity" );
977 
978  /* optical depth to electron scattering */
979  double dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
980  /* allow larger dr when constant temperature,
981  * to prevent some ct models from taking too many cells */
983  dRTaue *= 3.;
984  drChoice.insert( dRTaue, "optical depth to electron scattering" );
985 
986  /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*/
987  if( dense.flong != 0. )
988  {
989  double drFluc = PI2/10./2./dense.flong;
990  /* >>chng 04 sep 18, caused cautions that ionization jumps occurred.
991  * set to have the value */
992  drFluc /= 2.;
993  drChoice.insert( drFluc, "density fluctuations" );
994  }
995 
996  /* keep dr constant in first two zones, if it wants to increase,
997  * but always allow it to decrease.
998  * to guard against large increases in efrac for balmer cont photo dominated models,
999  */
1000  double drStart = ( nzone <= 1 ) ? radius.drad : BigRadius;
1001  drChoice.insert( drStart, "capped to old DR in first zone" );
1002 
1003  // lgSdrmaxRel true if sdrmax is relative to current radius, false if limit in cm
1004  double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
1005  drChoice.insert( rfacmax*radius.sdrmax, "sdrmax" );
1006 
1007  /* >>chng 05 aug 05, in case of thermal front, where temp is falling quickly and
1008  * conditions change very fast, the zone thickness does not really affect the change
1009  * in conditions and can cause zoning to become very very thin, which causes
1010  * an abort. this occurs between 200 and 1000K. if we are doing temp soln,
1011  * temp is between these values, and temp is changing rapidly, do not make sone
1012  * thickness much smaller */
1013  /* do not use thermal front logic in case of recombination time dep cloud */
1014  if( nzone >= 5 && !dynamics.lgTimeDependentStatic )
1015  {
1016  /* >>chng 05 aug 23, upper bound of thermal from from 1000K to 4000K */
1017  /*if( phycon.te > 200. && phycon.te < 1000. && */
1018  if( phycon.te > 200. && phycon.te < 3000. &&
1019  /* >>chng 05 aug 23, from > 10% in zone to to two zones > 5%,
1020  * to fix leiden v3 with large H2 */
1021  (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
1022  (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
1023  (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
1024  {
1025  /* the 0.91 is to make dr unique, so that print statement that
1026  * follows will identify this reason */
1027  double drThermalFront = radius.drad * 0.91;
1028  drChoice.clear();
1029  drChoice.insert( drThermalFront, "thermal front logic" );
1030  }
1031  }
1032 
1033  ASSERT( drChoice.size() > 0 );
1034 
1035  /* dr = zero is a logical mistake */
1036  if( drChoice.begin()->first <= 0. )
1037  {
1038  drList::const_iterator ptr = drChoice.begin();
1039  fprintf( ioQQQ, " radius_next finds insane drNext: %.2e\n", ptr->first );
1040  fprintf( ioQQQ, " all drs follow:\n" );
1041  while( ptr != drChoice.end() )
1042  {
1043  fprintf( ioQQQ, " %.2e %s\n", ptr->first, ptr->second.c_str() );
1044  ++ptr;
1045  }
1047  }
1048 
1049  /* option to force min drad relative to depth */
1050  /* see google group thread "[cloudy-dev] sdrmin_rel_depth oscillation" ~2012/6/13 for discussion
1051  * of validity of this test */
1052  if( !radius.lgFixed && drChoice.begin()->first < radius.depth*radius.sdrmin_rel_depth )
1053  {
1054  drChoice.clear();
1055  drChoice.insert( radius.depth*radius.sdrmin_rel_depth, "sdrmin_rel_depth" );
1056  }
1057 
1058  /* option to force min drad */
1059  double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
1060  if( drChoice.begin()->first < rfacmin*radius.sdrmin )
1061  {
1062  drChoice.clear();
1063  drChoice.insert( rfacmin*radius.sdrmin, "sdrmin" );
1064  }
1065 
1066  /* this is general code that prevents zone thickness drNext from
1067  * becoming too thin, something that can happen for various bad reasons
1068  * HOWEVER we do not want to do this test for certain density laws,
1069  * for which very small zone thicknesses are unavoidable
1070  * the special cases are:
1071  * special density law,
1072  * globule density law,
1073  * carbon +-0 i front
1074  * the fluctuations command
1075  * drMinimum was set in radius_first to either sdrmin (set drmin) or
1076  * some fraction of the initial radius - it is always set
1077  * to something non-trivial.
1078  * sdrmin is set with the "set dr" command - is SMALLFLOAT by default */
1079  if( strcmp(dense.chDenseLaw,"DLW1") != 0 &&
1080  strcmp(dense.chDenseLaw ,"GLOB") != 0 &&
1081  dense.flong == 0. &&
1082  // stopping on depth to go is not a fault
1083  drChoice.begin()->first != DepthToGo &&
1084  // stopping on Av to go is not a fault
1085  drChoice.begin()->first != AV_to_go )
1086  {
1087  fixit( "drMinimum does not get reevaluated in each iteration,"
1088  " so time-dependent solutions just use the first drad and Hden."
1089  " Re-eval here for now to keep this minimum from blowing." );
1091 
1092  /* don't let dr get smaller than drMinimum, if this resets drNext
1093  * then code stops with warning that zones got too thin
1094  * this can happen due to numerical oscillations, which the code
1095  * tries to damp out by making the zones thinner.
1096  * scale by radius.Radius/radius.rinner to stop very spherical
1097  * simulations from false trigger on dr too small */
1098  /* drMinimum is drad * hden, to make it proportional to optical depth.
1099  * This avoids false trigger across thermal fronts. */
1101  {
1103  /* leaving this at true will cause the model to stop with a warning
1104  * that the zone thickness is too small */
1105  radius.lgDrMinUsed = true;
1106  /* set abort handler */
1107  lgAbort = true;
1108  /* must decrement nzone, since we will not complete this zone, and will not have
1109  * valid structure data for it */
1110  --nzone;
1111  fprintf( ioQQQ,
1112  "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
1113  "This is zone %ld iteration %ld. The thickness was %.2e\n",
1114  nzone,
1115  iteration,
1116  radius.drNext);
1117  fprintf( ioQQQ,
1118  " If this simulation seems reasonable you can override this limit "
1119  "with the command SET DRMIN %.2e\n\n",
1120  radius.drNext*2);
1121  /* set abort flag */
1122  lgAbort = true;
1123  return 1;
1124  }
1125  }
1126 
1127  /* factor to allow for slop in floating numbers */
1128  const double Z = 1.0001;
1129 
1130  /* following is to make thickness of model exact
1131  * n.b., on last zone, drNext can be NEGATIVE!!
1132  * DEPTH was incremented at start of zone calc in newrad,
1133  * has been outer edge of zone all throughout */
1134  double drOuterRadius = (iterations.StopThickness[iteration-1]-radius.depth)*Z;
1135  drChoice.insert( drOuterRadius, "outer radius" );
1136 
1137  if( cosmology.lgDo )
1138  {
1139  drChoice.clear();
1140  double drHz = SPEEDLIGHT / GetHubbleFactor( cosmology.redshift_current );
1141  drChoice.insert( drHz, "c over H(z)" );
1142  }
1143 
1144  // choose the smallest dR as the next choice
1145  radius.drNext = drChoice.begin()->first;
1146 
1147  /* set choice flag if provided -- in practice, just for rt.lgMaserSetDR */
1148  drChoice.begin()->second.select( );
1149 
1150  /* all this is to only save on last iteration
1151  * the save dr command is not really a save command, making this necessary
1152  * lgDRon is set true if "save dr" entered */
1153  /* lgDRPLst was set true if "save" had "last" on it */
1154  bool lgDoPun = ( save.lgDROn && !( save.lgDRPLst && !iterations.lgLastIt ) );
1155  if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
1156  {
1157  fprintf( save.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t%s\n", nzone, radius.depth,
1158  radius.drNext, radius.Depth2Go, drChoice.begin()->second.c_str() );
1159  }
1160 
1161  ASSERT( radius.drNext > 0. );
1162 
1163  if( trace.lgTrace )
1164  {
1165  fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was %.4e\n",
1166  radius.drNext, radius.drad );
1167  }
1168 
1169  return 0;
1170 }
1171 
1172 /*ContRate called by radius_next to find energy of maximum continuum-gas interaction */
1173 STATIC void ContRate(double *freqm,
1174  double *opacm)
1175 {
1176  long int i,
1177  ipHi,
1178  iplow,
1179  limit;
1180  double FreqH,
1181  Freq_nonIonizing,
1182  Opac_Hion,
1183  Opac_nonIonizing,
1184  Rate_max_Hion,
1185  Rate_max_nonIonizing;
1186 
1187  DEBUG_ENTRY( "ContRate()" );
1188 
1189  /*
1190  * find maximum continuum interaction rate,
1191  * these should be reset in following logic without exception,
1192  * if they are still zero at the end we have a logical error
1193  */
1194  Rate_max_nonIonizing = 0.;
1195  Freq_nonIonizing = 0.;
1196  Opac_nonIonizing = 0.;
1197 
1198  /* this must be reset to val >= 0 */
1199  *opacm = -1.;
1200  *freqm = -1.;
1201 
1202  /* do up to carbon photo edge if carbon is turned on */
1203  /* >>>chng 00 apr 07, add test for whether element is turned on */
1204  if( dense.lgElmtOn[ipCARBON] )
1205  {
1206  /* carbon is turned on, use carbon 1 edge */
1207  ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
1208  }
1209  else
1210  {
1211  /* carbon turned off, use hydrogen Balmer continuum */
1212  ipHi = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon-1;
1213  }
1214 
1215  /* start at plasma frequency */
1216  for( i=rfield.ipPlasma; i < ipHi; i++ )
1217  {
1218  /* this does not have grain opacity since grains totally passive
1219  * at energies smaller than CI edge */
1220  if( rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*(opac.opacity_abs[i] -
1221  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1222  {
1223  Rate_max_nonIonizing = rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*
1225  Freq_nonIonizing = rfield.anu(i);
1226  Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1227  }
1228  }
1229 
1230  /* not every continuum extends beyond C edge-this whole loop can add to zero
1231  * use total opacity here
1232  * test is to put in fir continuum if free free heating is important */
1233  if( safe_div(CoolHeavy.brems_heat_total,thermal.htot,0.0) > 0.05 )
1234  {
1235  /* this is index for energy where cloud free free optical depth is unity,
1236  * is zero if no freq are opt thin */
1237  iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1238  }
1239  else
1240  {
1241  /* >>>chng 00 apr 07, from Heavy.ipHeavy[0][5] to ipHi defined above, since
1242  * would crash if element not defined */
1243  iplow = ipHi;
1244  }
1245  /* do not go below plasma frequency */
1246  iplow = MAX2( rfield.ipPlasma , iplow );
1247 
1248  /* this energy range from carbon edge to hydrogen edge */
1249  limit = MIN2(rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1);
1250  for( i=iplow; i < limit; i++ )
1251  {
1252  if( rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*(opac.opacity_abs[i] -
1253  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
1254  {
1255  Rate_max_nonIonizing = rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*
1257  Freq_nonIonizing = rfield.anu(i);
1258  Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1259  }
1260  }
1261 
1262  /* variables to check continuum interactions over Lyman continuum */
1263  Rate_max_Hion = 0.;
1264  Opac_Hion = 0.;
1265  FreqH = 0.;
1266 
1267  /* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
1268  for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1; i < rfield.nflux; i++ )
1269  {
1270  if( rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*(opac.opacity_abs[i] -
1271  gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
1272  {
1273  /* Rate_max_Hion = anu(i)*flux(i)/widflx(i)*opac(i) */
1274  Rate_max_Hion = rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*
1276  FreqH = rfield.anu(i);
1277  Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
1278  }
1279  }
1280 
1281  /* use Lyman continuum if its opacity is larger than non-h ion */
1282  if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
1283  {
1284  /* this happens for laser source - use Lyman continuum */
1285  *opacm = Opac_Hion;
1286  *freqm = FreqH;
1287  }
1288  /* >>chng 05 aug 03, add last test on Opac_Hion for case where we go very
1289  * deep and very little radiation is left */
1290  else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/SDIV(Rate_max_nonIonizing) > 1e-10 && Opac_Hion > SMALLFLOAT )
1291  {
1292  /* use Lyman continuum */
1293  *opacm = Opac_Hion;
1294  *freqm = FreqH;
1295  }
1296  else
1297  {
1298  /* not much rate in the Lyman continuum, stick with low energy */
1299  *opacm = Opac_nonIonizing;
1300  *freqm = Freq_nonIonizing;
1301  }
1302 
1303 # if 0
1304  /*>>chng 06 aug 12, do not let zones become optically thick to
1305  * peak of dust emission - one of NA's vastly optically thick dust clouds
1306  * did not conserve total energy - very opticallly thick to ir dust emission
1307  * so ir was absorbed and reemitted many times
1308  * this makes sure the cells remain optically thin to dust emission */
1309  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
1310  {
1311  double fluxGrainPeak = -1.;
1312  long int ipGrainPeak = -1;
1313  long int ipGrainPeak2 = -1;
1314 
1315  i = 0;
1316  while( rfield.anu(i) < 0.0912 && i < (rfield.nflux-2) )
1317  {
1318  /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1319  * check on optical depth for wavelengths greater than 1 micron */
1320  if( gv.GrainEmission[i]*rfield.anu(i)*opac.opacity_abs[i] > fluxGrainPeak )
1321  {
1322  ipGrainPeak = i;
1323  fluxGrainPeak = gv.GrainEmission[i]*rfield.anu(i)*opac.opacity_abs[i];
1324  }
1325  ++i;
1326  }
1327  ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
1328 
1329  /* >>chng 06 aug 23. We have just found the wavelength and flux at the peak of the IR emission.
1330  * Now we want to find the wavelength, short of the peak, which corresponds to 5% of the
1331  * peak emission. This wavelength will be where the code checks to make sure the zone has
1332  * not become optically thick. Since dust opacity generally decreases with wavelength, this assures that
1333  * the optical depths are small for all wavelengths where the flux is appreciable. Tests show that
1334  * this allows for flux/luminosity conservation to within ~5% for an AV of 1e4 mag and at least 2 iterations*/
1335  i = ipGrainPeak;
1336  /* >>chng 06 aug 23. Only want to do this test for the IR dust continuum, therefore only
1337  * check on opacities for wavelengths shortward of the peak and greater than 1 micron
1338  * this routine can be called with the dust emission totally zero - it is only evaluated
1339  * late to save time. don't do the following tests if peak dust emission is zero */
1340  if( fluxGrainPeak > 0. )
1341  {
1342  while( rfield.anu(i) < 0.0912 && i < (rfield.nflux-2) )
1343  {
1344  /* find wavelength where flux = 5% of peak flux, shortward of the peak */
1345  if( gv.GrainEmission[i]*rfield.anu(i)*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
1346  {
1347  /* This will be the array number and flux at 10% of the peak */
1348  ipGrainPeak2 = i;
1349  }
1350  ++i;
1351  }
1352  ASSERT( ipGrainPeak2>=0 );
1353 
1354  /* use this as a limit to the dr if opacity is larger */
1355  if( opac.opacity_abs[ipGrainPeak2] > *opacm )
1356  {
1357  /* scale opacity by factor of 10, which further decreases the zone thickness*/
1358  *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
1359  *freqm = rfield.anu(ipGrainPeak2);
1360  /*fprintf(ioQQQ,"DEBUT opac peak %.2e %.2e \n",
1361  *opacm , *freqm );*/
1362  }
1363  }
1364  }
1365 # endif
1366 
1367  {
1368  /* following should be set true to print contributors */
1369  enum {DEBUG_LOC=false};
1370  if( DEBUG_LOC )
1371  {
1372  fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
1373  Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
1374  Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
1375  );
1376 
1377  }
1378  }
1379 
1380  /* these were set to -1 at start, and must have been reset in one of the
1381  * two loops. Logic error if still <0. */
1382  /* >>chng 05 aug 03, change logic to -1 on entry and check at least zero
1383  * here - will be zero if NO radiation field exists, perhaps since totally
1384  * attenuated */
1385  ASSERT( *opacm >= 0. && *freqm >= 0. );
1386  return;
1387 }
1388 
1389 /*GrainRateDr called by radius_next to find grain heating rate dr */
1390 STATIC void GrainRateDr(double *grfreqm,
1391  double *gropacm)
1392 {
1393  long int i,
1394  iplow;
1395  double xMax;
1396 
1397  DEBUG_ENTRY( "GrainRateDr()" );
1398 
1399  /* in all following changed from anu2 to anu july 25 95
1400  *
1401  * find maximum continuum interaction rate */
1402 
1403  /* not every continuum extends beyond C edge-this whole loop can add to zero
1404  * use total opacity here
1405  * test is to put in fir continuum if free free heating is important */
1406  if( safe_div(CoolHeavy.brems_heat_total,thermal.htot,0.0) > 0.05 )
1407  {
1408  /* this is pointer to energy where cloud free free optical depth is unity,
1409  * is zero if no freq are opt thin */
1410  iplow = MAX2(1 , rfield.ipEnergyBremsThin );
1411  }
1412  else
1413  {
1414  /* do up to carbon photo edge if carbon is turned on */
1415  /* >>>chng 00 apr 07, add test for whether element is turned on */
1416  if( dense.lgElmtOn[ipCARBON] )
1417  {
1418  /* carbon is turned on, use carbon 1 edge */
1419  iplow = Heavy.ipHeavy[ipCARBON][0];
1420  }
1421  else
1422  {
1423  /* carbon truned off, use hydrogen balmer continuum */
1424  iplow = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon;
1425  }
1426  }
1427 
1428  xMax = -1.;
1429  /* integrate up to H edge */
1430  for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
1431  {
1432  if( rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*opac.opacity_abs[i] > xMax )
1433  {
1434  xMax = rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*
1435  opac.opacity_abs[i];
1436  *grfreqm = rfield.anu(i);
1437  *gropacm = opac.opacity_abs[i];
1438  }
1439  }
1440  /* integrate up to heii edge if he is turned on,
1441  * this logic will not make sense if grains on but he off, which in itself makes no sense*/
1442  if( dense.lgElmtOn[ipHELIUM] )
1443  {
1444  for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
1445  {
1446  if( rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*opac.opacity_abs[i] > xMax )
1447  {
1448  xMax = rfield.anu(i)*rfield.flux[0][i]/rfield.widflx(i)*
1449  opac.opacity_abs[i];
1450  *grfreqm = rfield.anu(i);
1451  *gropacm = opac.opacity_abs[i];
1452  }
1453  }
1454  }
1455 
1456  /* possible that there is NO ionizing radiation, in extreme cases,
1457  * if so then xMax will still be negative */
1458  if( xMax <= 0. )
1459  {
1460  *gropacm = 0.;
1461  *grfreqm = 0.;
1462  }
1463  return;
1464 }
realnum col_h2
Definition: stopcalc.h:74
multimap< double, drChoiceItem >::const_iterator const_iterator
Definition: radius_next.cpp:76
bool lgContRadPresOn
Definition: pressure.h:65
t_mole_global mole_global
Definition: mole.cpp:7
double Radius
Definition: radius.h:31
realnum colnut
Definition: stopcalc.h:69
double depth
Definition: radius.h:31
vector< double > dstab
Definition: grainvar.h:527
double htot
Definition: thermal.h:169
realnum drMinimum
Definition: radius.h:178
t_thermal thermal
Definition: thermal.cpp:6
long int iptnu
Definition: stopcalc.h:29
realnum & opacity() const
Definition: emission.h:638
t_colden colden
Definition: colden.cpp:5
double * opacity_abs
Definition: opacity.h:103
double exp10(double x)
Definition: cddefines.h:1383
realnum dr_ionfrac_limit
Definition: struc.h:84
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
realnum * depth_last
Definition: struc.h:25
long int ipEnergyBremsThin
Definition: rfield.h:228
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double widflx(size_t i) const
Definition: mesh.h:147
t_opac opac
Definition: opacity.cpp:5
int num_calc
Definition: mole.h:362
realnum ** flux
Definition: rfield.h:70
t_struc struc
Definition: struc.cpp:6
t_Heavy Heavy
Definition: heavy.cpp:5
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
const realnum SMALLFLOAT
Definition: cpu.h:246
multimap< double, drChoiceItem > m_map
Definition: radius_next.cpp:74
double opac_mag_V_point
Definition: rfield.h:267
realnum AV_extended
Definition: stopcalc.h:89
long int mas_species
Definition: rt.h:208
const char * c_str() const
Definition: radius_next.cpp:66
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
realnum * ednstr
Definition: struc.h:25
realnum drChange
Definition: radius.h:188
bool lgDrBug
Definition: trace.h:61
realnum ** molecules
Definition: struc.h:71
t_phycon phycon
Definition: phycon.cpp:6
realnum col_monoxco
Definition: stopcalc.h:86
size_t size() const
Definition: radius_next.cpp:95
double sdrmax
Definition: radius.h:158
realnum column
Definition: mole.h:422
void clear()
Definition: radius_next.cpp:83
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
STATIC void ContRate(double *freqm, double *opacm)
bool lgTemperatureConstant
Definition: thermal.h:44
vector< double > StopThickness
Definition: iterations.h:77
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
realnum FillFac
Definition: geometry.h:29
long int nzone
Definition: cddefines.cpp:14
double HeatH2Dexc_used
Definition: hmi.h:140
Definition: mole.h:378
t_dynamics dynamics
Definition: dynamics.cpp:42
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:807
bool lgDo
Definition: cosmology.h:44
double PresTotlCurr
Definition: pressure.h:46
double anu(size_t i) const
Definition: mesh.h:111
realnum colpls
Definition: stopcalc.h:69
realnum glbpow
Definition: radius.h:133
bool lgDrMinUsed
Definition: radius.h:185
molezone * findspecieslocal_validate(const char buf[])
t_dense dense
Definition: global.cpp:15
bool associated() const
Definition: transition.h:54
double sdrmin_rel_depth
Definition: radius.h:161
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nflux_with_check
Definition: rfield.h:51
realnum glbdst
Definition: radius.h:133
bool lgMaserSetDR
Definition: rt.h:201
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int iteration
Definition: cddefines.cpp:16
t_trace trace
Definition: trace.cpp:5
double drad
Definition: radius.h:31
double rinner
Definition: radius.h:31
double sdrmin
Definition: radius.h:157
t_abund abund
Definition: abund.cpp:5
realnum pinzon
Definition: pressure.h:69
t_geometry geometry
Definition: geometry.cpp:5
t_dark_matter dark
Definition: dark_matter.cpp:5
realnum HColStop
Definition: stopcalc.h:69
const_iterator end() const
Definition: radius_next.cpp:91
realnum redshift_current
Definition: cosmology.h:26
double sound_speed_isothermal
Definition: timesc.h:55
double RhoGravity
Definition: pressure.h:82
bool lgLeidenHack
Definition: mole.h:334
vector< realnum > GrainEmission
Definition: grainvar.h:580
bool lgStrongDLimbo
Definition: pressure.h:135
bool lgSdrminRel
Definition: radius.h:165
const int ipH1s
Definition: iso.h:29
double & heat() const
Definition: collision.h:224
double HeatH2Dish_used
Definition: hmi.h:140
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
realnum & dampXvel() const
Definition: emission.h:598
long int ipPlasma
Definition: rfield.h:434
EmissionList::reference Emis() const
Definition: transition.h:447
double Depth2Go
Definition: radius.h:31
t_mole_local mole
Definition: mole.cpp:8
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
const_iterator begin() const
Definition: radius_next.cpp:87
const TransitionProxy FndLineHt(long int *level)
float realnum
Definition: cddefines.h:124
valarray< class molezone > species
Definition: mole.h:468
#define EXIT_FAILURE
Definition: cddefines.h:168
realnum AccelTotalOutward
Definition: wind.h:52
double H0_ov_Tspin
Definition: colden.h:52
realnum col_species
Definition: stopcalc.h:134
realnum & Pesc() const
Definition: emission.h:568
realnum AV_point
Definition: stopcalc.h:89
bool lgSdrmaxRel
Definition: radius.h:166
vector< diatomics * > diatoms
Definition: h2.cpp:8
string chSpeciesColumn
Definition: stopcalc.h:133
realnum flong
Definition: dense.h:262
#define cdEXIT(FAIL)
Definition: cddefines.h:484
realnum DirectionalCosin
Definition: geometry.h:25
DepthTable DLW
Definition: dense.h:198
bool lgDRPLst
Definition: save.h:447
double tabval(double r0, double depth) const
Definition: depth_table.cpp:8
double brems_heat_total
Definition: coolheavy.h:36
realnum col_H0_ov_Tspin
Definition: stopcalc.h:83
void insert(double val, const string &str, bool *flag=NULL)
Definition: radius_next.cpp:77
bool lgGrainPhysicsOn
Definition: grainvar.h:481
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
t_iterations iterations
Definition: iterations.cpp:6
int gravity_symmetry
Definition: pressure.h:84
double column(const genericState &gs)
realnum * H2_abund
Definition: struc.h:71
t_radius radius
Definition: radius.cpp:5
double heating(long nelem, long ion)
Definition: thermal.h:186
t_timesc timesc
Definition: timesc.cpp:7
double dRad
Definition: dynamics.h:144
bool lgElmtOn[LIMELM]
Definition: dense.h:160
double extin_mag_V_point
Definition: rfield.h:260
realnum gas_phase[LIMELM]
Definition: dense.h:76
realnum glbrad
Definition: radius.h:133
long int mas_lo
Definition: rt.h:208
STATIC void GrainRateDr(double *grfreqm, double *gropacm)
static const double SAFETY
#define ASSERT(exp)
Definition: cddefines.h:617
bool lgLastIt
Definition: iterations.h:47
double Tspin21cm
Definition: hyperfine.h:56
char chDenseLaw[5]
Definition: dense.h:176
const int ipH2s
Definition: iso.h:30
double extin_mag_V_extended
Definition: rfield.h:264
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
double opac_mag_V_extended
Definition: rfield.h:267
double den
Definition: mole.h:421
CollisionProxy Coll() const
Definition: transition.h:463
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double & PopOpc() const
Definition: emission.h:658
t_cosmology cosmology
Definition: cosmology.cpp:8
long int mas_hi
Definition: rt.h:208
const int ipHELIUM
Definition: cddefines.h:349
realnum xMassDensity
Definition: dense.h:101
realnum * drad_last
Definition: struc.h:25
realnum ** gas_phase
Definition: struc.h:75
double H2_total
Definition: hmi.h:25
int radius_next(void)
bool lgDROn
Definition: save.h:447
double eden
Definition: dense.h:201
realnum ** TauAbsGeo
Definition: opacity.h:90
realnum tauend
Definition: stopcalc.h:23
bool lgSonicPoint
Definition: pressure.h:128
#define MAX2(a, b)
Definition: cddefines.h:828
realnum & damp() const
Definition: emission.h:608
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
MoleculeList list
Definition: mole.h:365
double dense_fabden(double radius, double depth)
bool lgRecom
Definition: dynamics.h:108
realnum GetHubbleFactor(realnum z)
Definition: cosmology.cpp:10
bool lgDenFlucOn
Definition: dense.h:255
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
double pow(double x, int i)
Definition: cddefines.h:782
long int n_initial_relax
Definition: dynamics.h:132
const int ipCARBON
Definition: cddefines.h:353
realnum dTauMase
Definition: rt.h:198
double dense_parametric_wind(double rad)
bool lgStatic(void) const
Definition: wind.h:24
drChoiceItem(const string &str, bool *flag)
Definition: radius_next.cpp:60
realnum * testr
Definition: struc.h:25
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
#define fixit(a)
Definition: cddefines.h:416
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
double eden_from_metals
Definition: dense.h:235
realnum glbden
Definition: radius.h:133
double lgFixed
Definition: radius.h:159
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
long int nzonePreviousIteration
Definition: struc.h:22
FILE * ipDRout
Definition: save.h:446
bool lgDustOn() const
Definition: grainvar.h:477
const int ipHYDROGEN
Definition: cddefines.h:348
long int nflux
Definition: rfield.h:48
realnum colden[NCOLD]
Definition: colden.h:32
realnum *** xIonDense
Definition: struc.h:64
long int mas_ion
Definition: rt.h:208
realnum col_h2_nut
Definition: stopcalc.h:80
void select() const
Definition: radius_next.cpp:61
realnum windv
Definition: wind.h:18
realnum & TauIn() const
Definition: emission.h:458
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
bool lgAbort
Definition: cddefines.cpp:10
double drNext
Definition: radius.h:66
t_rt rt
Definition: rt.cpp:5
double & pump() const
Definition: emission.h:518