cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_tau_init.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 /*RT_tau_init set initial outward optical depths at start of first iteration,
4  * it is only called by cloudy one time per complete calculation, just after
5  * continuum set up and before start of convergence attempts.
6  *
7  * RT_tau_reset after first iteration, updates the optical depths, mirroring this
8  * routine but with the previous iteration's variables */
9 #include "cddefines.h"
10 #define TAULIM 1e8
11 #include "taulines.h"
12 #include "doppvel.h"
13 #include "iso.h"
14 #include "h2.h"
15 #include "lines_service.h"
16 #include "rfield.h"
17 #include "dense.h"
18 #include "opacity.h"
19 #include "thermal.h"
20 #include "geometry.h"
21 #include "stopcalc.h"
22 #include "ipoint.h"
23 #include "abund.h"
24 #include "conv.h"
25 #include "atomfeii.h"
26 #include "rt.h"
27 #include "trace.h"
28 
29 void RT_tau_init(void)
30 {
31  long int i,
32  nelem,
33  ipISO,
34  ipHi,
35  ipLo,
36  nHi;
37  long lgHit; /* used for logic check */
38 
39  double AbunRatio,
40  balc,
41  coleff;
42 
43  realnum f, BalmerTauOn;
44 
45  DEBUG_ENTRY( "RT_tau_init()" );
46 
47  ASSERT( dense.eden > 0. );
48 
49  /* Zero lines first */
50  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
51  {
52  for( nelem=ipISO; nelem < LIMELM; nelem++ )
53  {
54  if( dense.lgElmtOn[nelem] )
55  {
56  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
57  {
58  if( iso.lgDielRecom[ipISO] )
59  TransitionZero( &SatelliteLines[ipISO][nelem][ipHi] );
60 
61  for( ipLo=0; ipLo < ipHi; ipLo++ )
62  {
63  TransitionZero( &Transitions[ipISO][nelem][ipHi][ipLo] );
64  }
65  }
66  for( ipHi=2; ipHi <iso.nLyman[ipISO]; ipHi++ )
67  {
68  TransitionZero( &ExtraLymanLines[ipISO][nelem][ipHi] );
69  }
70  }
71  }
72  }
73 
74  /* initialize heavy element line array */
75  for( i=1; i <= nLevel1; i++ )
76  {
77  TransitionZero( &TauLines[i] );
78  }
79  /* this is a dummy optical depth array for non-existant lines
80  * when this goes over to struc, make sure all are set to zero here since
81  * init in grid may depend on it */
83 
84  /* lines in cooling function with Mewe approximate collision strengths */
85  for( i=0; i < nWindLine; i++ )
86  {
87  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
88  {
89  /* inward optical depth */
90  TransitionZero( &TauLine2[i] );
91  }
92  }
93 
94  /* inner shell lines */
95  for( i=0; i < nUTA; i++ )
96  {
97  /* these are line optical depth arrays
98  * inward optical depth */
99  /* heat is special for this array - it is heat per pump */
100  double hsave = UTALines[i].Coll.heat;
101  TransitionZero( &UTALines[i] );
102  UTALines[i].Coll.heat = hsave;
103  }
104 
105  /* hyperfine structure lines */
106  for( i=0; i < nHFLines; i++ )
107  {
108  TransitionZero( &HFLines[i] );
109  }
110 
111  /*Atoms & Molecules*/
112  for( i=0; i < linesAdded2; i++)
113  {
114  TransitionZero( atmolEmis[i].tran );
115  }
116 
117  /* CO carbon monoxide lines */
118  for( i=0; i < nCORotate; i++ )
119  {
120  /* >>chng 03 feb 14, use main zero function */
123  }
124 
125  /* initialize optical depths in H2 */
126  H2_LineZero();
127 
128  /* initialize large atom FeII arrays */
129  FeII_LineZero();
130 
131  /* ==================================================================*/
132  /* end setting lines to zero */
133 
134  /* >>chng 02 feb 08, moved to here from opacitycreateall, which was called later.
135  * bug inhibited convergence in some models. Caught by PvH */
136  /* this is option to stop at certain optical depth */
137  if( StopCalc.taunu > 0. )
138  {
141  }
142  else
143  {
145  /* >>chng 04 dec 21, remove from here and init to 1e30 in zero */
146  /*StopCalc.tauend = 1e30f;*/
147  }
148 
149  /* if taunu was set with a stop optical depth command, then iptnu must
150  * have also been set to a continuum index before this code is reaced */
151  ASSERT( StopCalc.taunu == 0. || StopCalc.iptnu >= 0 );
152 
153  /* set initial and total optical depths in arrays
154  * TAUNU is set when lines read in, sets stopping radius */
155  if( StopCalc.taunu > 0. )
156  {
157  /* an optical depth has been specified */
159  {
160  /* at ionizing energies */
161  for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
162  {
163  /* taumin set to 1e-20, can be reset with taumin command */
164  opac.TauAbsGeo[1][i] = opac.taumin;
165  opac.TauScatGeo[1][i] = opac.taumin;
166  opac.TauTotalGeo[1][i] = opac.taumin;
167  }
168 
169  for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
170  {
171  /* TauAbsGeo(i,2) = tauend * (anu(i)/anu(iptnu))**(-2.43) */
172  opac.TauAbsGeo[1][i] = StopCalc.tauend;
173  opac.TauScatGeo[1][i] = opac.taumin;
174  opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
175  }
176  }
177 
178  else
179  {
180  /* not specified at ionizing energies */
181  for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
182  {
183  opac.TauAbsGeo[1][i] = StopCalc.tauend;
184  opac.TauScatGeo[1][i] = StopCalc.tauend;
186  }
187 
188  for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
189  {
190  opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f));
191  opac.TauScatGeo[1][i] = opac.taumin;
192  opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
193  }
194 
195  }
196  }
197 
198  else
199  {
200  /* ending optical depth not specified, assume 1E8 at 1 Ryd */
201  for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
202  {
203  opac.TauAbsGeo[1][i] = opac.taumin;
204  opac.TauScatGeo[1][i] = opac.taumin;
205  opac.TauTotalGeo[1][i] = opac.taumin;
206  }
207 
208  for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
209  {
210  opac.TauAbsGeo[1][i] = (realnum)(TAULIM*pow(rfield.anu[i],(realnum)-2.43f));
211  opac.TauScatGeo[1][i] = opac.taumin;
212  opac.TauTotalGeo[1][i] = opac.TauAbsGeo[1][i] + opac.TauScatGeo[1][i];
213  }
214  }
215 
216  /* if lgSphere then double outer, set inner to half
217  * assume will be optically thin at sub-ionizing energies */
218  if( geometry.lgSphere || opac.lgCaseB )
219  {
220  for( i=0; i < (iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - 1); i++ )
221  {
222  opac.TauAbsGeo[0][i] = opac.taumin;
223  opac.TauAbsGeo[1][i] = opac.taumin*2.f;
224  opac.TauScatGeo[0][i] = opac.taumin;
225  opac.TauScatGeo[1][i] = opac.taumin*2.f;
226  opac.TauTotalGeo[0][i] = 2.f*opac.taumin;
227  opac.TauTotalGeo[1][i] = 4.f*opac.taumin;
228  }
229 
230  for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nupper; i++ )
231  {
232  opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i];
233  opac.TauAbsGeo[1][i] *= 2.;
234  opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i];
235  opac.TauScatGeo[1][i] *= 2.;
236  opac.TauTotalGeo[0][i] = opac.TauTotalGeo[1][i];
237  opac.TauTotalGeo[1][i] *= 2.;
238  }
239 
240  if( StopCalc.taunu > 0. )
241  {
242  /* ending optical depth specified, and lgSphere also */
243  StopCalc.tauend *= 2.;
244  }
245  }
246 
247  /* fix escape prob for He, metals, first set log of Te, needed by RECEFF
248  * do not do this if temperature has been set by command */
250  {
251  double TeNew;
252  /* this is a typical temperature for the H+ zone, and will use it is
253  * the line widths & opt depth in the following. this is not meant to be the first
254  * temperature during the search phase. */
255  TeNew = 2e4;
256  /* >>chng 05 jul 19, in PDR models the guess of the optical depth in Lya could
257  * be very good since often limiting column density is set, but must use
258  * the temperature of H0 gas. in a dense BLR this is roughly 7000K and
259  * closer to 100K for a low-density PDR - use these guesses */
260  if( dense.gas_phase[ipHYDROGEN] >= 1e9 )
261  {
262  /* this is a typical BLR H0 temp */
263  TeNew = 7000.;
264  }
265  else if( dense.gas_phase[ipHYDROGEN] <= 1e5 )
266  {
267  /* this is a typical PDR H0 temp */
268  TeNew = 100.;
269  }
270  else
271  {
272  /* power law interpolation between them */
273  TeNew = 0.5012 * pow( (double)dense.gas_phase[ipHYDROGEN], 0.46 );
274  }
275 
276  /* propagate this temperature through the code */
277  /* must not call tfidle at this stage since not all vectors have been allocated */
278  TempChange(TeNew );
279  }
280 
281  /* set inward optical depths for hydrogenic ions to small number proportional to abundance */
282  for( nelem=0; nelem < LIMELM; nelem++ )
283  {
284  if( dense.lgElmtOn[nelem] )
285  {
286  /* now get actual optical depths */
287  AbunRatio = abund.solar[nelem]/abund.solar[0];
288  for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
289  {
290  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
291  {
292  /* set all inward optical depths to taumin, regardless of abundance
293  * this is a very small number, 1e-20 */
294  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin;
295  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin;
296  }
297  }
298 
299  /* La may be case B, tlamin set to 1e9 by default with case b command */
300  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin;
301 
302  /* scale factor so that all other lyman lines are appropriate for this Lya optical depth*/
303  f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
304  fixit(); /* this appears to be redundant to code below. */
305 
306  for( nHi=3; nHi<=iso.n_HighestResolved_max[ipH_LIKE][nelem]; nHi++ )
307  {
308  ipHi = iso.QuantumNumbers2Index[ipH_LIKE][nelem][nHi][1][2];
309  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin,
310  f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity );
311  }
312  for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
313  {
314  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2( opac.taumin,
315  f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity );
316  }
317 
318  /* after this set of if's the total lya optical depth will be known,
319  * common code later on will set rest of lyman lines
320  * if case b then set total lyman to twice inner */
321  if( opac.lgCaseB )
322  {
323  /* force outer optical depth to twice inner if case B */
324  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot =
325  (realnum)(2.*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn);
326  /* force off Balmer et al optical depths */
327  BalmerTauOn = 0.f;
328  }
329 
330  else
331  {
332  /* usual case for H LyA; try to guess outer optical depth */
333  if( StopCalc.colnut < 6e29 )
334  {
335  /* neutral column is defined */
336  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.colnut*
337  rt.DoubleTau*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity/
338  DoppVel.doppler[nelem]*AbunRatio);
339  ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
340  }
341  /* has optical depth at some energy where we want to stop been specified?
342  * taunu is energy where
343  * this has been specified - this checks on Lyman continuum, taunu = 1 */
344  else if( StopCalc.taunu < 3. && StopCalc.taunu >= 0.99 )
345  {
346  /* Lyman continuum optical depth defined */
347  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(StopCalc.tauend*
348  1.2e4*1.28e6/DoppVel.doppler[nelem]*rt.DoubleTau*AbunRatio);
349  ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
350  }
351  else if( StopCalc.HColStop < 6e29 )
352  {
353  /* neutral column not defined, guess from total col and U */
354  coleff = StopCalc.HColStop -
355  MIN2(MIN2(rfield.qhtot/dense.eden,1e24)/2.6e-13,0.6*StopCalc.HColStop);
356 
357  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(coleff*
358  7.6e-14*AbunRatio);
359  ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
360  }
361  else
362  {
363  /* no way to estimate 912 optical depth, set to large number */
364  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot = (realnum)(1e20*
365  AbunRatio);
366  ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 0. );
367  }
368  /* allow Balmer et al. optical depths */
369  BalmerTauOn = 1.f;
370  }
371 
372  /* Lya total optical depth now known, is it optically thick?*/
373  if( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot > 1. )
374  {
375  rt.TAddHLya = (realnum)MIN2(1.,Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/
376  1e4);
377  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn += rt.TAddHLya;
378  }
379 
380  else
381  {
383  }
384 
385  /* this scale factor is to set other lyman lines, given the Lya optical depth */
386  f = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot/
387  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
388 
389  ipISO = ipH_LIKE;
390  ASSERT( ipISO<NISO && nelem < LIMELM );
391  for( nHi=3; nHi<=iso.n_HighestResolved_max[ipISO][nelem]; nHi++ )
392  {
393  ipHi = iso.QuantumNumbers2Index[ipISO][nelem][nHi][1][2];
394  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin,
395  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f );
396 
397  /* increment inward optical depths by rt all lyman lines, inward
398  * optical depth was set above, this adds to it. this was originally
399  * set some place above */
400  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya*
401  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/
402  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
403 
404  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2(
405  opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn );
406  }
407  for( ipHi=iso.numLevels_max[ipH_LIKE][nelem] - iso.nCollapsed_max[ipH_LIKE][nelem]; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
408  {
409  /* set total optical depth for higher lyman lines */
410  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot = MAX2( opac.taumin,
411  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity * f );
412 
413  /* increment inward optical depths by rt all lyman lines, inward
414  * optical depth was set above, this adds to it. this was originally
415  * set some place above */
416  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn += rt.TAddHLya*
417  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity/
418  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
419 
420  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn = MAX2(
421  opac.taumin, Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn );
422  }
423 
424  /* try to guess what Balmer cont optical guess,
425  * first branch is case where we will stop at a balmer continuum optical
426  * depth - taunu is energy where tauend was set */
427  if( StopCalc.taunu > 0.24 && StopCalc.taunu < 0.7 )
428  {
429  Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot = (realnum)(StopCalc.tauend*
430  3.7e4*BalmerTauOn*AbunRatio + 1e-20);
431 
432  Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend*
433  3.7e4*BalmerTauOn*AbunRatio + 1e-20);
434 
435  Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot = (realnum)(StopCalc.tauend*
436  3.7e4*BalmerTauOn*AbunRatio + 1e-20);
437  }
438 
439  else
440  {
441  /* this is a guess based on Ferland&Netzer 1979, but it gets very large */
442  balc = rfield.qhtot*2.1e-19*BalmerTauOn*AbunRatio + 1e-20;
443 
444  Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot =
445  (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
446 
447  Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot =
448  (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
449 
450  Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot =
451  (realnum)(MIN2(2e5, balc*3.7e4*BalmerTauOn+1e-20));
452 
453  ASSERT( Transitions[ipH_LIKE][nelem][ipH3p][ipH2s].Emis->TauTot > 0.);
454  ASSERT( Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot > 0.);
455  ASSERT( Transitions[ipH_LIKE][nelem][ipH3d][ipH2p].Emis->TauTot > 0.);
456  }
457 
458  /* 2s - 1s strongly forbidden */
459  Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauTot = 2.f*opac.taumin;
460  Transitions[ipH_LIKE][nelem][ipH2s][ipH1s].Emis->TauIn = opac.taumin;
461 
462  /* transitions down to 2s are special since 'alpha' (2s-2p) has
463  * small optical depth, so use 3 to 2p instead */
464  f = Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->TauTot/
465  Transitions[ipH_LIKE][nelem][ipH3s][ipH2p].Emis->opacity* BalmerTauOn;
466 
467  ipLo = ipH2s;
468  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
469  {
470  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
471  continue;
472 
473  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin +
474  f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity;
475  ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.);
476  }
477 
478  /* this is for rest of lines, scale from opacity */
479  for( ipLo=ipH2p; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
480  {
481  long ipISO = ipH_LIKE, ipNS, ipNPlusOneP;
482 
483  /* scale everything with same factor we use for n+1 P -> nS */
484  ipNS = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo) ][0][2];
485  if( ( N_(ipLo) + 1 ) <=
487  {
488  ipNPlusOneP = iso.QuantumNumbers2Index[ipH_LIKE][nelem][ N_(ipLo)+1 ][1][2];
489  }
490  else
491  {
492  ipNPlusOneP = ipNS + 1;
493  }
494 
495 #if 1 /* old way */
496  ipNS = ipLo;
497  ipNPlusOneP = ipNS + 1;
498 #endif
499 
500  if( Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].ipCont <= 0 )
501  {
502  f = SMALLFLOAT;
503  }
504  else
505  {
506  f = Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->TauTot/
507  Transitions[ipH_LIKE][nelem][ipNPlusOneP][ipNS].Emis->opacity*
508  BalmerTauOn;
509  }
510 
511  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
512  {
513  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
514  continue;
515 
516  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot = opac.taumin +
517  f* Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->opacity;
518  ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot > 0.);
519  }
520  }
521 
522  /* this loop is over all possible H lines, do some final cleanup */
523  for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
524  {
525  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
526  {
527  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
528  continue;
529 
530  /* TauCon is line optical depth to inner face used for continuum pumping rate,
531  * not equal to TauIn in case of static sphere since TauCon will never
532  * count the far side line optical depth */
533  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauCon =
534  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn;
535 
536  /* make sure inward optical depth is not larger than half total */
537  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn =
538  MIN2( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn ,
539  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot/2.f );
540  ASSERT(Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn > 0.f);
541 
542  /* this is fraction of line that goes inward, not known until second iteration*/
543  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->FracInwd = 0.5;
544  }
545  }
546  }
547  }
548 
549  /* initialize all he-like optical depths */
550  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
551  {
552  if( dense.lgElmtOn[nelem] )
553  {
554  for( ipLo=0; ipLo < (iso.numLevels_max[ipHE_LIKE][nelem] - 1); ipLo++ )
555  {
556  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
557  {
558  /* set all inward optical depths to taumin, regardless of abundance
559  * this is a very small number, 1e-20 */
560  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauIn = opac.taumin;
561  Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauTot = 2.f*opac.taumin;
562  }
563  }
564  }
565  }
566 
567  /* now do helium like sequence if case b */
568  if( opac.lgCaseB )
569  {
570  for( nelem=1; nelem < LIMELM; nelem++ )
571  {
572  if( dense.lgElmtOn[nelem] )
573  {
574  double Aprev;
575  realnum ratio;
576  /* La may be case B, tlamin set to 1e9 by default with case b command */
577  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin;
578 
579  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon =
580  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
581 
582  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot =
583  2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
584 
585  ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity;
586 
587  /* this will be the trans prob of the previous lyman line, will use this to
588  * find the next one up in the series */
589  Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul;
590 
591  i = ipHe2p1P+1;
592  /* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess
593  * which are which - this will work for any number of levels */
594  for( i=ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
595  {
596  if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
597  continue;
598 
599  /* >>chng 02 mar 15, add test for collapsed levels - As are very much
600  * smaller at boundary between collapsed/resolved so this test is needed*/
601  if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. ||
602  StatesElem[ipHE_LIKE][nelem][i].S < 0 )
603  {
604  Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul;
605 
606  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn =
607  ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity;
608  /* reset line optical depth to continuum source */
609  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon =
610  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
611  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot =
612  2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
613  }
614  }
615  }
616  }
617  }
618 
619  /* begin sanity check, total greater than 1/0.9 time inward */
620  lgHit = false;
621  for( nelem=0; nelem < LIMELM; nelem++ )
622  {
623  if( dense.lgElmtOn[nelem] )
624  {
625  for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
626  {
627  for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
628  {
629  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
630  continue;
631 
632  if( Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot <
633  0.9f*Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn )
634  {
635  fprintf(ioQQQ,
636  "RT_tau_init insanity for h line, Z=%li lo=%li hi=%li tot=%g in=%g \n",
637  nelem , ipLo, ipHi , Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauTot ,
638  Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->TauIn );
639  lgHit = true;
640  }
641  }
642  }
643  }
644  }
645  if( lgHit )
646  {
647  fprintf( ioQQQ," stopping due to insanity in RT_tau_init\n");
648  ShowMe();
649  cdEXIT(EXIT_FAILURE);
650  }
651  /*end sanity check */
652 
653  /* fix offset for effective column density optical depth */
654  rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
655 
656  /* this is flag detected by dest prob routines to see whether ots rates are
657  * oscaillating - damp them out if so */
658  conv.lgOscilOTS = false;
659 
660  /* now that optical depths have been incremented, do escape prob, this
661  * is located here instead on in cloudy.c (where it belongs) because
662  * information generated by RT_line_all is needed for the following printout */
663 
664  RT_line_all(true , false );
665 
666  /* rest of routine is printout in case of trace */
667  if( trace.lgTrace )
668  {
669  /* default is h-like */
670  ipISO = ipH_LIKE;
672  ipISO = ipHE_LIKE;
673 
674  if( trace.lgIsoTraceFull[ipISO] )
675  {
676  fprintf( ioQQQ, "\n\n up Transitions[ipISO][hi][lo].TauTot array as set in RT_tau_init ipZTrace (nelem)= %ld\n",
678  for( ipHi=2; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
679  {
680  fprintf( ioQQQ, " %3ld", ipHi );
681  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
682  {
683  if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
684  continue;
685 
686  fprintf( ioQQQ,PrintEfmt("%9.2e",
687  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauTot ));
688  }
689  fprintf( ioQQQ, "\n" );
690  }
691 
692  fprintf( ioQQQ, "\n\n TauIn array\n" );
693  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
694  {
695  fprintf( ioQQQ, " %3ld", ipHi );
696  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
697  {
698  if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
699  continue;
700 
701  fprintf( ioQQQ,PrintEfmt("%9.2e",
702  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->TauIn ));
703  }
704  fprintf( ioQQQ, "\n" );
705  }
706 
707  fprintf( ioQQQ, "\n\n Aul As array\n" );
708  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
709  {
710  fprintf( ioQQQ, " %3ld", ipHi );
711  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
712  {
713  if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
714  continue;
715 
716  fprintf( ioQQQ,PrintEfmt("%9.2e",
717  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul) );
718  }
719  fprintf( ioQQQ, "\n" );
720  }
721 
722  fprintf( ioQQQ, "\n\n Aul*esc array\n" );
723  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
724  {
725  fprintf( ioQQQ, " %3ld", ipHi );
726  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
727  {
728  if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
729  continue;
730 
731  fprintf( ioQQQ,PrintEfmt("%9.2e",
732  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul*
733  (Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pdest +
734  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pelec_esc +
735  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Pesc) ));
736  }
737  fprintf( ioQQQ, "\n" );
738  }
739 
740  fprintf( ioQQQ, "\n\n H opac array\n" );
741  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][trace.ipIsoTrace[ipISO]]; ipHi++ )
742  {
743  fprintf( ioQQQ, " %3ld", ipHi );
744  for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
745  {
746  if( Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->Aul <= iso.SmallA )
747  continue;
748 
749  fprintf( ioQQQ,PrintEfmt("%9.2e",
750  Transitions[ipISO][trace.ipIsoTrace[ipISO]][ipHi][ipLo].Emis->opacity ));
751  }
752  fprintf( ioQQQ, "\n" );
753  }
754  }
755 
756  else
757  {
758  fprintf( ioQQQ, " RT_tau_init called.\n" );
759  }
760  }
761 
762  ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauIn> 0. );
763  ASSERT( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc>= 0. );
764  return;
765 }

Generated for cloudy by doxygen 1.8.3.1