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