cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
zone_startend.cpp
Go to the documentation of this file.
1
/* This file is part of Cloudy and is copyright (C)1978-2008 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 "
phycon.h
"
9
#include "
opacity.h
"
10
#include "
rfield.h
"
11
#include "
struc.h
"
12
#include "
thermal.h
"
13
#include "
dense.h
"
14
#include "
h2.h
"
15
#include "
geometry.h
"
16
#include "
conv.h
"
17
#include "
dynamics.h
"
18
#include "
radius.h
"
19
#include "
zones.h
"
20
#include "
doppvel.h
"
21
/* this is number of zones to include in guess for next temperatures */
22
#define IOFF 3
23
24
void
ZoneStart
(
const
char
*chMode)
25
{
26
bool
lgNeOscil,
27
lgTeOscil;
28
long
int
i;
29
30
double
dt1 , de1,
31
/* this is correction factor for dilution of radiation
32
* across current zone, set in ZoneStart to correct
33
* flux for center of current zone, and undone in ZoneDone */
34
DilutionCorrec ,
35
drFac ,
36
dTauThisZone,
37
outwrd,
38
ratio,
39
rm,
40
rad_middle_zone,
41
vin,
42
vout;
43
44
static
double
DTaver , DEaver,
45
dt2 , de2;
46
47
DEBUG_ENTRY
(
"ZoneStart()"
);
48
50
/* this is a turbulent velocity power law. */
51
if
(
DoppVel
.
lgTurbLawOn
)
52
{
53
DoppVel
.
TurbVel
=
DoppVel
.
TurbVelZero
*
54
pow( (
realnum
)(
radius
.
Radius
/
radius
.
rinner
),
DoppVel
.
TurbVelLaw
);
55
}
56
57
/* this sub can be called in two ways, 'incr' means increment
58
* radius variables, while 'init' means only initialize rest */
59
/* called at start of a zone to update all variables
60
* having to do with starting this zone */
61
62
/* first establish current filling factor (normally 1) since used within
63
* following branches */
64
geometry
.
FillFac
= (
realnum
)(
geometry
.
fiscal
*
65
pow(
radius
.
Radius
/
radius
.
rinner
,(
double
)
geometry
.
filpow
));
66
geometry
.
FillFac
= (
realnum
)
MIN2
(1.,
geometry
.
FillFac
);
67
68
if
( strcmp(chMode,
"init"
) == 0 )
69
{
70
/* initialize the variables at start of calculation, after this exits
71
* radius is equal to the initial radius, not outer edge of zone */
72
73
/* depth to illuminated face */
74
radius
.
depth
= 1e-30;
75
76
/* integral of depth times filling factor */
77
radius
.
depth_x_fillfac
=
radius
.
depth
*
geometry
.
FillFac
;
78
/* effective radius */
79
radius
.
drad_x_fillfac
=
radius
.
drad
*
geometry
.
FillFac
;
80
81
/* reset radius to inner radius, at this point equal to illuminated face radius */
82
radius
.
Radius
=
radius
.
rinner
;
83
radius
.
Radius_mid_zone
=
radius
.
rinner
+
radius
.
drad
/2.;
84
85
/* thickness of next zone */
86
radius
.
drNext
=
radius
.
drad
;
87
88
/* depth to illuminated face */
89
radius
.
depth_mid_zone
=
radius
.
drad
/2.;
90
91
/* depth variable for globule */
92
radius
.
glbdst
=
radius
.
glbrad
;
93
94
/* this is case where outer radius is set */
95
if
(
radius
.
router
[
iteration
-1] < 5e29 )
96
{
97
radius
.
Depth2Go
=
radius
.
router
[
iteration
-1];
98
}
99
else
if
(
iteration
> 1 )
100
{
101
/* this case second or later iteration, use previous thickness */
102
radius
.
Depth2Go
=
radius
.
router
[
iteration
-2];
103
}
104
else
105
{
106
/* first iteration and depth not set, so we cannot estimate it */
107
radius
.
Depth2Go
= 0.;
108
}
109
}
110
else
if
( strcmp(chMode,
"incr"
) == 0 )
111
{
112
/* update radius variables - called by cloudy at start of this zone's calcs */
113
radius
.
drad
=
radius
.
drNext
;
114
/* >>chng 06 mar 21, remember largest and smallest dr over this iteration
115
* for use in recombination time dep logic */
116
radius
.
dr_min_last_iter
=
MIN2
(
radius
.
dr_min_last_iter
,
radius
.
drNext
);
117
radius
.
dr_max_last_iter
=
MAX2
(
radius
.
dr_max_last_iter
,
radius
.
drNext
);
118
119
ASSERT
(
radius
.
drad
> 0. );
120
radius
.
depth
+=
radius
.
drad
;
121
radius
.
depth_mid_zone
=
radius
.
depth
-
radius
.
drad
/2.;
122
123
/* effective radius */
124
radius
.
drad_x_fillfac
=
radius
.
drad
*
geometry
.
FillFac
;
125
126
/* integral of depth times filling factor */
127
radius
.
depth_x_fillfac
+=
radius
.
drad_x_fillfac
;
128
129
/* decrement Depth2Go but do not let become negative */
130
radius
.
Depth2Go
=
MAX2
( 0.,
radius
.
Depth2Go
-
radius
.
drad
);
131
132
/* increment the radius, during the calculation Radius is the
133
* outer radius of the current zone*/
134
radius
.
Radius
+=
radius
.
drad
*
radius
.
dRadSign
;
135
136
/* Radius is now outer edge of this zone since incremented above,
137
* so need to add drad to it */
138
radius
.
Radius_mid_zone
=
radius
.
Radius
-
radius
.
drad
/2.;
139
140
/***********************************************************************
141
*
142
* attenuate rfield to center of this zone
143
*
144
***********************************************************************/
145
146
/* radius was just incremented above, so this resets continuum to
147
* flux at center of zone we are about to compute. radius will be
148
* incremented immediately following this. this correction is undone
149
* when ZoneDone called */
150
151
/* this will be the optical thickness of the next zone
152
* AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command,
153
* option for illumination of slab at an angle */
154
155
/* radius.dRNeff is next dr with filling factor, this will only be
156
* used to get approximate correction for attenuation
157
* of radiation in this zone,
158
* equations derived in hazy, Netzer&Ferland 83, for factor accounting
159
* any changes here should be checked with both sphonly.in and pphonly*/
160
/* honlyotspp seems most sensitive to the 1.35 */
161
drFac =
radius
.
drad
*
geometry
.
FillFac
*
geometry
.
DirectionalCosin
*1.35;
162
163
/* dilute radiation so that rfield is in center of zone, this
164
* goes into tmn at start of zone here, but is later divided out
165
* when ZoneEnd called */
166
DilutionCorrec = 1./
POW2
(
167
(
radius
.
Radius
-
radius
.
drad
/2.)/(
radius
.
Radius
-
radius
.
drad
) );
168
169
/* note that this for loop is through <= nflux, to include the [nflux]
170
* unit integration verification token. The token was set in
171
* in MadeDiffuse and propagated out in metdif. the opacity at that energy is
172
* zero, so only the geometrical part of tmn will affect things. The final
173
* sum is verified in PrtComment */
174
for
( i=0; i <=
rfield
.
nflux
; i++ )
175
{
176
dTauThisZone =
opac
.
opacity_abs
[i]*drFac;
177
/* TMN is array of scale factors which account for attenuation
178
* of continuum across the zone (necessary to conserve energy
179
* at the 1% - 5% level.) sphere effects in
180
* drNext was set by NEXTDR and will be next dr */
181
182
if
( dTauThisZone < 1e-4 )
183
{
184
/* this small tau limit is the most likely branch, about 60% in parispn.in */
185
opac
.
tmn
[i] = 1.f;
186
}
187
else
if
( dTauThisZone < 5. )
188
{
189
/* this middle tau limit happens in the remaining 40% */
190
opac
.
tmn
[i] = (
realnum
)((1. - exp(-dTauThisZone))/(dTauThisZone));
191
}
192
else
193
{
194
/* this happens almost not at all */
195
opac
.
tmn
[i] = (
realnum
)(1./(dTauThisZone));
196
}
197
198
/* now add on correction for dilution across this zone */
199
opac
.
tmn
[i] *= (
realnum
)DilutionCorrec;
200
201
rfield
.
flux_beam_const
[i] *=
opac
.
tmn
[i];
202
rfield
.
flux_beam_time
[i] *=
opac
.
tmn
[i];
203
rfield
.
flux_isotropic
[i] *=
opac
.
tmn
[i];
204
rfield
.
flux
[i] =
rfield
.
flux_beam_const
[i] +
rfield
.
flux_beam_time
[i] +
205
rfield
.
flux_isotropic
[i];
206
207
/* actually do the corrections
208
rfield.flux[i] *= opac.tmn[i]; */
209
210
/* >>chng 03 nov 08, update SummedCon here since flux changed */
211
rfield
.
SummedCon
[i] =
rfield
.
flux
[i] +
rfield
.
SummedDif
[i];
212
}
213
214
/* following is distance to globule, usually does not matter */
215
radius
.
glbdst
-= (
realnum
)
radius
.
drad
;
216
217
/* if gt 5th zone, and not constant pressure, and not oscillating te
218
* then guess next temperature : option to predict next temperature
219
* NZLIM is dim of structure variables saving temp, do data if nzone>NZLIM
220
* IOFF is number of zones to look over, set set to 3 in the define above */
221
/* >>chng 01 mar 12, add option to not do this, set with no tepred command */
222
if
( (strcmp(
dense
.
chDenseLaw
,
"CPRE"
) != 0) &&
223
thermal
.
lgPredNextTe
&&
224
(
nzone
>
IOFF
+ 1) )
225
{
226
phycon
.
TeInit
=
phycon
.
te
;
227
phycon
.
EdenInit
=
dense
.
eden
;
228
lgTeOscil =
false
;
229
lgNeOscil =
false
;
230
dt1 = 0.;
231
dt2 = 0.;
232
de1 = 0.;
233
de2 = 0.;
234
DTaver = 0.;
235
DEaver = 0.;
236
for
( i=
nzone
-
IOFF
-1; i < (
nzone
- 1); i++ )
237
{
238
dt1 = dt2;
239
de1 = de2;
240
/* this will get the average change in temperature for the
241
* past IOFF zones */
242
dt2 =
struc
.
testr
[i-1] -
struc
.
testr
[i];
243
de2 =
struc
.
ednstr
[i-1] -
struc
.
ednstr
[i];
244
DTaver += dt2;
245
DEaver += de2;
246
if
( dt1*dt2 < 0. )
247
{
248
lgTeOscil =
true
;
249
}
250
if
( de1*de2 < 0. )
251
{
252
lgNeOscil =
true
;
253
}
254
}
255
256
/* option to guess next electron temperature if no oscillations */
257
if
( !lgTeOscil )
258
{
259
DTaver /= (double)(
IOFF
);
260
/* don't want to over correct, use smaller of two */
261
dt2 = fabs(dt2);
262
rm = fabs(DTaver);
263
DTaver =
sign
(
MIN2
(rm,dt2),DTaver);
264
/* do not let it change by more than 5% of current temp */
265
/* >>chng 03 mar 18, from 5% to 1% - convergence is much better
266
* now, and throwing the next Te off by 5% could easily disturb
267
* the solvers - primal.in was a good case */
268
DTaver =
sign
(
MIN2
(rm,0.01*
phycon
.
te
),DTaver);
269
/* this actually changes the temperature */
270
TempChange
(
phycon
.
te
- DTaver ,
true
);
271
}
272
else
273
{
274
/* temp was oscillating - too dangerous to do anything */
275
DTaver = 0.;
276
}
277
/* this is used to remember the proposed temperature, for output
278
* in the punch predictor command */
279
phycon
.
TeProp
=
phycon
.
te
;
280
281
/* option to guess next electron density if no oscillations */
282
if
( !lgNeOscil )
283
{
284
DEaver /=
IOFF
;
285
de2 = fabs(de2);
286
rm = fabs(DEaver);
287
/* don't want to over correct, use smaller of two */
288
DEaver =
sign
(
MIN2
(rm,de2),DEaver);
289
/* do not let it change by more than 5% of current temp */
290
DEaver =
sign
(
MIN2
(rm,0.05*
dense
.
eden
),DEaver);
291
/* this actually changes the temperature */
292
dense
.
eden
-= DEaver;
293
}
294
else
295
{
296
/* temp was oscillating - too dangerous to do anything */
297
DEaver = 0.;
298
}
299
/* this is used to remember the proposed temperature, for output
300
* in the punch predictor command */
301
phycon
.
EdenProp
=
dense
.
eden
;
302
/* must call TempChange since ionization has changed, there are some
303
* terms that affect collision rates (H0 term in electron collision) */
304
TempChange
(
phycon
.
te
,
false
);
305
}
306
}
307
308
else
309
{
310
fprintf(
ioQQQ
,
" PROBLEM ZoneStart called with insane argument, %4.4s\n"
,
311
chMode );
312
/* TotalInsanity exits after announcing a problem */
313
TotalInsanity
();
314
}
315
316
/* find the mean dr for the past few zones - this var is used
317
* in setting the optical scale in induced processes, and large
318
* changes in it will drive oscillations in the H+ fraction */
319
radius
.
drad_x_fillfac_mean
=
radius
.
drad_x_fillfac
;
320
i = 1;
321
while
( i < 5 && nzone-i-1>=0 )
322
{
323
radius
.
drad_x_fillfac_mean
+=
struc
.
drad_x_fillfac
[
nzone
-i-1];
324
++i;
325
}
326
radius
.
drad_x_fillfac_mean
/= i;
327
328
/* do advection if enabled */
329
if
(
dynamics
.
lgAdvection
)
330
DynaStartZone
();
331
332
/* clear flag indicating the ionization convergence failures
333
* have ever occurred in current zone
334
conv.lgConvIonizThisZone = false; */
335
336
/* this will say how many times the large H2 molecule has been called in this zone -
337
* if not called (due to low H2 abundance) then not need to update its line arrays */
338
h2
.
nCallH2_this_zone
= 0;
339
340
/* check whether zone thickness is too small relative to radius */
341
if
( strcmp(
dense
.
chDenseLaw
,
"GLOB"
) == 0 )
342
{
343
ratio =
radius
.
drad
/(
radius
.
glbdst
+
radius
.
drad
);
344
/* this only matters for globule command */
345
if
( ratio < 1e-14 )
346
{
347
radius
.
lgdR2Small
=
true
;
348
}
349
else
350
{
351
radius
.
lgdR2Small
=
false
;
352
}
353
}
354
355
/* area factor, ratio of radius to out edge of this zone
356
* relative to radius of illuminated face of cloud */
357
/*radius.r1r0sq = (realnum)POW2(
358
(radius.Radius - radius.drad*radius.dRadSign/2.)/radius.rinner);*/
359
/*>>>chng 99 jan 25, definition of r1r0sq is now outer edge of current zone
360
* relative to inner edge of cloud */
361
radius
.
r1r0sq
=
POW2
(
radius
.
Radius
/
radius
.
rinner
);
362
363
/* following only approximate, used for center of next zone */
364
radius
.
dRNeff
=
radius
.
drNext
*
geometry
.
FillFac
;
365
366
/* this is the middle of the zone */
367
rad_middle_zone =
radius
.
Radius
-
radius
.
drad
/2.;
368
369
/* Radius is outer radius of this zone, so radius - drad is inner radius
370
* rinner is inner radius of nebula; dVeff is the effective vol element
371
* dVeff is vol rel to inner radius, so has units cm
372
* for plane parallel dVeff is dReff */
373
if
(
radius
.
drad
/
radius
.
Radius
> 0.01 )
374
{
375
vin = ((
radius
.
Radius
-
radius
.
drad
)/
radius
.
rinner
)*
376
(
MIN2
((
radius
.
Radius
-
radius
.
drad
),
radius
.
CylindHigh
)/
radius
.
rinner
)*
377
(
radius
.
Radius
-
radius
.
drad
);
378
379
vout = (
radius
.
Radius
/
radius
.
rinner
)*
380
(
MIN2
(
radius
.
Radius
,
radius
.
CylindHigh
)/
radius
.
rinner
)*
381
radius
.
Radius
;
382
383
/* default of iEmissPower is 2, observing the full sphere */
384
/* this is the usual case the full volume, just the difference in the two vols */
385
/* this will become the true vol when multiplied by the square of the inner radius */
386
radius
.
dVeff
= (vout - vin)/3.*
geometry
.
FillFac
;
387
}
388
else
389
{
390
/* thin cell limit */
391
/* rad_middle_zone is the middle of the zone;
392
* radius is always the OUTER edge of the current zone */
393
radius
.
dVeff
= (rad_middle_zone/
radius
.
rinner
)*
radius
.
drad
*
geometry
.
FillFac
*
394
(
MIN2
(rad_middle_zone,
radius
.
CylindHigh
)/
radius
.
rinner
);
395
}
396
397
/* this is possible correction for slit projected onto resolved object -
398
* this only happens when aperture command is entered.
399
* NB not clear what happens when slit and cylinder are both used -
400
* this would be very orientation specific */
401
402
/* default of iEmissPower is 2, set to 0 is just aperture beam,
403
* and to 1 if aperture slit set */
404
if
(
geometry
.
iEmissPower
== 1 )
405
{
406
/* this is long slit option, remove one power of radius */
407
radius
.
dVeff
/= (rad_middle_zone/
radius
.
rinner
);
408
}
409
else
if
(
geometry
.
iEmissPower
== 0 )
410
{
411
/* this is short slit option, remove two powers of radius */
412
radius
.
dVeff
/=
POW2
(rad_middle_zone/
radius
.
rinner
);
413
}
414
415
/* covering factor, default is unity */
416
radius
.
dVeff
*=
geometry
.
covgeo
;
417
418
/* these are needed for line intensities in outward beam
419
* lgSphere set, geometry.covrt usually 1, 0 when not lgSphere
420
* so outward is 1 when lgSphere set 1/2 when not set, */
421
outwrd = (1. +
geometry
.
covrt
)/2.;
422
423
/*>>>chng 99 apr 23, from above to below */
424
/*radius.dVolOutwrd = outwrd*POW2( (radius.Radius-radius.drad_x_fillfac/2.)/radius.Radius) *
425
radius.drad;*/
426
/* this includes covering fact, the effective vol,, and 1/r^2 dilution,
427
* dReff includes filling factor. the covering factor part is 1 for sphere,
428
* 1/2 for open */
429
/*radius.dVolOutwrd = outwrd*radius.Radius*radius.drad_x_fillfac/(radius.Radius +
430
2.*radius.drad);*/
431
/* dReff from above, includes filling factor */
432
radius
.
dVolOutwrd
= outwrd*
radius
.
drad_x_fillfac
;
433
ASSERT
(
radius
.
dVolOutwrd
> 0. );
434
435
/* following will be used to "uncorrect" for this in lines when predicting continua
436
radius.GeoDil = radius.Radius/(radius.Radius + 2.*radius.drad);*/
437
438
/* this should multiply the line to become the net inward. geo part is 1/2 for
439
* open geometry, 0 for closed. for case of isotropic emitter only this and dVolOutwrd
440
* above are used */
441
radius
.
dVolReflec
= (1. - outwrd)*
radius
.
drad_x_fillfac
*
radius
.
r1r0sq
;
442
443
if
(
geometry
.
lgSphere
)
444
{
445
/* inward beams do not go in when lgSphere set since geometry symmetric */
446
radius
.
BeamInIn
= 0.;
447
radius
.
BeamInOut
=
radius
.
Radius
*
radius
.
drad_x_fillfac
/(
radius
.
Radius
+
448
2.*
radius
.
drad
);
449
}
450
else
451
{
452
radius
.
BeamInIn
=
radius
.
drad_x_fillfac
*
radius
.
r1r0sq
;
453
454
/* inward beams do not go out */
455
radius
.
BeamInOut
= 0.;
456
}
457
/* factor for outwardly directed part of line */
458
radius
.
BeamOutOut
=
radius
.
Radius
*
radius
.
drad_x_fillfac
/(
radius
.
Radius
+
459
2.*
radius
.
drad
);
460
return
;
461
}
462
463
void
ZoneEnd
(
void
)
464
{
465
long
i;
466
467
DEBUG_ENTRY
(
"ZoneEnd()"
);
468
469
/***********************************************************************
470
*
471
* correct rfield for attenuation from center of zone to inner edge
472
*
473
***********************************************************************/
474
475
/* radius is outer radius of this zone, this resets continuum to
476
* flux at illuminated face of zone we have already computed. */
477
478
/* opac.tmn defined in ZoneStart */
479
/* undilute radiation so that rfield is at face of zone */
480
/* NB, upper limit of sum includes [nflux] to test unit verification cell */
481
for
( i=0; i <=
rfield
.
nflux
; i++ )
482
{
483
rfield
.
flux_beam_const
[i] /=
opac
.
tmn
[i];
484
rfield
.
flux_beam_time
[i] /=
opac
.
tmn
[i];
485
rfield
.
flux_isotropic
[i] /=
opac
.
tmn
[i];
486
rfield
.
flux
[i] =
rfield
.
flux_beam_const
[i] +
rfield
.
flux_beam_time
[i] +
487
rfield
.
flux_isotropic
[i];
488
/*rfield.flux[i] /= opac.tmn[i];*/
489
/* >>chng 03 nov 08, update SummedCon here since flux changed */
490
rfield
.
SummedCon
[i] =
rfield
.
flux
[i] +
rfield
.
SummedDif
[i];
491
}
492
493
/* do advection if enabled */
494
if
(
dynamics
.
lgAdvection
)
495
DynaEndZone
();
496
return
;
497
}
Generated for cloudy by
1.8.4