cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_line_all.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_line_all do escape and destruction probabilities for all lines in code.
4  * Called with false arg by ionize, to only redo destruction probabilities,
5  * and with true by cloudy to do both escape and destruction */
6 #include "cddefines.h"
7 #include "taulines.h"
8 #include "atomfeii.h"
9 #include "dense.h"
10 #include "conv.h"
11 #include "atoms.h"
12 #include "rfield.h"
13 #include "wind.h"
14 #include "iso.h"
15 #include "h2.h"
16 #include "opacity.h"
17 #include "trace.h"
18 #include "lines_service.h"
19 #include "atmdat.h"
20 #include "hydrogenic.h"
21 #include "rt.h"
22 
24  /* this is true if we want to do both escape and destruction probabilities,
25  * and false if only destruction probabilities are needed */
26  bool lgDoEsc ,
27  /* flag saying whether to update fine opacities */
28  bool lgUpdateFineOpac )
29 {
30  long int i,
31  ion,
32  ipISO,
33  nelem;
34  long ipHi , ipLo;
35  double factor,
36  coloi;
37  double SaveLyaPesc[NISO][LIMELM] ,
38  SaveLyaPdest[NISO][LIMELM];
39  /* this flag says whether to update the line RT */
40  bool lgRT_update = lgUpdateFineOpac || lgDoEsc || conv.lgIonStageTrimed;
41 
42  DEBUG_ENTRY( "RT_line_all()" );
43 
44  if( trace.lgTrace )
45  fprintf( ioQQQ, " RT_line_all called\n" );
46 
47  /* skip line transfer if requested, and not first call in this zone */
49  {
50  return;
51  }
52 
53  /* this array is huge and takes significant time to zero out or update,
54  * only do so when needed, */
55  rfield.lgFine_opac_update = lgUpdateFineOpac;
57  {
58  /* zero out fine opacity array */
59  memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
60 
61  /* this is offset within fine opacity array caused by Doppler shift
62  * between first zone, the standard of rest, and current position
63  * in case of acceleration away from star in outflowing wind
64  * dWind is positive, this means that the rest frame original
65  * velocity is red shifted to lower energy */
66  realnum dWind = wind.windv - wind.windv0;
67  /* will add ipVelShift to index of original energy so redshift has
68  * to be negative */
69  rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
70  }
71 
72 # if 0
73  /* this code is a copy of the code within line_one that does cloaking
74  * within this zone. it can be used to see how a particular line
75  * is being treated. */
76  {
77 #include "doppvel.h"
78  double dTau , aa;
79  ipISO = 0; nelem = 0;ipLo = 0;
80  ipHi = 23;
81  dTau = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc *
82  Transitions[ipISO][nelem][ipHi][ipLo].Emis->opacity /
83  DoppVel.doppler[nelem] + opac.opacity_abs[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
85  aa = log(1. + dTau ) / SDIV(dTau);
86  fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
88  aa);
89  }
90 # endif
91 
92  /* find Stark escape probabilities for hydrogen itself */
93  RT_stark();
94 
95  /*if( lgUpdateFineOpac )
96  fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
97  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
98  {
99  /* loop over all iso-electronic sequences */
100  for( nelem=ipISO; nelem<LIMELM; ++nelem )
101  {
102  /* parent ion stage, for H is 1, for He is 1 for He-like and
103  * 2 for H-like */
104  ion = nelem+1-ipISO;
105 
106  /* element turned off */
107  if( !dense.lgElmtOn[nelem] )
108  continue;
109  /* need we consider this ion? */
110  if( ion <=dense.IonHigh[nelem] )
111  {
112  /* save escape and destruction probs for each Lya since
113  * these are special cases, and quantities are not damped
114  * in routine that evaluates each, test is to avoid
115  * evaluating iso sequence that does not exist for an element,
116  * as in He-like hydrogen */
117  if( ipISO<=nelem )
118  {
119  SaveLyaPesc[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc;
120  SaveLyaPdest[ipISO][nelem] = Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest;
121  }
122 
123  /* convert pops to per unit vol rather than per ion */
124  if( dense.xIonDense[nelem][ion] > 1e-30 )
125  {
126  factor = dense.xIonDense[nelem][ion];
127  }
128  else
129  {
130  /* case where almost no parent ion - this will make
131  * very large line opacity, so background dest small */
132  factor = 1.;
133  }
134 
135  /* loop over all lines */
136  for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ++ipHi )
137  {
138  for( ipLo=0; ipLo < ipHi; ++ipLo )
139  {
140  /* negative ipCont means this is not a real line, so do not
141  * transfer it */
142  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 0 )
143  continue;
144 
145  double SavePopOpc = Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc;
146  /* must temporarily make ipLnPopOpc physical */
147  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc *= factor;
148  /*if( ipISO==0 && nelem==0 && ipHi==12 && ipLo==11 )
149  fprintf(ioQQQ,"DEBUG rt loop %li %li %li %li %.3e %.3e %.3e\n",
150  ipISO, nelem, ipHi , ipLo ,
151  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc ,
152  Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn,
153  Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot);*/
154 
155  /* generate escape prob, pumping rate, destruction prob,
156  * inward outward fracs */
157  RT_line_one(&Transitions[ipISO][nelem][ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true);
158 
159 
160  {
161  /* set true to print pump rates*/
162  enum {DEBUG_LOC=false};
163  if( DEBUG_LOC && nelem==1&& ipLo==0 /*&& iteration==2*/ )
164  {
165  fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
166  ipHi ,
167  fnzone,
168  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest);
169  }
170  }
171 
172  /* go back to original units so that final correction ok */
173  Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = SavePopOpc;
174  }
175  }
176 
177  if( ipISO > 0 && lgDoEsc )
178  {
179  /* do not let too much Lya escape outward since so important
180  * H Lya is special case done above */
181  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc *=
182  opac.ExpmTau[Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].ipCont-1];
183  }
184 
185  /* now update two photon induced rates */
186  atmdat_2phot_rate(ipISO , nelem);
187 
188  /* this is option to not do destruction probabilities
189  * for case b no pdest option */
190  if( opac.lgCaseB_no_pdest )
191  {
192  ipLo = 0;
193  /* okay to let this go to numLevels_max. */
194  for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi )
195  {
196  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
197  continue;
198 
199  /* hose the previously computed destruction probability */
200  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest = SMALLFLOAT;
201  }
202  }
203 
204  ipLo = 0;
205  /* these are the extra lyman lines for the iso sequences */
206  /*for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )*/
207  /* only update if significant abundance and need to update fine opac */
208  if( dense.xIonDense[nelem][ion] > 1e-30 && lgUpdateFineOpac )
209  {
210  /* we just want the population of the ground state */
211  factor = StatesElem[ipISO][nelem][0].Pop * dense.xIonDense[nelem][ion];
212  for( ipHi=StatesElem[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1].n; ipHi < iso.nLyman[ipISO]; ipHi++ )
213  {
214  /* must make ipLnPopOpc physical */
215  ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc = factor;
216  ExtraLymanLines[ipISO][nelem][ipHi].Lo->Pop =
217  StatesElem[ipISO][nelem][ipLo].Pop;
218 
219  /* actually do the work */
220  RT_line_one(&ExtraLymanLines[ipISO][nelem][ipHi] , lgDoEsc , lgUpdateFineOpac,true);
221 
222  /* reset to funny population of ground state */
223  ExtraLymanLines[ipISO][nelem][ipHi].Emis->PopOpc =
224  StatesElem[ipISO][nelem][0].Pop;
225  }
226  }
227 
228  }/* if nelem if ion <=dense.IonHigh */
229  }/* loop over nelem */
230  }/* loop over ipISO */
231 
232 
233  /* add on Stark escape probabilities for H itself */
234  for( ipISO=0; ipISO<NISO; ipISO++ )
235  {
236  for( nelem=ipISO; nelem<LIMELM; nelem++ )
237  {
238  if( nelem < 2 || dense.lgElmtOn[nelem] )
239  {
240  /* note that pesc and dest are updated no matter what escprob logic we
241  * specify, and that these are not updated if we have overrun the
242  * optical depth scale. only update here in that case */
243  if( lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) )
244  {
245  if( lgDoEsc )
246  {
247  for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
248  {
249  for( ipLo=0; ipLo < ipHi; ipLo++ )
250  {
251  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
252  continue;
253 
254  /* >>chng 03 jun 07, do not clobber esp prob when line is masing -
255  * this had effect of preventing total escape prob from getting larger than 1 */
256  if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc<1. )
257  {
258  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc = MIN2(1.f,
259  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc+
260  (realnum)iso.pestrk[ipISO][nelem][ipLo][ipHi]);
261  }
262  }
263  }
264  }
265  /* always add on stark broadening term to Lya*/
266  else if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc < 1. )
267  {
268  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc = MIN2(1.f,
269  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc+
270  (realnum)iso.pestrk[ipISO][nelem][0][iso.nLyaLevel[ipISO]]);
271  }
272  }
273  }
274  }
275  }
276 
277  /* note that pesc and dest are updated no matter what escprob logic we
278  * specify, and that these are not updated if we have overrun the
279  * optical depth scale. only update here in that case */
281  {
282  /* do the 8446 problem */
283  atom_oi_calc(&coloi);
284 
286  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul);
287 
288  if( trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE] )
289  {
290  fprintf( ioQQQ, " RT_line_all calls P8446 who found pmph31=%10.2e\n",
291  atoms.pmph31 );
292  }
293 
294  /* add on destruction of hydrogen Lya by FeII
295  * now add in FeII deexcitation via overlap,
296  * but only as loss, not photoionization, source
297  * dstfe2lya is Ly-alpha deexcit by FeII overlap - conversion into FeII em */
298  /* find FeII overlap destruction rate,
299  * this does NOTHING when large FeII atom is turned on,
300  * in this case evaluation is done in call to FeIILevelPops */
302  /*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya ,
303  hydro.HLineWidth);*/
304 
305  /* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */
306  /* >>chng 01 apr 01, add test for tout >= 0.,
307  * must not add to Pdest when it has not been refreshed here */
309  }
310 
311  /* take mean of old and new to damp oscillations in Lya (only) */
312  if( nzone > 1 )
313  {
314  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
315  {
316  /* loop over all iso-electronic sequences */
317  for( nelem=ipISO; nelem<LIMELM; ++nelem )
318  {
319  ion = nelem+1-ipISO;
320  if( ion <=dense.IonHigh[nelem] && ipISO<=nelem )
321  {
322  /* now take mean of old and new escape dest probs for Lya */
323  /* only ipLY_A redist did not already take average */
324  if( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->iRedisFun==ipLY_A &&
325  lgTauGood( &Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0] ) )
326  {
327  /* the Lya routine is different in that escape and dest
328  * probs are always updated, unless we have overrun the
329  * optical depth scale */
330  /*lint -e644 SaveLyaPdest not init */
331  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest =
332  ((realnum)SaveLyaPdest[ipISO][nelem] +
333  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pdest) / 2.f;
334  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc =
335  ((realnum)SaveLyaPesc[ipISO][nelem] +
336  Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Emis->Pesc) / 2.f;
337  /*lint +e644 SaveLyaPdest not init */
338  }
339  }
340  }
341  }
342  }
343  /*5986 fprintf(ioQQQ,"DEBUG he rt_all\t%li\t%li\t%.4e\t%.4e\n",
344  iteration, nzone,
345  Transitions[1][1][10][5].Emis->TauIn,
346  Transitions[1][1][10][5].Emis->Pesc);*/
347 
348  /* is continuum pumping of H Lyman lines included? yes, but turned off
349  * with atom h-like Lyman pumping off command */
350  if( !hydro.lgLymanPumping )
351  {
352  ipISO = ipH_LIKE;
353  nelem = ipHYDROGEN;
354  ipLo = 0;
355  for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipISO][nelem]; ++ipHi )
356  {
357  /* negative ipCont means this is not a real line, so do not
358  * transfer it */
359  Transitions[ipISO][nelem][ipHi][ipLo].Emis->pump = 0.;
360  }
361  }
362 
363  {
364  /* following should be set true to print ots contributors for he-like lines*/
365  /*@-redef@*/
366  enum {DEBUG_LOC=false};
367  /*@+redef@*/
368  if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
369  {
370  /* option to dump a line */
371  DumpLine(&Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0] );
372  }
373  }
374 
375  /* level 1 lines */
376  for( i=1; i <= nLevel1; i++ )
377  {
378  RT_line_one(&TauLines[i] , lgDoEsc , lgUpdateFineOpac,true);
379  }
380  /* co carbon monoxide */
381  for( i=0; i < nCORotate; i++ )
382  {
383  RT_line_one(&C12O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true);
384  RT_line_one(&C13O16Rotate[i] , lgDoEsc , lgUpdateFineOpac,true);
385  }
386  for( i=0; i < nHFLines; i++ )
387  {
388  RT_line_one(&HFLines[i] , lgDoEsc , lgUpdateFineOpac,true);
389  }
390  /* Database Atoms & Molecules*/
391  for( i=0; i < linesAdded2; i++)
392  {
393  RT_line_one(atmolEmis[i].tran , lgDoEsc , lgUpdateFineOpac,true);
394  }
395  /* this is a major time sink for this routine */
396  if( lgRT_update )
397  {
398  for( i=0; i < nUTA; i++ )
399  {
400  /* Aul = -1 means this line is superceded by better calculations in
401  * different data set */
402  if( UTALines[i].Emis->Aul > 0. )
403  {
404  /* these are not defined in cooling routines so must be set here */
406  UTALines[i].Lo->Pop = dense.xIonDense[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1];
407  UTALines[i].Hi->Pop = 0.;
408  RT_line_one(&UTALines[i] , lgDoEsc , lgUpdateFineOpac,true);
409  }
410  }
411  }
412 
413  /* h-like ions only have one electron, no satellite lines. */
414  for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
415  {
416  /* loop over all iso-electronic sequences */
417  for( nelem=ipISO; nelem<LIMELM; ++nelem )
418  {
419  if( dense.lgElmtOn[nelem] )
420  {
421  for( ipLo=1; ipLo < (iso.numLevels_local[ipISO][nelem]-1); ++ipLo )
422  {
423  RT_line_one(&SatelliteLines[ipISO][nelem][ipLo], lgDoEsc, lgUpdateFineOpac, true);
424  }
425  }
426  }
427  }
428 
429  /* >>chng 03 aug 22, must call RT_line_one every single time, since that routine
430  * establishes the fine opacities for the level 2 lines */
431  /* >>chng 01 aug 11, the level 2 lines were only updated when lgDoEsc was true,
432  * this meant that dest probs were not updated but once, causing very high-Z models
433  * to develop numerical oscialltions. removed if, now evaluate level 2 always */
434  /* only update their dest/esc probs one time per zone */
435  /* lgLevel2_OTS_Imp set true in dimacool if ots rates were significant */
436  /* level 2 heavy element lines in cooling with only g-bar,*/
437  for( i=0; i < nWindLine; i++ )
438  {
439  /* >>chng 02 sug 11, do not include ions done in iso-sequences */
440  /*if( TauLine2[i].IonStg <= TauLine2[i].nelem-1 )*/
441  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
442  {
443  RT_line_one(&TauLine2[i] , lgDoEsc , lgUpdateFineOpac,true);
444  }
445  }
446 
447  /* the large H2 molecule */
448  H2_RTMake( lgDoEsc , lgUpdateFineOpac);
449 
450  /* The large model FeII atom */
451  FeII_RT_Make( lgDoEsc , lgUpdateFineOpac);
452  return;
453 }

Generated for cloudy by doxygen 1.8.3.1