cloudy trunk
Loading...
Searching...
No Matches
rt_tau_reset.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_reset after first iteration, updates the optical depths, mirroring this
4 * routine but with the previous iteration's variables */
5#include "cddefines.h"
6#include "taulines.h"
7#include "trace.h"
8#include "iso.h"
9#include "rfield.h"
10#include "opacity.h"
11#include "h2.h"
12#include "mole.h"
13#include "geometry.h"
14#include "dense.h"
15#include "atomfeii.h"
16#include "colden.h"
17#include "rt.h"
18
19/* ====================================================================== */
20/*RT_tau_reset update total optical depth scale,
21 * called by cloudy after iteration is complete */
22void RT_tau_reset(void)
23{
24 long int i,
25 ipHi,
26 ipISO,
27 nelem,
28 ipLo;
29
30 DEBUG_ENTRY( "RT_tau_reset()" );
31
32 if( trace.lgTrace )
33 {
34 fprintf( ioQQQ, " UPDATE estimating new optical depths\n" );
35 if( trace.lgHBug && trace.lgIsoTraceFull[ipH_LIKE] )
36 {
37 fprintf( ioQQQ, " New Hydrogen outward optical depths:\n" );
38 for( ipHi=1; ipHi < iso_sp[ipH_LIKE][ trace.ipIsoTrace[ipH_LIKE] ].numLevels_max; ipHi++ )
39 {
40 fprintf( ioQQQ, "%3ld", ipHi );
41 for( ipLo=0; ipLo < ipHi; ipLo++ )
42 {
43 if( iso_sp[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]].trans(ipHi,ipLo).ipCont() <= 0 )
44 fprintf( ioQQQ, "%10.2e", 1e-30 );
45 else
46 fprintf( ioQQQ, "%10.2e",
47 iso_sp[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]].trans(ipHi,ipLo).Emis().TauIn() );
48 }
49 fprintf( ioQQQ, "\n" );
50 }
51 }
52 }
53
54 /* save column densities of H species */
55 for( i=0; i<NCOLD; ++i )
56 {
57 colden.colden_old[i] = colden.colden[i];
58 }
59 for( i=0; i < mole_global.num_calc; i++ )
60 {
61 mole.species[i].column_old = mole.species[i].column;
62 }
63
64 opac.telec = opac.taumin;
65 opac.thmin = opac.taumin;
66
67 /* all iso sequences */
68 for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
69 {
70 for( nelem=ipISO; nelem < LIMELM; nelem++ )
71 {
72 if( dense.lgElmtOn[nelem] )
73 {
74 if( iso_ctrl.lgDielRecom[ipISO] )
75 {
76 for( ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
77 {
78 RT_line_one_tau_reset(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipHi]]);
79 }
80 }
81
82 for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
83 {
84 /* the bound-bound transitions within the model atoms */
85 for( ipLo=0; ipLo < ipHi; ipLo++ )
86 {
87 enum {DEBUG_LOC=false};
88 if( DEBUG_LOC )
89 {
90 if( ipISO==1 && nelem==1 && ipHi==iso_ctrl.nLyaLevel[ipISO] && ipLo==0 )
91 fprintf(ioQQQ,"DEBUG rt before loop %li %li %li %li tot %.3e in %.3e\n",
92 ipISO, nelem, ipHi , ipLo ,
93 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
94 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
95 }
96 /*RT_line_one_tau_reset computes average of old and new optical depths
97 * for new scale at end of iter */
98 RT_line_one_tau_reset(iso_sp[ipISO][nelem].trans(ipHi,ipLo));
99 if( DEBUG_LOC )
100 {
101 if( ipISO==1 && nelem==1 && ipHi==iso_ctrl.nLyaLevel[ipISO] && ipLo==0 )
102 fprintf(ioQQQ,"DEBUG rt after loop %li %li %li %li tot %.3e in %.3e\n",
103 ipISO, nelem, ipHi , ipLo ,
104 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauTot() ,
105 iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Emis().TauIn());
106 }
107 }
108 }
109 /* the extra Lyman lines */
110 for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
111 {
112 /* fully transfer all of the extra lines even though
113 * have not solved for their upper level populations */
114 RT_line_one_tau_reset(ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]]);
115 }
116 }
117 }
118 }
119
120 /* >>>chng 99 nov 11 did not have case b for hydrogenic species on second and
121 * higher iterations */
122 /* option to clobber these taus for Lyman lines, if case b is set */
123 if( opac.lgCaseB )
124 {
125 for( nelem=0; nelem < LIMELM; nelem++ )
126 {
127 if( dense.lgElmtOn[nelem] )
128 {
129 realnum f;
130 /* La may be case B, tlamin set to 1e9 by default with case b command */
131 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn() = opac.tlamin;
132 /* >>>chng 99 nov 22, did not reset TauCon */
133 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauCon() = iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
134 iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauTot() =
135 2.f*iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
136 f = opac.tlamin/iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().opacity();
137
138 for( ipHi=3; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_max; ipHi++ )
139 {
140 if( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).ipCont() <= 0 )
141 continue;
142
143 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauIn() =
144 f*iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().opacity();
145 /* reset line optical depth to continuum source */
146 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauCon() = iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauIn();
147 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipH1s).Emis().TauTot() =
148 2.f*iso_sp[ipH_LIKE][nelem].trans(ipH2p,ipH1s).Emis().TauIn();
149 }
150 }
151 }
152
153 /* now do helium like sequence - different since collapsed levels
154 * all go to ground */
155 for( nelem=1; nelem < LIMELM; nelem++ )
156 {
157 if( dense.lgElmtOn[nelem] )
158 {
159 double Aprev;
160 realnum ratio;
161 /* La may be case B, tlamin set to 1e9 by default with case b command */
162 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn() = opac.tlamin;
163
164 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauCon() = iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn();
165
166 iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauTot() =
167 2.f*iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().TauIn();
168
169 ratio = opac.tlamin/iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().opacity();
170
171 /* this will be the trans prob of the previous lyman line, will use this to
172 * find the next one up in the series */
173 Aprev = iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).Emis().Aul();
174
175 i = ipHe2p1P+1;
176 /* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess
177 * which are which - this will work for any number of levels */
178 for( i = ipHe2p1P+1; i < iso_sp[ipHE_LIKE][nelem].numLevels_max; i++ )
179 {
180 if( iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).ipCont() <= 0 )
181 continue;
182
183 /* >>chng 02 mar 19 use proper test for resonance collapsed lines */
184 if( iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().Aul()> Aprev/10. ||
185 iso_sp[ipHE_LIKE][nelem].st[i].S() < 0 )
186 {
187 Aprev = iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().Aul();
188 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn() =
189 ratio*iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().opacity();
190 /* reset line optical depth to continuum source */
191 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauCon() =
192 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn();
193 iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauTot() =
194 2.f*iso_sp[ipHE_LIKE][nelem].trans(i,ipHe1s1S).Emis().TauIn();
195 }
196 }
197 }
198 }
199 }
200
201 /* all heavy element lines */
202 for( i=1; i <= nLevel1; i++ )
203 {
205 /* >>chng 96 jul 06 following sanity check added */
206 ASSERT( fabs(TauLines[i].Emis().TauIn()) <= fabs(TauLines[i].Emis().TauTot()) );
207 }
208
209 /* zero out fine opacity fine grid fine mesh array */
210 memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
211 // Make sure stale zone opacities don't wrap around between
212 // iterations. If this /does/ have an effect, it means that the
213 // opacities being used in are one zone stale. Likely touch points
214 // are line overlap & radiation pressure calculations.
215 // So zero this isn't ideal, just better than using whatever the
216 // value happens to be in the last zone of the previous iteration...
217 memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
218
219 /* all level 2 heavy element lines */
220 for( i=0; i < nWindLine; i++ )
221 {
222 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
223 {
225 }
226 }
227
228 /* all hyperfine structure lines */
229 for( i=0; i < nHFLines; i++ )
230 {
232 }
233
234 /* external database lines */
235 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
236 {
237 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
238 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
239 {
240 RT_line_one_tau_reset((*em).Tran());
241 }
242 }
243
244 /* inner shell lines */
245 for( i=0; i < nUTA; i++ )
246 {
247 /* these are line optical depth arrays
248 * inward optical depth */
249 /* heat is special for this array - it is heat per pump */
250 double hsave = UTALines[i].Coll().heat();
252 UTALines[i].Coll().heat() = hsave;
253 }
254
255 /* the large H2 molecule */
256 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
257 (*diatom)->H2_RT_tau_reset();
258
259 /* large FeII atom */
261
262 if( opac.lgCaseB )
263 {
264 for( i=0; i < rfield.nupper; i++ )
265 {
266 /* DEPABS and SCT are abs and sct optical depth for depth only
267 * we will not change total optical depths, just reset inner to half
268 * TauAbsGeo(i,2) = 2.*TauAbsFace(i) */
269 opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
270 /* TauScatGeo(i,2) = 2.*TauScatFace(i) */
271 opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
272 }
273 }
274 else if( geometry.lgSphere )
275 {
276 /* closed or spherical case */
277 for( i=0; i < rfield.nupper; i++ )
278 {
279 /* [1] is total optical depth from previous iteration,
280 * [0] is optical depth at current position */
281 opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
282 opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
283 opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
284 opac.TauScatGeo[0][i] = opac.TauScatFace[i];
285 opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
286 opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
287 }
288 }
289 else
290 {
291 /* open geometry */
292 for( i=0; i < rfield.nupper; i++ )
293 {
294 opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
295 opac.TauTotalGeo[0][i] = opac.taumin;
296 opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
297 opac.TauAbsGeo[0][i] = opac.taumin;
298 opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
299 opac.TauScatGeo[0][i] = opac.taumin;
300 }
301 }
302
303 /* same for open and closed geometries */
304 for( i=0; i < rfield.nupper; i++ )
305 {
306 /* total optical depth across computed shell */
307 opac.TauAbsTotal[i] = opac.TauAbsFace[i];
308 /* e2( tau across shell), optical depth from ill face to shielded face of cloud
309 * not that opac.TauAbsFace is reset to small number just after this */
310 opac.E2TauAbsTotal[i] = (realnum)e2( opac.TauAbsTotal[i] );
311 /* TauAbsFace and TauScatFace are abs and sct optical depth to ill face */
312 opac.TauScatFace[i] = opac.taumin;
313 opac.TauAbsFace[i] = opac.taumin;
314 }
315
316 /* this is optical depth at x-ray point defining effective optical depth */
317 rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
318 return;
319}
void FeII_RT_tau_reset(void)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
double e2(double t)
Definition service.cpp:299
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
EmissionProxy::iterator iterator
Definition emission.h:317
t_colden colden
Definition colden.cpp:5
#define NCOLD
Definition colden.h:9
t_dense dense
Definition dense.cpp:24
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 ipH1s
Definition iso.h:27
const int ipHe2p1P
Definition iso.h:49
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
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_one_tau_reset(const TransitionProxy &t)
void RT_tau_reset(void)
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)
t_trace trace
Definition trace.cpp:5