cloudy trunk
Loading...
Searching...
No Matches
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
29void 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 */
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();
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}
void FeII_LineZero(void)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
long nWindLine
Definition cdinit.cpp:19
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
realnum & TauIn() const
Definition emission.h:423
realnum & opacity() const
Definition emission.h:593
realnum & Pdest() const
Definition emission.h:543
realnum & TauTot() const
Definition emission.h:433
TransitionProxy::iterator iterator
Definition transition.h:280
EmissionList::reference Emis() const
Definition transition.h:408
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long ipoint(double energy_ryd)
t_conv conv
Definition conv.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
void SumDensities(void)
Definition dense.cpp:200
realnum GetDopplerWidth(realnum massAMU)
t_geometry geometry
Definition geometry.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipHe1s1S
Definition iso.h:41
const int ipH3p
Definition iso.h:31
const int ipH1s
Definition iso.h:27
const int ipHe2p1P
Definition iso.h:49
#define N_(A_)
Definition iso.h:20
const int ipHE_LIKE
Definition iso.h:63
const int ipH3s
Definition iso.h:30
const int ipH2p
Definition iso.h:29
const int ipH3d
Definition iso.h:32
const int ipH2s
Definition iso.h:28
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
#define S(I_, J_)
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
void RT_line_all(void)
#define TAULIM
void RT_tau_init(void)
t_StopCalc StopCalc
Definition stopcalc.cpp:5
long int nSpecies
Definition taulines.cpp:21
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
vector< vector< TransitionList > > ExtraLymanLines
Definition taulines.cpp:25
TransitionList UTALines("UTALines", &AnonStates)
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nUTA
Definition taulines.cpp:26
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
multi_arr< int, 3 > ipExtraLymanLines
Definition taulines.cpp:24
TransitionList TauLines("TauLines", &AnonStates)
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5