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