cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
zone_startend.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 /*ZoneEnd last routine called after all zone calculations, before iter_end_check,
4  * upon exit radiation field is for outer edge of current zone */
5 /*ZoneStart set variables that change with each zone, like radius, depth,
6  * upon exit flux will be flux at center of zone about to be computed */
7 #include "cddefines.h"
8 #include "zones.h"
9 #include "phycon.h"
10 #include "opacity.h"
11 #include "rfield.h"
12 #include "struc.h"
13 #include "thermal.h"
14 #include "dense.h"
15 #include "h2.h"
16 #include "geometry.h"
17 #include "iterations.h"
18 #include "conv.h"
19 #include "save.h"
20 #include "dynamics.h"
21 #include "radius.h"
22 #include "doppvel.h"
23 #include "mole.h"
24 #include "cosmology.h"
25 /* this is number of zones to include in guess for next temperatures */
26 #define IOFF 3
27 
28 void ZoneStart(const char *chMode)
29 {
30  bool lgNeOscil,
31  lgTeOscil;
32  long int i;
33 
34  double
35  /* this is correction factor for dilution of radiation
36  * across current zone, set in ZoneStart to correct
37  * flux for center of current zone, and undone in ZoneDone */
38  DilutionCorrec ,
39  drFac ,
40  dTauThisZone,
41  outwrd,
42  ratio,
43  rm,
44  rad_middle_zone,
45  vin,
46  vout;
47 
48  static double DTaver , DEaver,
49  dt2 , de2;
50 
51  DEBUG_ENTRY( "ZoneStart()" );
52 
54  /* this is a turbulent velocity power law. */
55  if( DoppVel.lgTurbLawOn )
56  {
59  }
60 
61  /* this sub can be called in two ways, 'incr' means increment
62  * radius variables, while 'init' means only initialize rest */
63  /* called at start of a zone to update all variables
64  * having to do with starting this zone */
65 
66  /* first establish current filling factor (normally 1) since used within
67  * following branches */
71 
72  if( strcmp(chMode,"init") == 0 )
73  {
74  /* initialize the variables at start of calculation, after this exits
75  * radius is equal to the initial radius, not outer edge of zone */
76 
77  /* depth to illuminated face */
78  radius.depth = 1e-30;
79 
80  /* integral of depth times filling factor */
82  /* effective radius */
84 
85  /* reset radius to inner radius, at this point equal to illuminated face radius */
88 
89  /* thickness of next zone */
91 
92  /* depth to illuminated face */
94 
95  /* depth variable for globule */
97 
98  /* this is case where outer radius is set */
99  if( iterations.StopThickness[iteration-1] < 5e29 )
100  {
102  }
103  else if( iteration > 1 )
104  {
105  /* this case second or later iteration, use previous thickness */
107  }
108  else
109  {
110  /* first iteration and depth not set, so we cannot estimate it */
111  radius.Depth2Go = 0.;
112  }
113  }
114  else if( strcmp(chMode,"incr") == 0 )
115  {
116  /* update radius variables - called by cloudy at start of this zone's calcs */
119  /* >>chng 06 mar 21, remember largest and smallest dr over this iteration
120  * for use in recombination time dep logic */
123 
124  ASSERT( radius.drad > 0. );
125  radius.depth += radius.drad;
127 
128  /* effective radius */
130 
131  /* integral of depth times filling factor */
133 
134  /* decrement Depth2Go but do not let become negative */
136 
137  /* increment the radius, during the calculation Radius is the
138  * outer radius of the current zone*/
140 
141  /* Radius is now outer edge of this zone since incremented above,
142  * so need to add drad to it */
144 
145  /***********************************************************************
146  *
147  * attenuate rfield to center of this zone
148  *
149  ***********************************************************************/
150 
151  /* radius was just incremented above, so this resets continuum to
152  * flux at center of zone we are about to compute. radius will be
153  * incremented immediately following this. this correction is undone
154  * when ZoneDone called */
155 
156  /* this will be the optical thickness of the next zone
157  * AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command,
158  * option for illumination of slab at an angle */
159 
160  /* radius.dRNeff is next dr with filling factor, this will only be
161  * used to get approximate correction for attenuation
162  * of radiation in this zone,
163  * equations derived in hazy, Netzer&Ferland 83, for factor accounting
164  * any changes here should be checked with both sphonly.in and pphonly*/
165  /* honlyotspp seems most sensitive to the 1.35 */
167 
168  if( cosmology.lgDo )
169  drFac = 0.;
170 
171  /* dilute radiation so that rfield is in center of zone, this
172  * goes into tmn at start of zone here, but is later divided out
173  * when ZoneEnd called */
174  DilutionCorrec = 1./POW2(
176 
177  /* note that this for loop is through <= nflux, to include the [nflux]
178  * unit integration verification token. The token was set in
179  * in MadeDiffuse and propagated out in metdif. the opacity at that energy is
180  * zero, so only the geometrical part of tmn will affect things. The final
181  * sum is verified in PrtComment */
182  for( i=0; i <= rfield.nflux; i++ )
183  {
184  dTauThisZone = opac.opacity_abs[i]*drFac;
185  /* TMN is array of scale factors which account for attenuation
186  * of continuum across the zone (necessary to conserve energy
187  * at the 1% - 5% level.) sphere effects in
188  * drNext was set by NEXTDR and will be next dr */
189 
190  if( dTauThisZone < 1e-4 )
191  {
192  /* this small tau limit is the most likely branch, about 60% in parispn.in */
193  opac.tmn[i] = 1.f;
194  }
195  else if( dTauThisZone < 5. )
196  {
197  /* this middle tau limit happens in the remaining 40% */
198  opac.tmn[i] = (realnum)((1. - exp(-dTauThisZone))/(dTauThisZone));
199  }
200  else
201  {
202  /* this happens almost not at all */
203  opac.tmn[i] = (realnum)(1./(dTauThisZone));
204  }
205 
206  /* now add on correction for dilution across this zone */
207  opac.tmn[i] *= (realnum)DilutionCorrec;
208 
209  rfield.flux_beam_const[i] *= opac.tmn[i];
210  rfield.flux_beam_time[i] *= opac.tmn[i];
211  rfield.flux_isotropic[i] *= opac.tmn[i];
214 
215  /* >>chng 03 nov 08, update SummedCon here since flux changed */
216  rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
218  }
219 
220  /* following is distance to globule, usually does not matter */
222 
223  /* if gt 5th zone, and not constant pressure, and not oscillating te
224  * then guess next temperature : option to predict next temperature
225  * NZLIM is dim of structure variables saving temp, do data if nzone>NZLIM
226  * IOFF is number of zones to look over, set set to 3 in the define above */
227  /* >>chng 01 mar 12, add option to not do this, set with no tepred command */
228  if( (strcmp(dense.chDenseLaw,"CPRE") != 0) &&
230  (nzone > IOFF + 1) )
231  {
234  lgTeOscil = false;
235  lgNeOscil = false;
236  dt2 = 0.;
237  de2 = 0.;
238  DTaver = 0.;
239  DEaver = 0.;
240  for( i=nzone - IOFF-1; i < (nzone - 1); i++ )
241  {
242  double dt1 = dt2;
243  double de1 = de2;
244  /* this will get the average change in temperature for the
245  * past IOFF zones */
246  dt2 = struc.testr[i-1] - struc.testr[i];
247  de2 = struc.ednstr[i-1] - struc.ednstr[i];
248  DTaver += dt2;
249  DEaver += de2;
250  if( dt1*dt2 < 0. )
251  {
252  lgTeOscil = true;
253  }
254  if( de1*de2 < 0. )
255  {
256  lgNeOscil = true;
257  }
258  }
259 
260  /* option to guess next electron temperature if no oscillations */
261  if( !lgTeOscil )
262  {
263  DTaver /= (double)(IOFF);
264  /* don't want to over correct, use smaller of two */
265  dt2 = fabs(dt2);
266  rm = fabs(DTaver);
267  DTaver = sign(MIN2(rm,dt2),DTaver);
268  /* do not let it change by more than 5% of current temp */
269  /* >>chng 03 mar 18, from 5% to 1% - convergence is much better
270  * now, and throwing the next Te off by 5% could easily disturb
271  * the solvers - primal.in was a good case */
272  DTaver = sign(MIN2(rm,0.01*phycon.te),DTaver);
273  /* this actually changes the temperature */
274  TempChange(phycon.te - DTaver , true);
275  }
276  else
277  {
278  /* temp was oscillating - too dangerous to do anything */
279  DTaver = 0.;
280  }
281  /* this is used to remember the proposed temperature, for output
282  * in the save predictor command */
284 
285  /* option to guess next electron density if no oscillations */
286  if( !lgNeOscil )
287  {
288  DEaver /= IOFF;
289  de2 = fabs(de2);
290  rm = fabs(DEaver);
291  /* don't want to over correct, use smaller of two */
292  DEaver = sign(MIN2(rm,de2),DEaver);
293  /* do not let it change by more than 5% of current temp */
294  DEaver = sign(MIN2(rm,0.05*dense.eden),DEaver);
295  /* this actually changes the temperature */
296  EdenChange( dense.eden - DEaver );
297  }
298  else
299  {
300  /* temp was oscillating - too dangerous to do anything */
301  DEaver = 0.;
302  }
303  /* this is used to remember the proposed temperature, for output
304  * in the save predictor command */
306  /* must call TempChange since ionization has changed, there are some
307  * terms that affect collision rates (H0 term in electron collision) */
308  TempChange(phycon.te , false);
309  }
310  }
311 
312  else
313  {
314  fprintf( ioQQQ, " PROBLEM ZoneStart called with insane argument, %4.4s\n",
315  chMode );
316  /* TotalInsanity exits after announcing a problem */
317  TotalInsanity();
318  }
319 
320  /* do advection if enabled */
322  DynaStartZone();
323 
324  /* clear flag indicating the ionization convergence failures
325  * have ever occurred in current zone
326  conv.lgConvIonizThisZone = false; */
328 
329  /* this will say how many times the large H2 molecule has been called in this zone -
330  * if not called (due to low H2 abundance) then not need to update its line arrays */
331  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom)
332  (*diatom)->nCall_this_zone = 0;
333 
334  /* check whether zone thickness is too small relative to radius */
335  if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
336  {
337  ratio = radius.drad/(radius.glbdst + radius.drad);
338  /* this only matters for globule command */
339  if( ratio < 1e-14 )
340  {
341  radius.lgdR2Small = true;
342  }
343  else
344  {
345  radius.lgdR2Small = false;
346  }
347  }
348 
349  /* area factor, ratio of radius to outer edge of this zone
350  * relative to radius of illuminated face of cloud
351  * is unity in intensity case, per unit area */
353  if( radius.lgPredLumin )
354  {
355  if( save.lgLuminosityOld )
356  {
357  /* use old style, save continuum luminosity per unit cloud area
358  * at inner radius, when SET SAVE LUMINOSITY OLD used */
360  radius.PI4_rinner_sq = 1.;
361  }
362  else
363  {
366  }
367  }
368  else
369  {
370  radius.PI4_Radius_sq = 1.;
371  radius.PI4_rinner_sq = 1.;
372  }
373 
374  /* following only approximate, used for center of next zone */
376 
377  /* this is the middle of the zone */
378  rad_middle_zone = radius.Radius - radius.dRadSign*radius.drad/2.;
379 
380  /* this is used for long slit integrals */
381  if( radius.drad/radius.Radius > 0.01 )
382  {
383  double Ropp = radius.Radius - radius.dRadSign*radius.drad;
385  }
386  else
387  radius.darea_x_fillfac = PI2*rad_middle_zone*radius.drad*geometry.FillFac;
388 
389  /* Radius is outer radius of this zone, so radius - drad is inner radius
390  * rinner is inner radius of nebula; dVeffVol is the effective vol element
391  * dVeffVol is vol rel to inner radius, so has units cm
392  * for plane parallel dVeffVol is dReff */
393  if( radius.drad/radius.Radius > 0.01 )
394  {
395  double r1 = radius.Radius - radius.dRadSign*radius.drad;
396  double rin_zone = min(r1,radius.Radius);
397  double rout_zone = max(r1,radius.Radius);
398  vin = pow2(rin_zone/radius.rinner)*rin_zone/3.;
399  if( rin_zone > radius.CylindHigh )
400  {
401  // the volume of the cap of a sphere is given here:
402  // http://en.wikipedia.org/wiki/Spherical_cap
403  // we need two of these...
404  double h = rin_zone-radius.CylindHigh;
405  double v2cap = pow2(h/radius.rinner)*(rin_zone - h/3.)/2.;
406  vin -= v2cap;
407  }
408  vout = pow2(rout_zone/radius.rinner)*rout_zone/3.;
409  if( rout_zone > radius.CylindHigh )
410  {
411  double h = rout_zone-radius.CylindHigh;
412  double v2cap = pow2(h/radius.rinner)*(rout_zone - h/3.)/2.;
413  vout -= v2cap;
414  }
415  /* this is the usual case the full volume, just the difference in the two vols */
416  /* this will become the true vol when multiplied by 4pi*rinner^2 */
417  radius.dVeffVol = (vout - vin)*geometry.FillFac;
418  /* this is correction for slit projected onto resolved object -
419  * this only happens when aperture command is entered.
420  * when slit and cylinder are both used it is assumed that the slit
421  * is oriented along the rotational axis of the cylinder
422  * the case where the slit is perpendicular to this axis is identical
423  * to the normal spherical case, so can be done by removing the cylinder command
424  * slits oriented at any other angle can be done by dividing the cylinder height
425  * by the cosine of the angle */
426  /* default of iEmissPower is 2, set to 0 is just aperture beam,
427  * and to 1 if aperture long slit set */
428  if( geometry.iEmissPower == 2 )
429  {
431  }
432  else if( geometry.iEmissPower == 1 )
433  {
434  double ain = (rin_zone/radius.rinner)*rin_zone/2.;
435  if( rin_zone > radius.CylindHigh )
436  {
437  // the area of a circular segment is given here:
438  // http://en.wikipedia.org/wiki/Circular_segment
439  // we need two of those...
440  double Theta = 2.*acos(min(radius.CylindHigh/rin_zone,1.));
441  ain *= 1. - max(Theta - sin(Theta),0.)/PI;
442  }
443  double aout = (rout_zone/radius.rinner)*rout_zone/2.;
444  if( rout_zone > radius.CylindHigh )
445  {
446  double Theta = 2.*acos(min(radius.CylindHigh/rout_zone,1.));
447  aout *= 1. - max(Theta - sin(Theta),0.)/PI;
448  }
449  radius.dVeffAper = (aout - ain)*geometry.FillFac;
450  }
451  else if( geometry.iEmissPower == 0 )
452  {
454  }
455  }
456  else
457  {
458  /* thin cell limit */
459  /* rad_middle_zone is the middle of the zone */
460  radius.dVeffVol = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac*
461  (min(rad_middle_zone,radius.CylindHigh)/radius.rinner);
462  if( geometry.iEmissPower == 2 )
463  {
465  }
466  else if( geometry.iEmissPower == 1 )
467  {
468  radius.dVeffAper = (rad_middle_zone/radius.rinner)*radius.drad*geometry.FillFac;
469  if( rad_middle_zone > radius.CylindHigh )
470  {
471  double Theta = 2.*acos(min(radius.CylindHigh/rad_middle_zone,1.));
472  double q = sqrt(max(1.-pow2(radius.CylindHigh/rad_middle_zone),0.))*rad_middle_zone;
473  radius.dVeffAper *= 1. - max(Theta - sin(Theta),0.)/PI - max(1. - cos(Theta),0.)*radius.CylindHigh/(PI*q);
474  }
475  }
476  else if( geometry.iEmissPower == 0 )
477  {
479  }
480  }
481 
482  /* covering factor, default is unity */
485 
486  /* these are needed for line intensities in outward beam
487  * lgSphere set, geometry.covrt usually 1, 0 when not lgSphere
488  * so outward is 1 when lgSphere set 1/2 when not set, */
489  outwrd = (1. + geometry.covrt)/2.;
490 
491  /*>>>chng 99 apr 23, from above to below */
492  /*radius.dVolOutwrd = outwrd*POW2( (radius.Radius-radius.drad_x_fillfac/2.)/radius.Radius) *
493  radius.drad;*/
494  /* this includes covering fact, the effective vol,, and 1/r^2 dilution,
495  * dReff includes filling factor. the covering factor part is 1 for sphere,
496  * 1/2 for open */
497  /*radius.dVolOutwrd = outwrd*radius.Radius*radius.drad_x_fillfac/(radius.Radius +
498  2.*radius.drad);*/
499  /* dReff from above, includes filling factor */
501  ASSERT( radius.dVolOutwrd > 0. );
502 
503  /* following will be used to "uncorrect" for this in lines when predicting continua
504  radius.GeoDil = radius.Radius/(radius.Radius + 2.*radius.drad);*/
505 
506  /* this should multiply the line to become the net inward. geo part is 1/2 for
507  * open geometry, 0 for closed. for case of isotropic emitter only this and dVolOutwrd
508  * above are used */
510 
511  if( geometry.lgSphere )
512  {
513  /* inward beams do not go in when lgSphere set since geometry symmetric */
514  radius.BeamInIn = 0.;
516  2.*radius.drad);
517  }
518  else
519  {
521 
522  /* inward beams do not go out */
523  radius.BeamInOut = 0.;
524  }
525  /* factor for outwardly directed part of line */
527  2.*radius.drad);
528  return;
529 }
530 
531 void ZoneEnd(void)
532 {
533  DEBUG_ENTRY( "ZoneEnd()" );
534 
535  /***********************************************************************
536  *
537  * correct rfield for attenuation from center of zone to inner edge
538  *
539  ***********************************************************************/
540 
541  // verify sanity of rfield.nPositive, ==0 is case with no light
542  //fprintf(ioQQQ,"DEBUGGG, %li %li %.3e\n",rfield.nPositive, rfield.nflux,
543  // rfield.SummedCon[rfield.nPositive-1]) ;
544  if( rfield.nPositive>0 )
545  {
547  if (rfield.FluxFaint > 0.0)
548  {
549  for( long i=rfield.nPositive; i<rfield.nflux; ++i )
550  ASSERT( rfield.SummedCon[i]==0.);
551  }
552  }
553  //fprintf(ioQQQ,"DEBUGGG, %li %li \n",rfield.nflux,rfield.nflux_with_check) ;
555  /* radius is outer radius of this zone, this resets continuum to
556  * flux at illuminated face of zone we have already computed. */
557 
558  /* opac.tmn defined in ZoneStart */
559  /* undilute radiation so that rfield is at face of zone */
560  /* NB, upper limit of sum includes [nflux] to test unit verification cell */
561  for( long i=0; i <= rfield.nflux; i++ )
562  {
563  rfield.flux_beam_const[i] /= opac.tmn[i];
564  rfield.flux_beam_time[i] /= opac.tmn[i];
565  rfield.flux_isotropic[i] /= opac.tmn[i];
568  /* >>chng 03 nov 08, update SummedCon here since flux changed */
569  rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
571  }
572 
574 
575  /* do advection if enabled */
577  DynaEndZone();
578 
579  if (0)
581 
582  return;
583 }
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
t_thermal thermal
Definition: thermal.cpp:6
realnum * flux_isotropic
Definition: rfield.h:73
bool lgTurbLawOn
Definition: doppvel.h:33
double * opacity_abs
Definition: opacity.h:103
double drad_mid_zone
Definition: radius.h:31
void setTrimming()
Definition: rfield.cpp:92
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
double PI4_rinner_sq
Definition: radius.h:31
t_opac opac
Definition: opacity.cpp:5
realnum ** flux
Definition: rfield.h:70
t_struc struc
Definition: struc.cpp:6
double dRNeff
Definition: radius.h:95
void ZoneEnd(void)
double CylindHigh
Definition: radius.h:125
void resetCountersZone()
Definition: conv.h:322
long int iEmissPower
Definition: geometry.h:71
bool lgAdvection
Definition: dynamics.h:66
bool lgTimeDependentStatic
Definition: dynamics.h:102
void DynaStartZone(void)
Definition: dynamics.cpp:375
t_conv conv
Definition: conv.cpp:5
realnum * ednstr
Definition: struc.h:25
T sign(T x, T y)
Definition: cddefines.h:846
realnum fiscal
Definition: geometry.h:29
t_phycon phycon
Definition: phycon.cpp:6
double * SummedCon
Definition: rfield.h:163
realnum TurbVel
Definition: doppvel.h:21
realnum * SummedOcc
Definition: rfield.h:165
vector< double > StopThickness
Definition: iterations.h:77
realnum covgeo
Definition: geometry.h:45
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum FillFac
Definition: geometry.h:29
long int nzone
Definition: cddefines.cpp:14
t_DoppVel DoppVel
Definition: doppvel.cpp:5
t_dynamics dynamics
Definition: dynamics.cpp:42
double TeProp
Definition: phycon.h:99
#define MIN2(a, b)
Definition: cddefines.h:807
bool lgDo
Definition: cosmology.h:44
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
realnum TurbVelZero
Definition: doppvel.h:25
void ZoneStart(const char *chMode)
double EdenInit
Definition: phycon.h:99
t_dense dense
Definition: global.cpp:15
bool lgPredNextTe
Definition: thermal.h:40
double dVolReflec
Definition: radius.h:103
long int nflux_with_check
Definition: rfield.h:51
double depth_x_fillfac
Definition: radius.h:79
void DynaEndZone(void)
Definition: dynamics.cpp:833
double dVeffVol
Definition: radius.h:86
realnum glbdst
Definition: radius.h:133
double PI4_Radius_sq
Definition: radius.h:31
bool lgSphere
Definition: geometry.h:34
long int iteration
Definition: cddefines.cpp:16
double drad
Definition: radius.h:31
double rinner
Definition: radius.h:31
t_geometry geometry
Definition: geometry.cpp:5
#define POW2
Definition: cddefines.h:983
double Depth2Go
Definition: radius.h:31
t_rfield rfield
Definition: rfield.cpp:9
realnum TurbVelLaw
Definition: doppvel.h:29
realnum * convoc
Definition: rfield.h:115
float realnum
Definition: cddefines.h:124
double BeamInIn
Definition: radius.h:107
bool lgLuminosityOld
Definition: save.h:405
void mole_rk_bigchange(void)
bool lgdR2Small
Definition: radius.h:117
vector< diatomics * > diatoms
Definition: h2.cpp:8
double dr_max_last_iter
Definition: radius.h:182
long max(int a, long b)
Definition: cddefines.h:821
double FluxFaint
Definition: rfield.h:58
realnum DirectionalCosin
Definition: geometry.h:25
realnum covrt
Definition: geometry.h:61
long min(int a, long b)
Definition: cddefines.h:766
double depth_mid_zone
Definition: radius.h:31
t_iterations iterations
Definition: iterations.cpp:6
#define IOFF
t_radius radius
Definition: radius.cpp:5
double dRadSign
Definition: radius.h:73
bool lgPredLumin
Definition: radius.h:144
realnum glbrad
Definition: radius.h:133
double dVolOutwrd
Definition: radius.h:102
double Radius_mid_zone
Definition: radius.h:31
double BeamInOut
Definition: radius.h:110
#define ASSERT(exp)
Definition: cddefines.h:617
realnum covaper
Definition: geometry.h:54
char chDenseLaw[5]
Definition: dense.h:176
T pow2(T a)
Definition: cddefines.h:985
double drad_x_fillfac
Definition: radius.h:76
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
t_cosmology cosmology
Definition: cosmology.cpp:8
double BeamOutOut
Definition: radius.h:113
double TeInit
Definition: phycon.h:99
double EdenProp
Definition: phycon.h:99
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
realnum * SummedDif
Definition: rfield.h:164
double pow(double x, int i)
Definition: cddefines.h:782
double darea_x_fillfac
Definition: radius.h:82
realnum * testr
Definition: struc.h:25
double dr_min_last_iter
Definition: radius.h:181
void EdenChange(double EdenNew)
Definition: eden_change.cpp:11
realnum * tmn
Definition: opacity.h:148
realnum * flux_beam_time
Definition: rfield.h:76
double r1r0sq
Definition: radius.h:31
realnum * flux_beam_const
Definition: rfield.h:76
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
long int nflux
Definition: rfield.h:48
long int nPositive
Definition: rfield.h:55
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
double dVeffAper
Definition: radius.h:92
realnum filpow
Definition: geometry.h:29
double drNext
Definition: radius.h:66