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