cloudy  trunk
 All Data Structures 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-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  {
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 */
65  pow( radius.Radius/radius.rinner ,(double)geometry.filpow));
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 */
78  /* effective radius */
80 
81  /* reset radius to inner radius, at this point equal to illuminated face radius */
84 
85  /* thickness of next zone */
87 
88  /* depth to illuminated face */
90 
91  /* depth variable for globule */
93 
94  /* this is case where outer radius is set */
95  if( radius.router[iteration-1] < 5e29 )
96  {
98  }
99  else if( iteration > 1 )
100  {
101  /* this case second or later iteration, use previous thickness */
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 */
114  /* >>chng 06 mar 21, remember largest and smallest dr over this iteration
115  * for use in recombination time dep logic */
118 
119  ASSERT( radius.drad > 0. );
120  radius.depth += radius.drad;
122 
123  /* effective radius */
125 
126  /* integral of depth times filling factor */
128 
129  /* decrement Depth2Go but do not let become negative */
131 
132  /* increment the radius, during the calculation Radius is the
133  * outer radius of the current zone*/
135 
136  /* Radius is now outer edge of this zone since incremented above,
137  * so need to add drad to it */
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 */
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(
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];
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 */
212  }
213 
214  /* following is distance to globule, usually does not matter */
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) &&
224  (nzone > IOFF + 1) )
225  {
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 */
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 */
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 */
320  i = 1;
321  while( i < 5 && nzone-i-1>=0 )
322  {
324  ++i;
325  }
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 */
362 
363  /* following only approximate, used for center of next zone */
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)*
377  (radius.Radius - radius.drad);
378 
379  vout = (radius.Radius/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 */
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 */
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 */
442 
443  if( geometry.lgSphere )
444  {
445  /* inward beams do not go in when lgSphere set since geometry symmetric */
446  radius.BeamInIn = 0.;
448  2.*radius.drad);
449  }
450  else
451  {
453 
454  /* inward beams do not go out */
455  radius.BeamInOut = 0.;
456  }
457  /* factor for outwardly directed part of line */
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];
488  /*rfield.flux[i] /= opac.tmn[i];*/
489  /* >>chng 03 nov 08, update SummedCon here since flux changed */
491  }
492 
493  /* do advection if enabled */
494  if( dynamics.lgAdvection )
495  DynaEndZone();
496  return;
497 }

Generated for cloudy by doxygen 1.8.1.1