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