cloudy  trunk
 All Data Structures 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-2008 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 */
22 void RT_tau_reset(void)
23 {
24  long int i,
25  ipHi,
26  ipISO,
27  nelem,
28  ipLo;
29 
30  double WeightNew;
31 
32  DEBUG_ENTRY( "RT_tau_reset()" );
33 
34  if( trace.lgTrace )
35  {
36  fprintf( ioQQQ, " UPDATE estimating new optical depths\n" );
38  {
39  fprintf( ioQQQ, " New Hydrogen outward optical depths:\n" );
40  for( ipHi=1; ipHi < iso.numLevels_max[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]]; ipHi++ )
41  {
42  fprintf( ioQQQ, "%3ld", ipHi );
43  for( ipLo=0; ipLo < ipHi; ipLo++ )
44  {
45  if( Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].ipCont <= 0 )
46  fprintf( ioQQQ, "%10.2e", 1e-30 );
47  else
48  fprintf( ioQQQ, "%10.2e",
49  Transitions[ipH_LIKE][trace.ipIsoTrace[ipH_LIKE]][ipHi][ipLo].Emis->TauIn );
50  }
51  fprintf( ioQQQ, "\n" );
52  }
53  }
54  }
55 
56  /* pumping of CaH */
57  opac.tpcah[1] = opac.tpcah[0];
58  opac.tpcah[0] = opac.taumin;
59 
60  /* save column densities of H species */
61  for( i=0; i<NCOLD; ++i )
62  {
64  }
65  for( i=0; i < mole.num_comole_calc; i++ )
66  {
67  COmole[i]->hevcol_old = COmole[i]->hevcol;
68  }
69 
70  /* ======================================================================== */
71  /* must take average of old and new optical depths - were old in place? */
72  if( iteration <= 1 )
73  {
74  /* this is first pass */
75  WeightNew = 1.;
76  }
77  else
78  {
79  WeightNew = 0.75;
80  }
81 
84 
85  /* all iso sequences */
86  for(ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
87  {
88  for( nelem=ipISO; nelem < LIMELM; nelem++ )
89  {
90  if( dense.lgElmtOn[nelem] )
91  {
92  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
93  {
94  if( iso.lgDielRecom[ipISO] )
95  RT_line_one_tau_reset(&SatelliteLines[ipISO][nelem][ipHi],WeightNew);
96 
97  /* the bound-bound transitions within the model atoms */
98  for( ipLo=0; ipLo < ipHi; ipLo++ )
99  {
100  enum {DEBUG_LOC=false};
101  if( DEBUG_LOC )
102  {
103  if( ipISO==1 && nelem==1 && ipHi==iso.nLyaLevel[ipISO] && ipLo==0 )
104  fprintf(ioQQQ,"DEBUG rt before loop %li %li %li %li tot %.3e in %.3e\n",
105  ipISO, nelem, ipHi , ipLo ,
106  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
107  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
108  }
109  /*RT_line_one_tau_reset computes average of old and new optical depths
110  * for new scale at end of iter */
111  RT_line_one_tau_reset(&Transitions[ipISO][nelem][ipHi][ipLo],WeightNew);
112  if( DEBUG_LOC )
113  {
114  if( ipISO==1 && nelem==1 && ipHi==iso.nLyaLevel[ipISO] && ipLo==0 )
115  fprintf(ioQQQ,"DEBUG rt after loop %li %li %li %li tot %.3e in %.3e\n",
116  ipISO, nelem, ipHi , ipLo ,
117  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauTot ,
118  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->TauIn);
119  }
120  }
121  }
122  /* the extra Lyman lines */
123  for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
124  {
125  /* fully transfer all of the extra lines even though
126  * have not solved for their upper level populations */
127  RT_line_one_tau_reset(&ExtraLymanLines[ipISO][nelem][ipHi],WeightNew);
128  }
129  }
130  }
131  }
132 
133  /* >>>chng 99 nov 11 did not have case b for hydrogenic species on second and
134  * higher iterations */
135  /* option to clobber these taus for Lyman lines, if case b is set */
136  if( opac.lgCaseB )
137  {
138  for( nelem=0; nelem < LIMELM; nelem++ )
139  {
140  if( dense.lgElmtOn[nelem] )
141  {
142  realnum f;
143  /* La may be case B, tlamin set to 1e9 by default with case b command */
144  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn = opac.tlamin;
145  /* >>>chng 99 nov 22, did not reset TauCon */
146  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn;
147  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauTot =
148  2.f*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn;
149  f = opac.tlamin/Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->opacity;
150 
151  for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
152  {
153  if( Transitions[ipH_LIKE][nelem][ipHi][ipH1s].ipCont <= 0 )
154  continue;
155 
156  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn =
157  f*Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->opacity;
158  /* reset line optical depth to continuum source */
159  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauCon = Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauIn;
160  Transitions[ipH_LIKE][nelem][ipHi][ipH1s].Emis->TauTot =
161  2.f*Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->TauIn;
162  }
163  }
164  }
165 
166  /* now do helium like sequence - different since collapsed levels
167  * all go to ground */
168  for( nelem=1; nelem < LIMELM; nelem++ )
169  {
170  if( dense.lgElmtOn[nelem] )
171  {
172  double Aprev;
173  realnum ratio;
174  /* La may be case B, tlamin set to 1e9 by default with case b command */
175  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn = opac.tlamin;
176 
177  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauCon = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
178 
179  Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauTot =
180  2.f*Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->TauIn;
181 
182  ratio = opac.tlamin/Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->opacity;
183 
184  /* this will be the trans prob of the previous lyman line, will use this to
185  * find the next one up in the series */
186  Aprev = Transitions[ipHE_LIKE][nelem][ipHe2p1P][ipHe1s1S].Emis->Aul;
187 
188  i = ipHe2p1P+1;
189  /* >>chng 02 jan 05, remove explicit list of lyman lines, use As to guess
190  * which are which - this will work for any number of levels */
191  for( i = ipHe2p1P+1; i < iso.numLevels_max[ipHE_LIKE][nelem]; i++ )
192  {
193  if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
194  continue;
195 
196  /* >>chng 02 mar 19 use proper test for resonance collapsed lines */
197  if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul> Aprev/10. ||
198  StatesElem[ipHE_LIKE][nelem][i].S < 0 )
199  {
200  Aprev = Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul;
201  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn =
202  ratio*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->opacity;
203  /* reset line optical depth to continuum source */
204  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauCon =
205  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
206  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauTot =
207  2.f*Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->TauIn;
208  /*fprintf(ioQQQ,"%li\t%li\t%.2e\t%.2e\n",nelem, i,
209  Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Aul, Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].TauIn );*/
210  }
211  }
212  }
213  }
214  }
215 
216  /* all heavy element lines */
217  for( i=1; i <= nLevel1; i++ )
218  {
219  RT_line_one_tau_reset(&TauLines[i],WeightNew);
220  /* >>chng 96 jul 06 following sanity check added */
221  ASSERT( fabs(TauLines[i].Emis->TauIn) <= fabs(TauLines[i].Emis->TauTot) );
222  }
223 
224  /* zero out fine opacity fine grid fine mesh array */
225  memset(rfield.fine_opt_depth , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
226 
227  /* all level 2 heavy element lines */
228  for( i=0; i < nWindLine; i++ )
229  {
230  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
231  {
232  RT_line_one_tau_reset(&TauLine2[i],WeightNew);
233  }
234  }
235 
236  /* all hyperfine structure lines */
237  for( i=0; i < nHFLines; i++ )
238  {
239  RT_line_one_tau_reset(&HFLines[i],WeightNew);
240  }
241 
242  /*Atoms & Molecules*/
243  for( i=0; i < linesAdded2; i++)
244  {
245  RT_line_one_tau_reset( atmolEmis[i].tran,WeightNew );
246  }
247 
248  /* inner shell lines */
249  for( i=0; i < nUTA; i++ )
250  {
251  if( UTALines[i].Emis->Aul > 0. )
252  {
253  /* these are line optical depth arrays
254  * inward optical depth */
255  /* heat is special for this array - it is heat per pump */
256  double hsave = UTALines[i].Coll.heat;
257  RT_line_one_tau_reset(&UTALines[i],WeightNew);
258  UTALines[i].Coll.heat = hsave;
259  }
260  }
261 
262  /* co carbon monoxide lines */
263  for( i=0; i < nCORotate; i++ )
264  {
265  RT_line_one_tau_reset(&C12O16Rotate[i],WeightNew);
266  RT_line_one_tau_reset(&C13O16Rotate[i],WeightNew);
267  }
268 
269  /* the large H2 molecule */
270  H2_RT_tau_reset();
271 
272  /* large FeII atom */
274 
275  if( opac.lgCaseB )
276  {
277  for( i=0; i < rfield.nupper; i++ )
278  {
279  /* DEPABS and SCT are abs and sct optical depth for depth only
280  * we will not change total optical depths, just reset inner to half
281  * TauAbsGeo(i,2) = 2.*TauAbsFace(i) */
282  opac.TauAbsGeo[0][i] = opac.TauAbsGeo[1][i]/2.f;
283  /* TauScatGeo(i,2) = 2.*TauScatFace(i) */
284  opac.TauScatGeo[0][i] = opac.TauScatGeo[1][i]/2.f;
285  }
286  }
287  else if( geometry.lgSphere )
288  {
289  /* closed or spherical case */
290  for( i=0; i < rfield.nupper; i++ )
291  {
292  /* [1] is total optical depth from previous iteration,
293  * [0] is optical depth at current position */
294  opac.TauAbsGeo[1][i] = 2.f*opac.TauAbsFace[i];
295  opac.TauAbsGeo[0][i] = opac.TauAbsFace[i];
296  opac.TauScatGeo[1][i] = 2.f*opac.TauScatFace[i];
297  opac.TauScatGeo[0][i] = opac.TauScatFace[i];
298  opac.TauTotalGeo[1][i] = opac.TauScatGeo[1][i] + opac.TauAbsGeo[1][i];
299  opac.TauTotalGeo[0][i] = opac.TauScatGeo[0][i] + opac.TauAbsGeo[0][i];
300  }
301  }
302  else
303  {
304  /* open geometry */
305  for( i=0; i < rfield.nupper; i++ )
306  {
307  opac.TauTotalGeo[1][i] = opac.TauTotalGeo[0][i];
308  opac.TauTotalGeo[0][i] = opac.taumin;
309  opac.TauAbsGeo[1][i] = opac.TauAbsGeo[0][i];
310  opac.TauAbsGeo[0][i] = opac.taumin;
311  opac.TauScatGeo[1][i] = opac.TauScatGeo[0][i];
312  opac.TauScatGeo[0][i] = opac.taumin;
313  }
314  }
315 
316  /* same for open and closed geometries */
317  for( i=0; i < rfield.nupper; i++ )
318  {
319  /* total optical depth across computed shell */
321  /* e2( tau across shell), optical depth from ill face to shielded face of cloud
322  * not that opac.TauAbsFace is reset to small number just after this */
324  /* TauAbsFace and TauScatFace are abs and sct optical depth to ill face */
326  opac.TauAbsFace[i] = opac.taumin;
327  }
328 
329  /* this is optical depth at x-ray point defining effective optical depth */
330  rt.tauxry = opac.TauAbsGeo[0][rt.ipxry-1];
331  return;
332 }

Generated for cloudy by doxygen 1.8.1.2