cloudy  trunk
rt_diffuse.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_diffuse evaluate local diffuse emission for this zone,
4  * fill in ConEmitLocal[depth][energy] with diffuse emission,
5  * called by Cloudy, this routine adds energy to the outward beam
6  * OTS rates for this zone were set in RT_OTS - not here */
7 #include "cddefines.h"
8 #include "physconst.h"
9 #include "taulines.h"
10 #include "grains.h"
11 #include "grainvar.h"
12 #include "iso.h"
13 #include "dense.h"
14 #include "opacity.h"
15 #include "trace.h"
16 #include "coolheavy.h"
17 #include "rfield.h"
18 #include "phycon.h"
19 #include "hmi.h"
20 #include "radius.h"
21 #include "atmdat.h"
22 #include "heavy.h"
23 #include "atomfeii.h"
24 #include "lines_service.h"
25 #include "h2.h"
26 #include "ipoint.h"
27 #include "rt.h"
28 #include "mole.h"
29 #include "conv.h"
30 
31 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
32 #pragma optimization_level 0
33 #endif
34 void RT_diffuse(void)
35 {
36  /* arrays used in this routine
37  * rfield.ConEmitLocal[depth][energy] local emission per unit vol
38  * rfield.DiffuseEscape is the spectrum of diffuse emission that escapes this zone,
39  * at end of this routine part is thrown into the outward beam
40  * by adding to rfield.ConInterOut
41  * units are photons s-1 cm-3
42  * one-time init done on first call */
43 
44  /* rfield.DiffuseEscape and rfield.ConEmitLocal are same except that
45  * rfield.ConEmitLocal is local emission, would be source function if div by opac
46  * rfield.DiffuseEscape is part that escapes so has RT built into it
47  * rfield.DiffuseEscape is used to define rfield.ConInterOut below as per this statement
48  * rfield.ConInterOut[nu] += rfield.DiffuseEscape[nu]*(realnum)radius.dVolOutwrd;
49  */
50  /* \todo 0 define only rfield.ConEmitLocal as it is now done,
51  * do not define rfield.DiffuseEscape at all
52  * at bottom of this routine use inward and outward optical depths to define
53  * local and escaping parts
54  * this routine only defines
55  * rfield.ConInterOut - set to rfield.DiffuseEscape times vol element
56  * so this is only var that
57  * needs to be set
58  */
59 
60  long int ip=-100000,
61  ipla=-100000,
62  limit=-100000,
63  nu=-10000;
64 
65  double EdenAbund,
66  difflya,
67  fac,
68  factor,
69  gamma,
70  gion,
71  gn,
72  photon;
73 
74  DEBUG_ENTRY( "RT_diffuse()" );
75 
76  /* many arrays were malloced to nupper, and we will add unit flux to [nflux] -
77  8 this must be true to work */
78  ASSERT(rfield.nflux < rfield.nupper );
79 
80  /* this routine evaluates the local diffuse fields
81  * it fills in all of the following vectors */
82  memset(rfield.DiffuseEscape , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
83  memset(rfield.ConEmitLocal[nzone] , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
84  memset(rfield.TotDiff2Pht , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
85  memset(rfield.DiffuseLineEmission , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
86 
87  /* must abort after setting all of above to zero because some may be
88  * used in various ways before abort is complete */
89  if( lgAbort )
90  {
91  /* quit if we are aborting */
92  return;
93  }
94 
95  /* loop over iso-sequences of all elements
96  * to add all recombination continua and lines*/
97  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
98  {
99  /* >>chng 01 sep 23, rewrote for iso sequences */
100  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
101  {
102  // calculate recombination spectra and cooling
103  RT_iso_integrate_RRC( ipISO, nelem, true );
104 
105  /* the product of the densities of the parent ion and electrons */
106  EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
107 
108  /* recombination continua for all iso seq -
109  * if this stage of ionization exists */
110  if( dense.IonHigh[nelem] >= nelem+1-ipISO )
111  {
112  t_iso_sp* sp = &iso_sp[ipISO][nelem];
113 
114  // add line emission from the model iso atoms
115  for( long ipHi=1; ipHi < sp->numLevels_local; ipHi++ )
116  {
117  for( long ipLo=0; ipLo < ipHi; ipLo++ )
118  {
119  // skip non-radiative transitions
120  if( sp->trans(ipHi,ipLo).ipCont() < 1 )
121  continue;
122 
123  /* number of photons in the line has not been defined up until now,
124  * do so now. this is redone in lines. */
125  sp->trans(ipHi,ipLo).Emis().phots() =
126  sp->trans(ipHi,ipLo).Emis().Aul()*
127  sp->st[ipHi].Pop()*
128  sp->trans(ipHi,ipLo).Emis().Pesc();
129 
130  // Would be better to enable checks (and remove argument) --
131  // present state is to ensure backwards compatibility with previous
132  // unchecked code.
133  // First argument is fraction of line not emitted by scattering --
134  // would be better to do this on the basis of line physics rather than
135  // fiat...
136  const bool lgDoChecks = false;
137  sp->trans(ipHi,ipLo).outline(1.0, lgDoChecks );
138  }
139  }
140 
141  /*Iso treatment of two photon emission. */
142  /* NISO could in the future be increased, but we want this assert to blow
143  * so that it is understood this may not be correct for other iso sequences,
144  * probably should break since will not be present */
145  ASSERT( ipISO <= ipHE_LIKE );
146 
147  /* upper limit to 2-phot is energy of 2s to ground */
148  for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
149  {
150  CalcTwoPhotonEmission( *tnu, rfield.lgInducProcess && iso_ctrl.lgInd2nu_On );
151 
152  for( nu=0; nu < tnu->ipTwoPhoE; nu++ )
153  {
154  /* information - only used in save output */
155  rfield.TotDiff2Pht[nu] += tnu->local_emis[nu];
156 
157  /* total local diffuse emission */
158  rfield.ConEmitLocal[nzone][nu] += tnu->local_emis[nu];
159 
160  /* this is escaping part of two-photon emission,
161  * as determined from optical depth to illuminated face */
162  rfield.DiffuseEscape[nu] += tnu->local_emis[nu] * opac.ExpmTau[nu];
163  }
164  enum {DEBUG_LOC=false};
165  if( DEBUG_LOC )
166  {
167  fprintf( ioQQQ, "Two-photon emission coefficients - ipISO, nelem = %2li, %2li\n", ipISO, nelem );
168  PrtTwoPhotonEmissCoef( *tnu, EdenAbund );
169  }
170  }
171  }
172  }
173  }
174 
175  /* add recombination continua for elements heavier than those done with iso seq */
176  for( long nelem=NISO; nelem < LIMELM; nelem++ )
177  {
178  // zero out all stages since dense.IonLow[nelem] may have been lower last time around
179  for( long ion=0; ion < nelem-NISO+1; ion++ )
180  {
181  Heavy.RadRecCon[nelem][ion] = 0.;
182  }
183 
184  /* do not include species with iso-sequence in following */
185  /* >>chng 03 sep 09, upper bound was wrong, did not include NISO */
186  for( long ion=dense.IonLow[nelem]; ion < nelem-NISO+1; ion++ )
187  {
188  if( dense.xIonDense[nelem][ion+1] > 0. )
189  {
190  long int ns, nshell,igRec , igIon,
191  iplow , iphi , ipop;
192 
193  ip = Heavy.ipHeavy[nelem][ion]-1;
194  ASSERT( ip >= 0 );
195 
196  /* nflux was reset upward in ConvInitSolution to encompass all
197  * possible line and continuum emission. this test should not
198  * possibly fail. It could if the ionization were to increase with depth
199  * although the continuum mesh is designed to deal with this.
200  * This test is important because the nflux cell in ConInterOut
201  * is used to carry out the unit integration, and if it gets
202  * clobbered by diffuse emission the code will declare
203  * insanity in PrtComment */
204  if( ip >= rfield.nflux )
205  continue;
206 
207  /* get shell number, stat weights for this species */
208  atmdat_outer_shell( nelem+1 , nelem+1-ion , &nshell, &igRec , &igIon );
209  gn = (double)igRec;
210  gion = (double)igIon;
211 
212  /* shell number */
213  ns = Heavy.nsShells[nelem][ion]-1;
214  ASSERT( ns == (nshell-1) );
215 
216  /* lower and upper energies, and offset for opacity stack */
217  iplow = opac.ipElement[nelem][ion][ns][0]-1;
218  iphi = opac.ipElement[nelem][ion][ns][1];
219  iphi = MIN2( iphi , rfield.nflux );
220  ipop = opac.ipElement[nelem][ion][ns][2];
221 
222  /* now convert ipop to the offset in the opacity stack from threshold */
223  ipop = ipop - iplow;
224 
225  EdenAbund = dense.eden*dense.xIonDense[nelem][ion+1];
226  gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
227 
228  /* this is ground state continuum from stored opacities */
229  if( rfield.ContBoltz[iplow] > SMALLFLOAT )
230  {
231  for( nu=iplow; nu < iphi; ++nu )
232  {
233  photon = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[iplow]*
234  rfield.widflx[nu]*opac.OpacStack[nu+ipop]*rfield.anu2[nu];
235  /* add heavy rec to ground in active beam,*/
240  rfield.ConEmitLocal[nzone][nu] += (realnum)photon*EdenAbund;
241  rfield.DiffuseEscape[nu] += (realnum)photon*EdenAbund*opac.ExpmTau[nu];
242 
243  // escaping RRC
244  Heavy.RadRecCon[nelem][ion] += rfield.anu[nu] *
245  emergent_line( photon*EdenAbund/2. , photon*EdenAbund/2. ,
246  // energy on fortran scale
247  nu+1 );
248  }
249  }
250  // units erg cm-3 s-1
251  Heavy.RadRecCon[nelem][ion] *= EN1RYD;
252 
253  /* now do the recombination Lya */
254  ipla = Heavy.ipLyHeavy[nelem][ion]-1;
255  ASSERT( ipla >= 0 );
256  /* xLyaHeavy is set to a fraction of the total rad rec in ion_recomb, includes eden */
257  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
258  rfield.DiffuseLineEmission[ipla] += (realnum)difflya;
259 
260  /* >>chng 03 jul 10, here and below, use outlin_noplot */
261  rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
262 
263  /* now do the recombination Balmer photons */
264  ipla = Heavy.ipBalHeavy[nelem][ion]-1;
265  ASSERT( ipla >= 0 );
266  /* xLyaHeavy is set to fraction of total rad rec in ion_recomb, includes eden */
267  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
268  rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
269  }
270  }
271  }
272 
273  /* free-free free free brems emission for all ions */
274  limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
275  for( nu=0; nu < limit; nu++ )
276  {
277  double TotBremsAllIons = 0., BremsThisIon;
278 
279  /* First add H- brems. Reaction is H(1s) + e -> H(1s) + e + hnu.
280  * OpacStack contains the ratio of the H- to H brems cross section.
281  * Multiply H brems by this and the population of H(1s). */
282  TotBremsAllIons += rfield.gff[1][nu] * opac.OpacStack[nu-1+opac.iphmra] * iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
283 
284  /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
285  * using gaunt factors from rfield.gff. */
286  for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
287  {
288  /* MAX2 occurs because we want to start at first ion (or above)
289  * and do not want atom */
290  for( long ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
291  {
292  /* eff. charge is ion, so first rfield.gff argument must be "ion". */
293  BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
294  TotBremsAllIons += BremsThisIon;
295  }
296  }
297 
298  /* add molecular ions */
299  for( long ipMol = 0; ipMol<mole_global.num_calc; ipMol++ )
300  {
301  if( !mole_global.list[ipMol]->isMonatomic() && mole_global.list[ipMol]->charge > 0 && mole_global.list[ipMol]->parentLabel.empty()
302  // H2+ and H3+ do not appear to be included above.
303  /* && mole_global.list[ipMol] != findspecies("H2+") &&
304  mole_global.list[ipMol] != findspecies("H3+") */ )
305  {
306  /* eff. charge is ion, so first rfield.gff argument must be "ion". */
307  long ion = mole_global.list[ipMol]->charge;
308  BremsThisIon = POW2( (double)ion )*mole.species[ipMol].den*rfield.gff[ion][nu];
309  TotBremsAllIons += BremsThisIon;
310  }
311  }
312 
314  /* >>chng 06 apr 05, no free free also turns off emission */
315  TotBremsAllIons *= dense.eden*1.032e-11*rfield.widflx[nu]*rfield.ContBoltz[nu]/rfield.anu[nu]/phycon.sqrte *
316  CoolHeavy.lgFreeOn;
317  ASSERT( TotBremsAllIons >= 0.);
318 
319  /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
320  * to outward beam - ConLocNoInter array removed as result
321  * if problems develop with very dense blr clouds, this may be reason */
322  /*rfield.ConLocNoInter[nu] += (realnum)fac;*/
323  /*rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;*/
324 
325  /* >>chng 05 feb 20, move into this test on brems opacity - should not be
326  * needed since would use expmtau to limit outward beam */
327  /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
328  * to outward beam - ConLocNoInter array removed as result
329  * if problems develop with very dense BLR clouds, this may be reason */
330  /*rfield.ConLocNoInter[nu] += (realnum)fac;*/
331  rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;
332 
333  /* do not add optically thick part to outward beam since self absorbed
334  * >>chng 96 feb 27, put back into outward beam since do not integrate
335  * over it anyway. */
336  /* >>chng 99 may 28, take back out of beam since DO integrate over it
337  * in very dense BLR clouds */
338  /* >>chng 01 jul 10, add here, in only one loop, where optically thin */
339  rfield.DiffuseEscape[nu] += (realnum)TotBremsAllIons;
340  }
341 
342  /* grain dust emission */
343  /* >>chng 01 nov 22, moved calculation of grain flux to qheat.c, PvH */
344  if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
345  {
346  /* this calculates diffuse emission from grains,
347  * and stores the result in gv.GrainEmission */
348  GrainMakeDiffuse();
349 
350  for( nu=0; nu < rfield.nflux; nu++ )
351  {
352  rfield.ConEmitLocal[nzone][nu] += gv.GrainEmission[nu];
353  rfield.DiffuseEscape[nu] += gv.GrainEmission[nu];
354  }
355  }
356 
357  /* hminus emission */
358  fac = dense.eden*(double)dense.xIonDense[ipHYDROGEN][0];
359  gn = 1.;
360  gion = 2.;
361  gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
362  /* >>chng 00 dec 15 change limit to -1 of H edge */
363  limit = MIN2(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1,rfield.nflux);
364 
365  if( rfield.ContBoltz[hmi.iphmin-1] > 0. )
366  {
367  for( nu=hmi.iphmin-1; nu < limit; nu++ )
368  {
369  /* H- flux photons cm-3 s-1
370  * ContBoltz is ratio of Boltzmann factor for each freq */
371  factor = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[hmi.iphmin-1]*rfield.widflx[nu]*
372  opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
373  rfield.anu2[nu]*fac;
374  rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
375  rfield.DiffuseEscape[nu] += (realnum)factor;
376  }
377  }
378  else
379  {
380  for( nu=hmi.iphmin-1; nu < limit; nu++ )
381  {
382  double arg = MAX2(0.,TE1RYD*(rfield.anu[nu]-0.05544)/phycon.te);
383  /* this is the limit sexp normally uses */
384  if( arg > SEXP_LIMIT )
385  break;
386  /* H- flux photons cm-3 s-1
387  * flux is in photons per sec per ryd */
388  factor = gamma*exp(-arg)*rfield.widflx[nu]*
389  opac.OpacStack[nu-hmi.iphmin+opac.iphmop]*
390  rfield.anu2[nu]*fac;
391  rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
392  rfield.DiffuseEscape[nu] += (realnum)factor;
393  }
394  }
395 
396  /* outward level 1 line photons, 0 is dummy line */
397  for( long i=1; i <= nLevel1; i++ )
398  TauLines[i].outline_resonance();
399 
400  for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
401  {
402  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
403  {
404  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
405  {
406  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_local; i++ )
407  {
408  const TransitionList::iterator& tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][i];
409  (*tr).Emis().phots() =
410  (*tr).Emis().Aul()*
411  (*(*tr).Hi()).Pop()*
412  ((*tr).Emis().Pesc()+
413  (*tr).Emis().Pelec_esc());
414 
415  (*tr).outline_resonance();
416  }
417  }
418  }
419  }
420 
421  /* outward level 2 line photons */
422  for( long i=0; i < nWindLine; i++ )
423  {
424  /* must not also do lines that were already done as part
425  * of the isoelectronic sequences */
426  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
427  {
428  {
429  enum {DEBUG_LOC=false};
430  if( DEBUG_LOC /*&& nzone > 10*/ && i==4821 )
431  {
432  /* set up to dump the Fe 9 169A line */
433  fprintf(ioQQQ,"DEBUG dump lev2 line %li\n", i );
434  DumpLine( TauLine2[i] );
435  fprintf(ioQQQ,"DEBUG dump %.3e %.3e %.3e\n",
436  rfield.outlin[0][TauLine2[i].ipCont()-1],
437  TauLine2[i].Emis().phots()*TauLine2[i].Emis().FracInwd()*radius.BeamInOut*opac.tmn[i]*TauLine2[i].Emis().ColOvTot(),
438  TauLine2[i].Emis().phots()*(1. - TauLine2[i].Emis().FracInwd())*radius.BeamOutOut* TauLine2[i].Emis().ColOvTot() );
439  }
440  }
441  TauLine2[i].outline_resonance();
442  /*if( i==2576 ) fprintf(ioQQQ,"DEBUG dump %.3e %.3e \n",
443  rfield.outlin[0][TauLine2[i].ipCont()-1] , rfield.outlin_noplot[TauLine2[i].ipCont()-1]);*/
444  }
445  }
446 
447  /* outward hyperfine structure line photons */
448  for( long i=0; i < nHFLines; i++ )
449  {
450  HFLines[i].outline_resonance();
451  }
452 
453  /* external database lines */
454  for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
455  {
456  if( dBaseSpecies[ipSpecies].lgActive )
457  {
458  for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
459  tr != dBaseTrans[ipSpecies].end(); ++tr)
460  {
461  int ipHi = (*tr).ipHi();
462  if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
463  continue;
464  (*tr).outline_resonance();
465  }
466  }
467  }
468 
469  /* H2 emission */
470  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
471  (*diatom)->H2_RT_diffuse();
472 
473  /* do outward parts of FeII lines, if large atom is turned on */
474  FeII_RT_Out();
478  if( trace.lgTrace )
479  fprintf( ioQQQ, " RT_diffuse returns.\n" );
480 
481  /* >>chng 02 jul 25, zero out all light below plasma freq */
482  for( nu=0; nu < rfield.ipPlasma-1; nu++ )
483  {
484  rfield.flux_beam_const[nu] = 0.;
485  rfield.flux_beam_time[nu] = 0.;
486  rfield.flux_isotropic[nu] = 0.;
487  rfield.flux[0][nu] = 0.;
488  rfield.ConEmitLocal[nzone][nu] = 0.;
489  rfield.otscon[nu] = 0.;
490  rfield.otslin[nu] = 0.;
491  rfield.outlin[0][nu] = 0.;
492  rfield.outlin_noplot[nu] = 0.;
493  rfield.reflin[0][nu] = 0.;
494  rfield.TotDiff2Pht[nu] = 0.;
495  rfield.ConInterOut[nu] = 0.;
496  }
497 
498  /* find occupation number, also assert that no continua are negative */
499  for( nu=0; nu < rfield.nflux; nu++ )
500  {
501  /* >>chng 00 oct 03, add diffuse continua */
502  /* local diffuse continua */
503  rfield.OccNumbDiffCont[nu] =
504  rfield.ConEmitLocal[nzone][nu]*rfield.convoc[nu];
505 
506  /* units are photons cell-1 cm-2 s-1 */
507  rfield.ConSourceFcnLocal[nzone][nu] =
508  /* units photons cm-3 s-1 cell-1, */
509  (realnum)safe_div( (double)rfield.ConEmitLocal[nzone][nu],
510  /* units cm-1 */
511  opac.opacity_abs[nu] );
512 
513  /* confirm that all are non-negative */
514  ASSERT( rfield.flux_beam_const[nu] >= 0.);
515  ASSERT( rfield.flux_beam_time[nu] >= 0.);
516  ASSERT( rfield.flux_isotropic[nu] >= 0.);
517  ASSERT( rfield.flux[0][nu] >= 0.);
518  ASSERT( rfield.ConEmitLocal[nzone][nu] >= 0.);
519  ASSERT( rfield.otscon[nu] >= 0.);
520  ASSERT( rfield.otslin[nu] >= 0.);
521  ASSERT( rfield.outlin[0][nu] >= 0.);
522  ASSERT( rfield.outlin_noplot[nu] >= 0.);
523  ASSERT( rfield.reflin[0][nu] >= 0.);
524  ASSERT( rfield.TotDiff2Pht[nu] >= 0.);
525  ASSERT( rfield.ConInterOut[nu] >= 0.);
526  }
527 
528  /* option to kill outward lines with no outward lines command*/
529  if( rfield.lgKillOutLine )
530  {
531  for( nu=0; nu < rfield.nflux; nu++ )
532  {
533  rfield.outlin[0][nu] = 0.;
534  rfield.outlin_noplot[nu] = 0.;
535  }
536  }
537 
538  /* option to kill outward continua with no outward continua command*/
539  if( rfield.lgKillOutCont )
540  {
541  for( nu=0; nu < rfield.nflux; nu++ )
542  {
543  rfield.ConInterOut[nu] = 0.;
544  }
545  }
546  return;
547 }
548 
549 void RT_iso_integrate_RRC( const long ipISO, const long nelem, const bool lgUpdateContinuum )
550 {
551  DEBUG_ENTRY( "RT_iso_integrate_RRC()" );
552 
553  // this array stores the last temperature at which cooling coefficients were evaluated
554  static double TeUsed[NISO][LIMELM]={
555  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
556  {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
557 
558  if( !lgUpdateContinuum && fp_equal( phycon.te, TeUsed[ipISO][nelem] ) && conv.nTotalIoniz )
559  return;
560 
561  ASSERT( nelem >= ipISO );
562  ASSERT( nelem < LIMELM );
563 
564  /* this will be the sum of recombinations to all excited levels */
565  double SumCaseB = 0.;
566 
567  /* the product of the densities of the parent ion and electrons */
568  double EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
569 
570  // recombination continua for all iso seq -
571  // if this stage of ionization exists
572  if( dense.IonHigh[nelem] >= nelem+1-ipISO )
573  {
574  t_iso_sp* sp = &iso_sp[ipISO][nelem];
575 
576  // loop over all levels to include recombination diffuse continua,
577  // pick highest energy continuum point that opacities extend to
578  long ipHi = rfield.nflux;
579  // >>chng 06 aug 17, should go to numLevels_local instead of _max.
580  for( long n=0; n < sp->numLevels_local; n++ )
581  {
582  double Sum1level = 0.;
583  sp->fb[n].RadRecCon = 0.;
584  sp->fb[n].RadRecCoolCoef = 0.;
585  // the number is (2 pi me k/h^2) ^ -3/2 * 8 pi/c^2 / ge - it includes
586  // the stat weight of the free electron in the demominator
587  double gamma = 0.5*MILNE_CONST*sp->st[n].g()/iso_ctrl.stat_ion[ipISO]/phycon.te/phycon.sqrte;
588 
589  // loop over all recombination continua
590  // escaping part of recombinations are added to rfield.ConEmitLocal
591  // added to ConInterOut at end of routine
592  for( long nu=sp->fb[n].ipIsoLevNIonCon-1; nu < ipHi; nu++ )
593  {
594  // dwid used to adjust where within WIDFLX exp is evaluated -
595  // weighted to lower energy due to exp(-energy/T)
596  double dwid = 0.2;
597 
598  // this is the term in the negative exponential Boltzmann factor
599  // for continuum emission
600  double arg = (rfield.anu[nu]-sp->fb[n].xIsoLevNIonRyd+
601  rfield.widflx[nu]*dwid)/phycon.te_ryd;
602  arg = MAX2(0.,arg);
603  // don't bother evaluating for this or higher energies if
604  // Boltzmann factor is tiny.
605  if( arg > SEXP_LIMIT )
606  break;
607 
608  /* photon is in photons cm^3 s^-1 per cell */
609  double photon = gamma*exp(-arg)*rfield.widflx[nu]*
610  opac.OpacStack[ nu-sp->fb[n].ipIsoLevNIonCon + sp->fb[n].ipOpac ] *
611  rfield.anu2[nu];
612 
613  Sum1level += photon;
614 
615  fixit(); // We need to include induced recombination in diffuse spectrum.
616  // Probably best just to do the whole thing here, then no need for induced terms in iso_photo.
617 
618  if( lgUpdateContinuum )
619  {
620  /* total local diffuse emission units photons cm-3 s-1 cell-1,*/
621  rfield.ConEmitLocal[nzone][nu] += (realnum)(photon*EdenAbund);
622 
623  // sp->fb[n].RadRecomb[ipRecEsc] is escape probability
624  // rfield.DiffuseEscape is local emission that escapes this zone
625  rfield.DiffuseEscape[nu] +=
626  (realnum)(photon*EdenAbund*sp->fb[n].RadRecomb[ipRecEsc]);
627  }
628 
629  // total RRC radiative recombination continuum
630  sp->fb[n].RadRecCon += rfield.anu[nu] *
631  emergent_line( photon*EdenAbund/2. , photon*EdenAbund/2. ,
632  // energy on fortran scale
633  nu+1 );
634 
635  double energyAboveThresh = rfield.anu[nu] - sp->fb[n].xIsoLevNIonRyd;
636  energyAboveThresh = MAX2( 0., energyAboveThresh );
637  sp->fb[n].RadRecCoolCoef += energyAboveThresh * photon *
638  sp->fb[n].RadRecomb[ipRecNetEsc];
639  }
640 
641  // convert to erg cm-3 s-1
642  sp->fb[n].RadRecCon *= EN1RYD;
643  sp->fb[n].RadRecCoolCoef *= EN1RYD;
644  /* this will be used below to confirm case B sum */
645  if( n > 0 )
646  {
647  /* SumCaseB will be sum to all excited */
648  SumCaseB += Sum1level;
649  }
650  }
651 
652  // no RRC emission from levels that do not exist
653  for( long n=sp->numLevels_local;n<sp->numLevels_max; n++ )
654  {
655  sp->fb[n].RadRecCon = 0.;
656  sp->fb[n].RadRecCoolCoef = 0.;
657  }
658 
659  /* this is check on self-consistency */
660  sp->CaseBCheck = MAX2(sp->CaseBCheck,
661  (realnum)(SumCaseB/sp->RadRec_caseB));
662  }
663 
664  TeUsed[ipISO][nelem] = phycon.te;
665 
666  return;
667 }
668 
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:578
t_Heavy::ipLyHeavy
long int ipLyHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:13
RT_diffuse
void RT_diffuse(void)
Definition: rt_diffuse.cpp:34
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
t_CoolHeavy::lgFreeOn
bool lgFreeOn
Definition: coolheavy.h:116
t_rfield::gff
realnum ** gff
Definition: rfield.h:227
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
t_dense::eden
double eden
Definition: dense.h:190
DumpLine
void DumpLine(const TransitionProxy &t)
Definition: transition.cpp:100
dense
t_dense dense
Definition: dense.cpp:24
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_rfield::flux
realnum ** flux
Definition: rfield.h:86
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
nLevel1
long int nLevel1
Definition: taulines.cpp:28
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:281
diatom_iter
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
realnum
float realnum
Definition: cddefines.h:103
conv.h
rfield.h
t_rfield::ConSourceFcnLocal
realnum ** ConSourceFcnLocal
Definition: rfield.h:152
nWindLine
long nWindLine
Definition: cdinit.cpp:19
t_iso_sp::CaseBCheck
realnum CaseBCheck
Definition: iso.h:510
mole.h
t_rfield::convoc
realnum * convoc
Definition: rfield.h:134
t_opac::iphmop
long int iphmop
Definition: opacity.h:226
t_rfield::flux_isotropic
realnum * flux_isotropic
Definition: rfield.h:89
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
diatoms
vector< diatomics * > diatoms
Definition: h2.cpp:8
t_rfield::otscon
realnum * otscon
Definition: rfield.h:195
t_rfield::ipPlasma
long int ipPlasma
Definition: rfield.h:453
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
t_mole_global::list
MoleculeList list
Definition: mole.h:317
SEXP_LIMIT
const double SEXP_LIMIT
Definition: cddefines.h:1476
ipoint.h
ProxyIterator
Definition: proxy_iterator.h:58
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:199
GrainVar::lgGrainPhysicsOn
bool lgGrainPhysicsOn
Definition: grainvar.h:479
t_iso_sp::TwoNu
vector< two_photon > TwoNu
Definition: iso.h:586
lines_service.h
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
t_iso_sp
Definition: iso.h:441
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
opac
t_opac opac
Definition: opacity.cpp:5
atmdat.h
POW2
#define POW2
Definition: cddefines.h:929
MIN2
#define MIN2
Definition: cddefines.h:761
nzone
long int nzone
Definition: cddefines.cpp:14
t_Heavy::RadRecCon
double RadRecCon[LIMELM][LIMELM]
Definition: heavy.h:18
radius
t_radius radius
Definition: radius.cpp:5
t_rfield::DiffuseEscape
realnum * DiffuseEscape
Definition: rfield.h:184
t_rfield::ConEmitLocal
realnum ** ConEmitLocal
Definition: rfield.h:149
t_radius::dVolOutwrd
double dVolOutwrd
Definition: radius.h:97
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
t_hmi::iphmin
long int iphmin
Definition: hmi.h:117
dense.h
coolheavy.h
mole
t_mole_local mole
Definition: mole.cpp:7
t_isoCTRL::lgDielRecom
bool lgDielRecom[NISO]
Definition: iso.h:365
trace
t_trace trace
Definition: trace.cpp:5
EmissionProxy::phots
double & phots() const
Definition: emission.h:503
cddefines.h
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:961
PrtTwoPhotonEmissCoef
void PrtTwoPhotonEmissCoef(const two_photon &tnu, const double &densityProduct)
Definition: two_photon.cpp:157
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
t_opac::OpacStack
double * OpacStack
Definition: opacity.h:151
MILNE_CONST
const UNUSED double MILNE_CONST
Definition: physconst.h:233
t_opac::iphmra
long int iphmra
Definition: opacity.h:223
t_isoCTRL::stat_ion
realnum stat_ion[NISO]
Definition: iso.h:362
radius.h
grains.h
t_rfield::nflux
long int nflux
Definition: rfield.h:43
heavy.h
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
hmi.h
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
t_rfield::lgKillOutCont
bool lgKillOutCont
Definition: rfield.h:437
t_opac::tmn
realnum * tmn
Definition: opacity.h:136
MAX2
#define MAX2
Definition: cddefines.h:782
RT_iso_integrate_RRC
void RT_iso_integrate_RRC(const long ipISO, const long nelem, const bool lgUpdateContinuum)
Definition: rt_diffuse.cpp:549
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_radius::BeamOutOut
double BeamOutOut
Definition: radius.h:108
TransitionProxy::outline
void outline(double nonScatteredFraction, bool lgDoChecks) const
Definition: transition.cpp:44
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
t_iso_sp::st
qList st
Definition: iso.h:453
t_rfield::OccNumbDiffCont
realnum * OccNumbDiffCont
Definition: rfield.h:141
FeII_RT_Out
void FeII_RT_Out(void)
Definition: atom_feii.cpp:2542
t_phycon::te_ryd
double te_ryd
Definition: phycon.h:17
emergent_line
double emergent_line(double emissivity_in, double emissivity_out, long int ipCont)
Definition: lines_service.cpp:335
t_rfield::nupper
long int nupper
Definition: rfield.h:46
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
t_rfield::anu2
realnum * anu2
Definition: rfield.h:79
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_rfield::otslin
realnum * otslin
Definition: rfield.h:193
TauLines
TransitionList TauLines("TauLines", &AnonStates)
t_rfield::reflin
realnum ** reflin
Definition: rfield.h:206
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:513
grainvar.h
rt.h
t_rfield::anu
double * anu
Definition: rfield.h:58
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
hmi
t_hmi hmi
Definition: hmi.cpp:5
physconst.h
SatelliteLines
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:38
t_rfield::lgKillOutLine
bool lgKillOutLine
Definition: rfield.h:434
conv
t_conv conv
Definition: conv.cpp:5
t_rfield::flux_beam_const
realnum * flux_beam_const
Definition: rfield.h:92
atmdat_outer_shell
void atmdat_outer_shell(long int iz, long int in, long int *imax, long int *ig0, long int *ig1)
Definition: atmdat_outer_shell.cpp:8
gv
GrainVar gv
Definition: grainvar.cpp:5
t_opac::ExpmTau
realnum * ExpmTau
Definition: opacity.h:132
t_rfield::flux_beam_time
realnum * flux_beam_time
Definition: rfield.h:92
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
fixit
void fixit(void)
Definition: service.cpp:991
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:164
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:200
GrainMakeDiffuse
void GrainMakeDiffuse(void)
Definition: grains_qheat.cpp:170
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
CoolHeavy
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
nHFLines
long int nHFLines
Definition: taulines.cpp:31
phycon.h
t_rfield::ContBoltz
double * ContBoltz
Definition: rfield.h:145
t_radius::BeamInOut
double BeamInOut
Definition: radius.h:105
ipH1s
const int ipH1s
Definition: iso.h:27
t_rfield::DiffuseLineEmission
realnum * DiffuseLineEmission
Definition: rfield.h:203
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
TE1RYD
const UNUSED double TE1RYD
Definition: physconst.h:183
ipSatelliteLines
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:37
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
h2.h
nSpecies
long int nSpecies
Definition: taulines.cpp:21
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
t_rfield::ipMaxBolt
long int ipMaxBolt
Definition: rfield.h:249
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:279
t_isoCTRL::lgInd2nu_On
bool lgInd2nu_On
Definition: iso.h:355
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:269
t_Heavy::ipHeavy
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
GrainVar::lgDustOn
bool lgDustOn() const
Definition: grainvar.h:471
t_rfield::lgInducProcess
bool lgInducProcess
Definition: rfield.h:252
t_Heavy::xLyaHeavy
realnum xLyaHeavy[LIMELM][LIMELM]
Definition: heavy.h:21
EN1RYD
const UNUSED double EN1RYD
Definition: physconst.h:179
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
atomfeii.h
t_Heavy::ipBalHeavy
long int ipBalHeavy[LIMELM][LIMELM-1]
Definition: heavy.h:15
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
t_rfield::TotDiff2Pht
realnum * TotDiff2Pht
Definition: rfield.h:187
CalcTwoPhotonEmission
void CalcTwoPhotonEmission(two_photon &tnu, bool lgDoInduced)
Definition: two_photon.cpp:125
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
HFLines
TransitionList HFLines("HFLines", &AnonStates)
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
Generated by   doxygen 1.8.17