cloudy  trunk
 All Data Structures Namespaces 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-2017 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 #include "cddefines.h"
5 #include "taulines.h"
6 #include "dense.h"
7 #include "conv.h"
8 #include "rfield.h"
9 #include "wind.h"
10 #include "iso.h"
11 #include "h2.h"
12 #include "opacity.h"
13 #include "trace.h"
14 #include "hydrogenic.h"
15 #include "rt.h"
16 #include "cosmology.h"
17 #include "doppvel.h"
18 #include "atmdat.h"
19 
20 /*RT_line_all do escape and destruction probabilities for all lines in code. */
22 {
23  DEBUG_ENTRY( "RT_line_all_escape()" );
24  /* this flag says whether to update the line escape probabilities */
25  bool lgPescUpdate = conv.lgFirstSweepThisZone || conv.lgIonStageTrimed;
26  /* find Stark escape probabilities for hydrogen itself */
27  if( lgPescUpdate )
28  RT_stark();
29 
30  bool lgCheck = (error != NULL);
31  vector<realnum> oldPesc, oldPdest, oldPelec_esc;
32  if (lgCheck)
33  {
34  for( vector<TransitionList>::iterator it = AllTransitions.begin(); it != AllTransitions.end(); ++it )
35  {
36  for (TransitionList::iterator tr=it->begin();
37  tr != it->end(); ++tr)
38  {
39  if ((*tr).ipCont() <= 0)
40  continue;
41  oldPdest.push_back((*tr).Emis().Pdest());
42  oldPesc.push_back((*tr).Emis().Pesc());
43  oldPelec_esc.push_back((*tr).Emis().Pelec_esc());
44  }
45  }
46  }
47 
48 
50 
52  {
53  for( long ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
54  {
55  /* loop over all iso-electronic sequences */
56  for( long nelem=ipISO; nelem < LIMELM; ++nelem )
57  {
58  /* parent ion stage, for H is 1, for He is 1 for He-like and
59  * 2 for H-like */
60  long ion = nelem+1-ipISO;
61 
62  /* element turned off */
63  if( !dense.lgElmtOn[nelem] )
64  continue;
65  /* need we consider this ion? */
66  if( ion <= dense.IonHigh[nelem] )
67  {
68  /* loop over all lines */
69  /* this is option to not do destruction probabilities
70  * for case b no pdest option */
71  long ipLo = 0;
72  /* okay to let this go to numLevels_max. */
73  for( long ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
74  {
75  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
76  continue;
77 
78  /* hose the previously computed destruction probability */
79  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest() = SMALLFLOAT;
80  }
81  }
82  }
83  }
84  }
85 
86  /* note that pesc and dest are updated no matter what escprob logic we
87  * specify, and that these are not updated if we have overrun the
88  * optical depth scale. only update here in that case */
89  // Logic must be kept consistent with mask for RT_line_escape in rt_line_one.cpp
92  {
93  /*fprintf(ioQQQ,"DEBUG fe2 %.2e %.2e\n", hydro.dstfe2lya ,
94  hydro.HLineWidth);*/
95 
96  /* >>chng 00 jan 06, let dest be large than 1 to desaturate the atom */
97  /* >>chng 01 apr 01, add test for tout >= 0.,
98  * must not add to Pdest when it has not been refreshed here */
100  }
101 
102  /* is continuum pumping of H Lyman lines included? yes, but turned off
103  * with atom h-like Lyman pumping off command */
104  if( !hydro.lgLymanPumping )
105  {
106  long ipISO = ipH_LIKE;
107  long nelem = ipHYDROGEN;
108  long ipLo = 0;
109  for( long ipHi=ipLo+1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ++ipHi )
110  {
111  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump() = 0.;
112  }
113  }
114 
115  {
116  /* following should be set true to print ots contributors for he-like lines*/
117  enum {DEBUG_LOC=false};
118  if( DEBUG_LOC && nzone>433 /*&& iteration==2*/ )
119  {
120  /* option to dump a line */
121  DumpLine(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0) );
122  }
123  }
124 # if 0
125  /* this code is a copy of the code within line_one that does cloaking
126  * within this zone. it can be used to see how a particular line
127  * is being treated. */
128  {
129 #include "doppvel.h"
130  double dTau , aa;
131  ipISO = 0; nelem = 0;ipLo = 0;
132  ipHi = 23;
133  dTau = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() *
134  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().opacity() /
135  GetDopplerWidth(dense.AtomicWeight[nelem]) + opac.opacity_abs[iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont()-1];
136  dTau *= radius.drad_x_fillfac;
137  aa = log(1. + dTau ) / SDIV(dTau);
138  fprintf(ioQQQ,"DEBUG dTau\t%.2f\t%.5e\t%.5e\t%.5e\n",fnzone,dTau,
140  aa);
141  }
142 # endif
143 
144  realnum bigerror=0.0;
145  TransitionList::iterator trworst;
146  long iworst=-1,cworst=-1;
147  if (lgCheck)
148  {
149  long ind=0;
150  for( vector<TransitionList>::iterator it = AllTransitions.begin(); it != AllTransitions.end(); ++it )
151  {
152  for (TransitionList::iterator tr=it->begin();
153  tr != it->end(); ++tr)
154  {
155  if ((*tr).ipCont() <= 0)
156  continue;
157 
158  const realnum abslim = 1e-3;
159 
160  realnum frac = abs(oldPdest[ind]-(*tr).Emis().Pdest())/max((*tr).Emis().Pdest(),abslim);
161  if (frac > bigerror)
162  {
163  bigerror = frac;
164  trworst = tr;
165  iworst = ind;
166  cworst = 1;
167  }
168  frac = abs(oldPesc[ind]-(*tr).Emis().Pesc())/max((*tr).Emis().Pesc(),abslim);
169  if (frac > bigerror)
170  {
171  bigerror = frac;
172  trworst = tr;
173  iworst = ind;
174  cworst = 2;
175  }
176  frac = abs(oldPelec_esc[ind]-(*tr).Emis().Pelec_esc())/max((*tr).Emis().Pelec_esc(),abslim);
177  if (frac > bigerror)
178  {
179  bigerror = frac;
180  trworst = tr;
181  iworst = ind;
182  cworst = 3;
183  }
184  ++ind;
185  }
186  }
187  if (0 && bigerror > 0.0)
188  {
189  fprintf(ioQQQ,"RTEscChange nz %4ld %6f ",nzone,bigerror);
190  fprintf(ioQQQ,"dst%c %6f=>%6f esc%c %6f=>%6f eesc%c %6f=>%6f -- %s\n",
191  cworst == 1 ? '*':' ',oldPdest[iworst],(*trworst).Emis().Pdest(),
192  cworst == 2 ? '*':' ',oldPesc[iworst],(*trworst).Emis().Pesc(),
193  cworst == 3 ? '*':' ',oldPelec_esc[iworst],(*trworst).Emis().Pelec_esc(),
194  (*trworst).chLabel().c_str());
195  }
196  *error = bigerror;
197  }
198 }
199 
200 void RT_line_all( linefunc line_one )
201 {
202  long int ion,
203  ipISO,
204  nelem;
205  long ipHi , ipLo;
206 
207  DEBUG_ENTRY( "RT_line_all()" );
208 
209  if( trace.lgTrace )
210  fprintf( ioQQQ, " RT_line_all called\n" );
211 
212  /* rfield.lgDoLineTrans - skip line transfer if requested with no line transfer
213  * conv.nPres2Ioniz is zero during search phase and on first call in this zone */
215  {
216  return;
217  }
218 
219  /*if( lgUpdateFineOpac )
220  fprintf(ioQQQ,"debuggg\tlgUpdateFineOpac in rt_line_all\n");*/
221 
222  static vector<realnum> DopplerWidth(LIMELM);
223  for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
224  {
225  DopplerWidth[nelem] = GetDopplerWidth(dense.AtomicWeight[nelem]);
226  }
227 
228  for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
229  {
230  /* loop over all iso-electronic sequences */
231  for( nelem=ipISO; nelem < LIMELM; ++nelem )
232  {
233  /* parent ion stage, for H is 1, for He is 1 for He-like and
234  * 2 for H-like */
235  ion = nelem+1-ipISO;
236 
237  /* element turned off */
238  if( !dense.lgElmtOn[nelem] )
239  continue;
240  /* need we consider this ion? */
241  if( ion <= dense.IonHigh[nelem] )
242  {
243  /* loop over all lines */
244  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ++ipHi )
245  {
246  for( ipLo=0; ipLo < ipHi; ++ipLo )
247  {
248  /* negative ipCont means this is not a real line, so do not
249  * transfer it */
250  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() < 0 )
251  continue;
252 
253  /* generate escape prob, pumping rate, destruction prob,
254  * inward outward fracs */
255  fixit("should this use pestrk_up or pestrk?");
256  line_one( iso_sp[ipISO][nelem].trans(ipHi,ipLo),
257  true,(realnum)iso_sp[ipISO][nelem].ex[ipHi][ipLo].pestrk_up,
258  DopplerWidth[nelem]);
259 
260  /* set true to print pump rates*/
261  enum {DEBUG_LOC=false};
262  if( DEBUG_LOC && nelem==1&& ipLo==0 /*&& iteration==2*/ )
263  {
264  fprintf(ioQQQ,"DEBUG pdest\t%3li\t%.2f\t%.3e\n",
265  ipHi ,
266  fnzone,
267  iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pdest());
268  }
269  }
270  }
271  }
272  }
273  }
274 
276  {
277  for( ipISO=ipH_LIKE; ipISO < NISO; ++ipISO )
278  {
279  /* loop over all iso-electronic sequences */
280  for( nelem=ipISO; nelem < LIMELM; ++nelem )
281  {
282  /* parent ion stage, for H is 1, for He is 1 for He-like and
283  * 2 for H-like */
284  ion = nelem+1-ipISO;
285 
286  /* element turned off */
287  if( !dense.lgElmtOn[nelem] )
288  continue;
289  /* need we consider this ion? */
290  if( ion <= dense.IonHigh[nelem] )
291  {
292  /* loop over all lines */
293  ipLo = 0;
294  /* these are the extra Lyman lines for the iso sequences */
295  /*for( ipHi=2; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
296  /* only update if significant abundance and need to update fine opac */
297  if( dense.xIonDense[nelem][ion] > 1e-30 )
298  {
299  for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
300  {
301  TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
302  /* we just want the population of the ground state */
303  (*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();
304  (*(*tr).Lo()).Pop() =
305  iso_sp[ipISO][nelem].st[ipLo].Pop();
306 
307  /* actually do the work */
308  line_one( *tr, true, 0.f, DopplerWidth[nelem]);
309  }
310  }
311  }/* if nelem if ion <=dense.IonHigh */
312  }/* loop over nelem */
313  }/* loop over ipISO */
314 
315  /* this is a major time sink for this routine - only evaluate on last
316  * sweep when fine opacities are updated since only effect of UTAs is
317  * to pump inner shell lines and add to total opacity */
318 
319  for( size_t i=0; i < UTALines.size(); i++ )
320  {
321  /* these are not defined in cooling routines so must be set here */
322  UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
323  (*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
324  (*UTALines[i].Hi()).Pop() = 0.;
325  line_one( UTALines[i], true,0.f,
326  DopplerWidth[(*UTALines[i].Hi()).nelem()-1] );
327  }
328  }
329 
330 
331  /* external database lines */
332  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
333  {
334  if( dBaseSpecies[ipSpecies].lgActive )
335  {
336  realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
337  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
338  tr != dBaseTrans[ipSpecies].end(); ++tr)
339  {
340  int ipHi = (*tr).ipHi();
341  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
342  continue;
343  if( (*tr).ipCont() > 0 )
344  {
345  line_one( *tr, true, 0.f, DopplerWidth );
346  }
347  }
348  }
349  }
350 
351  /* do satellite lines for iso sequences gt H-like
352  * H-like ions only have one electron, no satellite lines. */
353  for( ipISO=ipHE_LIKE; ipISO<NISO; ++ipISO )
354  {
355  /* loop over all iso-electronic sequences */
356  for( nelem=ipISO; nelem<LIMELM; ++nelem )
357  {
358  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
359  {
360  for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ++ipLo )
361  {
362  line_one( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], true,0.f,
363  DopplerWidth[nelem] );
364  }
365  }
366  }
367  }
368 
369  /* the large H2 molecule */
370  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
371  (*diatom)->H2_RTMake( line_one );
372 
373  for( size_t i=0; i < HFLines.size(); i++ )
374  {
375  line_one( HFLines[i], true,0.f, DopplerWidth[(*HFLines[i].Hi()).nelem()-1] );
376  }
377 
378  for( long i=0; i < nWindLine; i++ )
379  {
380  ion = (*TauLine2[i].Hi()).IonStg();
381  nelem = (*TauLine2[i].Hi()).nelem();
382 
383  if( (dense.lgIonChiantiOn[nelem-1][ion-1] && !atmdat.lgChiantiHybrid) ||
384  (dense.lgIonStoutOn[nelem-1][ion-1] && !atmdat.lgStoutHybrid) )
385  {
386  /* If a species uses Chianti or Stout and hybrid is off, skip the level 2 lines */
387  continue;
388  }
389  // iso sequence ions are done elsewhere
390  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO &&
391  // dense.maxWN is positive if hybrid chianti is turned on and this element is included
392  // in CloudyChianti.ini - zero otherwise
393  TauLine2[i].EnergyWN() >= dense.maxWN[nelem-1][ion-1])
394  {
395  line_one( TauLine2[i], true,0.f, DopplerWidth[(*TauLine2[i].Hi()).nelem()-1] );
396  }
397  }
398 
399  return;
400 }
401 
403 {
404  DEBUG_ENTRY( "RT_fine_init()" );
405  /* zero out fine opacity array */
406  /* this array is huge and takes significant time to zero out or update,
407  * only do so when needed, */
408  memset(rfield.fine_opac_zone , 0 , (unsigned long)rfield.nfine*sizeof(realnum) );
409 
410  if( 0 && cosmology.lgDo )
411  {
412  realnum dVel = (realnum)SPEEDLIGHT * ( 1.f - 1.f/POW2(1.f+cosmology.redshift_start-cosmology.redshift_current) );
413  rfield.ipFineConVelShift = -(long int)( dVel/ rfield.fine_opac_velocity_width + 0.5 );
414  }
415  else
416  {
417  /* this is offset within fine opacity array caused by Doppler shift
418  * between first zone, the standard of rest, and current position
419  * in case of acceleration away from star in outflowing wind
420  * dWind is positive, this means that the rest frame original
421  * velocity is red shifted to lower energy */
422  realnum dWind = wind.windv - wind.windv0;
423 
424  /* will add ipVelShift to index of original energy so redshift has
425  * to be negative */
426  rfield.ipFineConVelShift = -(long int)(dWind / rfield.fine_opac_velocity_width+0.5);
427  }
428 }
void DumpLine(const TransitionProxy &t)
Definition: transition.cpp:138
t_atmdat atmdat
Definition: atmdat.cpp:6
realnum & opacity() const
Definition: emission.h:638
bool lgStoutHybrid
Definition: atmdat.h:404
size_t size(void) const
Definition: transition.h:331
void RT_fine_clear()
double * opacity_abs
Definition: opacity.h:103
qList st
Definition: iso.h:482
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
const realnum SMALLFLOAT
Definition: cpu.h:246
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
bool lgFirstSweepThisZone
Definition: conv.h:148
const int NISO
Definition: cddefines.h:310
void RT_line_one_escape(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
long int IonHigh[LIMELM+1]
Definition: dense.h:130
realnum windv0
Definition: wind.h:11
realnum fine_opac_velocity_width
Definition: rfield.h:369
realnum redshift_start
Definition: cosmology.h:26
t_conv conv
Definition: conv.cpp:5
TransitionList HFLines("HFLines",&AnonStates)
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nzone
Definition: cddefines.cpp:14
TransitionList TauLine2("TauLine2",&AnonStates)
bool lgDo
Definition: cosmology.h:44
void RT_line_all_escape(realnum *error)
Definition: rt_line_all.cpp:21
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
bool lgIonStageTrimed
Definition: conv.h:184
bool lgIonStoutOn[LIMELM][LIMELM+1]
Definition: dense.h:143
bool lgCaseB_no_pdest
Definition: opacity.h:184
Wind wind
Definition: wind.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_trace trace
Definition: trace.cpp:5
bool lgIonChiantiOn[LIMELM][LIMELM+1]
Definition: dense.h:140
bool lgChiantiHybrid
Definition: atmdat.h:376
realnum redshift_current
Definition: cosmology.h:26
#define POW2
Definition: cddefines.h:983
long int nLyman[NISO]
Definition: iso.h:352
const int ipH1s
Definition: iso.h:29
long int nPres2Ioniz
Definition: conv.h:145
bool lgTrace
Definition: trace.h:12
EmissionList::reference Emis() const
Definition: transition.h:447
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
void(* linefunc)(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
Definition: h2_priv.h:10
t_rfield rfield
Definition: rfield.cpp:9
long & ipCont() const
Definition: transition.h:489
bool lgDielRecom[NISO]
Definition: iso.h:385
float realnum
Definition: cddefines.h:124
long int ipFineConVelShift
Definition: rfield.h:399
vector< diatomics * > diatoms
Definition: h2.cpp:8
long max(int a, long b)
Definition: cddefines.h:821
t_hydro hydro
Definition: hydrogenic.cpp:5
void RT_stark(void)
Definition: rt_stark.cpp:12
void RT_line_all(linefunc line_one)
realnum GetDopplerWidth(realnum massAMU)
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:35
long nWindLine
Definition: cdinit.cpp:19
t_radius radius
Definition: radius.cpp:5
bool lgDoLineTrans
Definition: rfield.h:101
bool lgElmtOn[LIMELM]
Definition: dense.h:160
realnum * fine_opac_zone
Definition: rfield.h:389
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
int nLyaLevel[NISO]
Definition: iso.h:397
const int ipH2p
Definition: iso.h:31
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
realnum & Pdest() const
Definition: emission.h:588
bool lgLymanPumping
Definition: hydrogenic.h:144
const int ipH_LIKE
Definition: iso.h:64
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
const int LIMELM
Definition: cddefines.h:307
long nfine
Definition: rfield.h:385
double drad_x_fillfac
Definition: radius.h:76
realnum dstfe2lya
Definition: hydrogenic.h:96
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double & PopOpc() const
Definition: emission.h:658
t_cosmology cosmology
Definition: cosmology.cpp:8
vector< species > dBaseSpecies
Definition: taulines.cpp:15
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
double maxWN[LIMELM][LIMELM+1]
Definition: dense.h:146
vector< TransitionList > AllTransitions
Definition: taulines.cpp:9
long int numLevels_max
Definition: iso.h:524
bool lgLastSweepThisZone
Definition: conv.h:150
#define fixit(a)
Definition: cddefines.h:416
bool lgTauGood(const TransitionProxy &t)
Definition: transition.h:644
double fnzone
Definition: cddefines.cpp:15
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
double frac(double d)
const int ipHYDROGEN
Definition: cddefines.h:348
long int numLevels_local
Definition: iso.h:529
realnum windv
Definition: wind.h:18
EmissionList & Emis()
Definition: transition.h:363
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
double & pump() const
Definition: emission.h:518