cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
rt_diffuse.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_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 
29 #define TwoS (1+ipISO)
30 
31 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
32 #pragma optimization_level 0
33 #endif
34 void RT_diffuse(void)
35 {
36  /* set this true to print a table of 2 photon emission coefficients
37  * comparable to Brown & Mathews (1971) */
38  static long lgPrt2PhtEmisCoef = false;
39 
40  /* arrays used in this routine
41  * ARRAY rfield.ConEmitLocal[depth][energy] local emission per unit vol
42  * ARRAY rfield.ConOTS_local_photons - local ots two-photon rates
43  * ARRAY DiffuseEscape is the spectrum of diffuse emission that escapes this zone,
44  * at end of this routine part is thrown into the outward beam
45  * by adding to rfield.ConInterOut
46  * units are photons s-1 cm-3
47  * one-time init done on first call */
48  static realnum *DiffuseEscape;
49 
50  /* DiffuseEscape and rfield.ConEmitLocal are same except that
51  * rfield.ConEmitLocal is local emission, would be source function if div by opac
52  * DiffuseEscape is part that escapes so has RT built into it
53  * DiffuseEscape is used to define rfield.ConInterOut below as per this statement
54  * rfield.ConInterOut[nu] += DiffuseEscape[nu]*(realnum)radius.dVolOutwrd;
55  */
56  /* \todo 0 define only rfield.ConEmitLocal as it is now done,
57  * do not define DiffuseEscape at all
58  * at bottom of this routine use inward and outward optical depths to define
59  * local and escaping parts
60  * this routine only defines
61  * rfield.ConInterOut - set to DiffuseEscape times vol element
62  * rfield.ConOTS_local_photons (two photon only)
63  * so this is only var that
64  * needs to be set
65  */
66  static bool lgMustInit=true;
67 
68  long int i=-100000,
69  ip=-100000,
70  ipISO=-100000,
71  ipHi=-100000,
72  ipLo=-100000,
73  ipla=-100000,
74  nelem=-100000,
75  ion=-100000,
76  limit=-100000,
77  n=-100000,
78  nu=-10000;
79 
80  double EnergyLimit;
81 
82  double EdenAbund,
83  difflya,
84  xInWrd,
85  arg,
86  fac,
87  factor,
88  gamma,
89  gion,
90  gn,
91  photon,
92  sum,
93  Sum1level,
94  SumCaseB;
95 
96  realnum Dilution , two_photon_emiss;
97 
98  DEBUG_ENTRY( "RT_diffuse()" );
99 
100  /* one time creation of space for ThowOut array */
101  if( lgMustInit )
102  {
103  lgMustInit = false;
104  DiffuseEscape = (realnum*)MALLOC((unsigned)rfield.nupper*sizeof(realnum));
105  }
106 
107  /* many arrays were malloced to nupper, and we will add unit flux to [nflux] -
108  8 this must be true to work */
110 
111  /* this routine evaluates the local diffuse fields
112  * it fills in all of the following vectors */
113  memset(DiffuseEscape , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
114  memset(rfield.ConEmitLocal[nzone] , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
115  memset(rfield.ConOTS_local_photons , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
116  memset(rfield.TotDiff2Pht , 0 , (unsigned)rfield.nupper*sizeof(realnum) );
117 
118  /* must abort after setting all of above to zero because some may be
119  * used in various ways before abort is complete */
120  if( lgAbort )
121  {
122  /* quit if we are aborting */
123  return;
124  }
125 
126  /* loop over iso-sequences of all elements
127  * to add all recombination continua and lines*/
128  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
129  {
130  /* >>chng 01 sep 23, rewrote for iso sequences */
131  for( nelem=ipISO; nelem < LIMELM; nelem++ )
132  {
133  /* this will be the sum of recombinations to all excited levels */
134  SumCaseB = 0.;
135 
136  /* the product of the densities of the parent ion and electrons */
137  EdenAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
138 
139  /* recombination continua for all iso seq -
140  * if this stage of ionization exists */
141  if( dense.IonHigh[nelem] >= nelem+1-ipISO )
142  {
143  /* loop over all levels to include recombination diffuse continua,
144  * pick highest energy continuum point that opacities extend to
145  * for ground continuum, go up to highest defined Boltzmann factor,
146  * at bottom of loop will be reset to ground state photo edge */
147  ipHi = rfield.nflux;
148  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
149  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
150  {
151  Sum1level = 0.;
152 
153  /* >>chng 02 nov 20 - pull the plug on old he treatment */
154  /* the number is (2 pi me k/h^2) ^ -3/2 * 8 pi/c^2 / ge - it includes
155  * the stat weight of the free electron in the demominator */
157  /*gamma = 2.0618684e11*StatesElem[ipISO][nelem][n].g/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte;*/
158  gamma = 0.5*MILNE_CONST*StatesElem[ipISO][nelem][n].g/iso.stat_ion[ipISO]/phycon.te/phycon.sqrte;
159 
160  /* loop over all recombination continua
161  * escaping part of recombinations are added to rfield.ConEmitLocal
162  * added to ConInterOut at end of routine */
163  for( nu=iso.ipIsoLevNIonCon[ipISO][nelem][n]-1; nu < ipHi; nu++ )
164  {
165  /* dwid used to adjust where within WIDFLX exp is evaluated -
166  * weighted to lower energy due to exp(-energy/T) */
167  double dwid = 0.2;
168 
169  /* this is the term in the negative exponential Boltzmann factor
170  * for continuum emission */
171  arg = (rfield.anu[nu]-iso.xIsoLevNIonRyd[ipISO][nelem][n]+
172  rfield.widflx[nu]*dwid)/phycon.te_ryd;
173  arg = MAX2(0.,arg);
174  /* don't bother evaluating for this or higher energies if
175  * Boltzmann factor is tiny. */
176  if( arg > SEXP_LIMIT )
177  break;
178 
179  /* photon is in photons cm^3 s^-1 per cell */
180  photon = gamma*exp(-arg)*rfield.widflx[nu]*
181  opac.OpacStack[ nu-iso.ipIsoLevNIonCon[ipISO][nelem][n] + iso.ipOpac[ipISO][nelem][n] ] *
182  rfield.anu2[nu];
183 
184  /*ASSERT( photon >= 0. );*/
185  Sum1level += photon;
186  /* total local diffuse emission */
187  rfield.ConEmitLocal[nzone][nu] += (realnum)(photon*EdenAbund);
188  /* iso.RadRecomb[ipISO][nelem][n][ipRecEsc] is escape probability
189  * DiffuseEscape is local emission that escapes this zone */
190  DiffuseEscape[nu] +=
191  (realnum)(photon*EdenAbund*iso.RadRecomb[ipISO][nelem][n][ipRecEsc]);
192  }
193  /* this will be used below to confirm case B sum */
194  if( n > 0 )
195  {
196  /* SumCaseB will be sum to all excited */
197  SumCaseB += Sum1level;
198  }
199 
200  /* on entry to this loop ipHi was set to the upper limit of the code,
201  * and ground state recombination continua for all energies was added. For
202  * excited states (which are the ones that will be done after
203  * first pass through) only go to ground state threshold,
204  * since that will be so much stronger than excited state recom*/
205  ipHi = iso.ipIsoLevNIonCon[ipISO][nelem][0]-1;
206  }
207  /* this is check on self-consistency */
208  iso.CaseBCheck[ipISO][nelem] = MAX2(iso.CaseBCheck[ipISO][nelem],
209  (realnum)(SumCaseB/iso.RadRec_caseB[ipISO][nelem]));
210 
211  /* this add line emission from the model atoms,
212  * do not include very highest level since disturbed by topoff */
213  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
214  for( ipLo=0; ipLo < (iso.numLevels_local[ipISO][nelem] - 2); ipLo++ )
215  {
216  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
217  for( ipHi=ipLo+1; ipHi < iso.numLevels_local[ipISO][nelem] - 1; ipHi++ )
218  {
219  /* must not include fake transitions (2s-1s, two photon, or truly
220  * forbidden transitions */
221  if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 1 )
222  continue;
223 
224  /* pointer to line energy in continuum array */
225  ip = Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1;
226 
227  /* number of photons in the line has not been defined up until now,
228  * do so now. this is redone in lines. */
229  Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots =
230  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul*
231  StatesElem[ipISO][nelem][ipHi].Pop*
232  Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc*
233  dense.xIonDense[nelem][nelem+1-ipISO];
234 
235  /* inward fraction of line */
236  xInWrd = Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots*
237  Transitions[ipISO][nelem][ipHi][ipLo].Emis->FracInwd;
238 
239  /* reflected part of line */
240  rfield.reflin[ip] += (realnum)(xInWrd*radius.BeamInIn);
241 
242  /* inward beam that goes out since sphere set */
243  /* in all this the ColOvTot term has been commented out,
244  * since this is not meaningful for something like the hydrogen atom,
245  * where most forms by recombination */
246  rfield.outlin[ip] += (realnum)(xInWrd*radius.BeamInOut*opac.tmn[ip]/*
247  Transitions[ipISO][nelem][ipHi][ipLo].ColOvTot*/);
248 
249  /* outward part of line */
250  rfield.outlin[ip] +=
251  (realnum)(Transitions[ipISO][nelem][ipHi][ipLo].Emis->phots*
252  (1. - Transitions[ipISO][nelem][ipHi][ipLo].Emis->FracInwd)*
253  radius.BeamOutOut*opac.tmn[ip]/*
254  Transitions[ipISO][nelem][ipHi][ipLo].ColOvTot*/);
255  }
256  }
257 
258  /*Iso treatment of two photon emission. */
259  /* NISO could in the future be increased, but we want this assert to blow
260  * so that it is understood this may not be correct for other iso sequences,
261  * probably should break since will not be present */
262  ASSERT( ipISO <= ipHE_LIKE );
263 
264  /* upper limit to 2-phot is energy of 2s to ground */
265  EnergyLimit = Transitions[ipISO][nelem][TwoS][0].EnergyWN/RYD_INF;
266 
267  /* this could fail when pops very low and Pop2Ion is zero */
268  ASSERT( iso.ipTwoPhoE[ipISO][nelem]>0 && EnergyLimit>0. );
269 
270  sum = 0.;
271  /* two photon emission, iso.ipTwoPhoE[ipISO][nelem] is
272  * continuum array index for Lya energy */
273  for( nu=0; nu<iso.ipTwoPhoE[ipISO][nelem]; nu++ )
274  {
275  /* We do not assert rfield.anu[nu]<=EnergyLimit because the maximum
276  * index could be set to point to the next highest bin. */
277  ASSERT( rfield.anu[nu]/EnergyLimit < 1.01 || rfield.anu[nu-1]<EnergyLimit);
278 
279  /* iso.As2nu[ipISO][nelem][nu] is transition probability A per bin
280  * So sum is the total transition probability - this sum should
281  * add up to the A two photon */
282  sum += iso.As2nu[ipISO][nelem][nu];
283 
284  /* flag saying whether to also do he triplets */
285  const bool lgDo2TripSToo = false;
286 
287  if( ipISO == ipHE_LIKE && lgDo2TripSToo )
288  {
289  two_photon_emiss = 2.f*dense.xIonDense[nelem][nelem+1-ipISO]*iso.As2nu[ipISO][nelem][nu]*
290  ((realnum)StatesElem[ipISO][nelem][TwoS].Pop+0.0001f*(realnum)StatesElem[ipISO][nelem][ipHe2s3S].Pop);
291  }
292  else
293  {
294  /* StatesElem[ipISO][nelem][TwoS].Pop is dimensionless and represents the population
295  * of the TwoS level relative to that of the ion. dense.xIonDense[nelem][nelem+1-ipISO]
296  * is the density of the current ion (cm^-3). The factor of 2 is for two photons per
297  * transition. Thus two_photon_emiss has dimension photons cm-3 s-1. */
298  two_photon_emiss = 2.f*(realnum)StatesElem[ipISO][nelem][TwoS].Pop*dense.xIonDense[nelem][nelem+1-ipISO]
299  *iso.As2nu[ipISO][nelem][nu];
300  }
301 
302  /* >>chng 03 mar 26, only do if induced processes turned on,
303  * otherwise inconsistent with rate solver treatment. */
304  /* >>chng 02 aug 14, add induced two photon emission */
305  /* product of occupation numbers for induced emission */
306  if( iso.lgInd2nu_On)
307  {
308  /* this is the total rate */
309  two_photon_emiss *= (1.f + rfield.SummedOcc[nu]) *
310  (1.f+rfield.SummedOcc[iso.ipSym2nu[ipISO][nelem][nu]-1]);
311  }
312 
313  /* information - only used in punch output */
314  rfield.TotDiff2Pht[nu] += two_photon_emiss;
315 
316  /* total local diffuse emission */
317  rfield.ConEmitLocal[nzone][nu] += two_photon_emiss;
318 
319  /* this is escaping part of two-photon emission,
320  * as determined from optical depth to illuminated face */
321  DiffuseEscape[nu] += two_photon_emiss*opac.ExpmTau[nu];
322 
323  /* save locally destroyed OTS two-photon continuum - this is only
324  * place this is set in this routine */
325  rfield.ConOTS_local_photons[nu] += two_photon_emiss*(1.f - opac.ExpmTau[nu]);
326  }
327 
328  /* a sanity check on the code, see Spitzer and Greenstein (1951), eqn 4. */
329  /* >>refer HI 2nu Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
330  ASSERT( fabs( 1.f - (realnum)sum/Transitions[ipISO][nelem][TwoS][0].Emis->Aul ) < 0.01f );
331  }
332 
333  /* option to print hydrogen and helium two-photon emission coefficients? */
334  if( lgPrt2PhtEmisCoef )
335  {
336  long yTimes20;
337  double y, E2nu;
338 
339  lgPrt2PhtEmisCoef = false;
340 
341  fprintf( ioQQQ, "\ny\tGammaNot(2q)\n");
342 
343  for( yTimes20=1; yTimes20<=10; yTimes20++ )
344  {
345  y = yTimes20/20.;
346 
347  fprintf( ioQQQ, "%.3e\t", (double)y );
348 
349  E2nu = Transitions[0][0][1][0].EnergyWN/RYD_INF;
350  i = ipoint(y*E2nu);
351  fprintf( ioQQQ, "%.3e\t",
352  8./3.*HPLANCK*StatesElem[0][0][1].Pop/dense.eden
353  *y*iso.As2nu[0][0][i]*E2nu/rfield.widflx[i] );
354 
355  E2nu = Transitions[1][1][2][0].EnergyWN/RYD_INF;
356  i = ipoint(y*E2nu);
357  fprintf( ioQQQ, "%.3e\n",
358  8./3.*HPLANCK*StatesElem[1][1][2].Pop/dense.eden
359  *y*iso.As2nu[1][1][i]*E2nu/rfield.widflx[i] );
360 
361  /*
362  fprintf( ioQQQ, "%.3e\t%.3e\n",
363  rfield.TotDiff2Pht[i]*rfield.anu[i]*EN1RYD
364  /E2nu/rfield.widflx[i]/dense.eden/dense.xIonDense[nelem][nelem+1-ipISO]/FR1RYD,
365  8./3.*HPLANCK*StatesElem[ipISO][nelem][TwoS].Pop/dense.eden
366  *y*iso.As2nu[ipISO][nelem][i]*E2nu/rfield.widflx[i]); */
367  }
368  fprintf( ioQQQ, "eden%.3e\n", dense.eden );
369  }
370  }
371  }
372 
373  /* add recombination continua for elements heavier than those done with iso seq */
374  for( nelem=NISO; nelem < LIMELM; nelem++ )
375  {
376  /* do not include species with iso-sequence in following */
377  /* >>chng 03 sep 09, upper bound was wrong, did not include NISO */
378  for( ion=dense.IonLow[nelem]; ion < nelem-NISO+1; ion++ )
379  {
380  if( dense.xIonDense[nelem][ion+1] > 0. )
381  {
382  long int ns, nshell,igRec , igIon,
383  iplow , iphi , ipop;
384 
385  ip = Heavy.ipHeavy[nelem][ion]-1;
386  ASSERT( ip >= 0 );
387 
388  /* nflux was reset upward in ConvInitSolution to encompass all
389  * possible line and continuum emission. this test should not
390  * possibly fail. It could if the ionization were to increase with depth
391  * although the continuum mesh is designed to deal with this.
392  * This test is important because the nflux cell in ConInterOut
393  * is used to carry out the unit integration, and if it gets
394  * clobbered by diffuse emission the code will declare
395  * insanity in PrtComment */
396  if( ip >= rfield.nflux )
397  continue;
398 
399  /* get shell number, stat weights for this species */
400  atmdat_outer_shell( nelem+1 , nelem+1-ion , &nshell, &igRec , &igIon );
401  gn = (double)igRec;
402  gion = (double)igIon;
403 
404  /* shell number */
405  ns = Heavy.nsShells[nelem][ion]-1;
406  ASSERT( ns == (nshell-1) );
407 
408  /* lower and upper energies, and offset for opacity stack */
409  iplow = opac.ipElement[nelem][ion][ns][0]-1;
410  iphi = opac.ipElement[nelem][ion][ns][1];
411  iphi = MIN2( iphi , rfield.nflux );
412  ipop = opac.ipElement[nelem][ion][ns][2];
413 
414  /* now convert ipop to the offset in the opacity stack from threshold */
415  ipop = ipop - iplow;
416 
417  gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte*dense.eden*dense.xIonDense[nelem][ion+1];
418 
419  /* this is ground state continuum from stored opacities */
420  if( rfield.ContBoltz[iplow] > SMALLFLOAT )
421  {
422  for( nu=iplow; nu < iphi; ++nu )
423  {
424  photon = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[iplow]*
425  rfield.widflx[nu]*opac.OpacStack[nu+ipop]*rfield.anu2[nu];
426  /* add heavy rec to ground in active beam,*/
431  rfield.ConEmitLocal[nzone][nu] += (realnum)photon;
432  /*DiffuseEscape[i] += (realnum)photon;*/
433  /*>>chng 03 sep 08, index had been i, should have been nu */
434  DiffuseEscape[nu] += (realnum)photon*opac.ExpmTau[nu];
435  }
436  }
437 
438  /* now do the recombination Lya */
439  ipla = Heavy.ipLyHeavy[nelem][ion]-1;
440  ASSERT( ipla >= 0 );
441  /* xLyaHeavy is set to a fraction of the total rad rec in ion_recomb, includes eden */
442  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
443  /* >>chng 03 jul 10, here and below, use outlin_noplot */
444  rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
445 
446  /* now do the recombination Balmer photons */
447  ipla = Heavy.ipBalHeavy[nelem][ion]-1;
448  ASSERT( ipla >= 0 );
449  /* xLyaHeavy is set to fraction of total rad rec in ion_recomb, includes eden */
450  difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
451  rfield.outlin_noplot[ipla] += (realnum)(difflya*radius.dVolOutwrd*opac.tmn[ipla]*opac.ExpmTau[ipla]);
452  }
453  }
454  }
455 
456  /* free-free free free brems emission for all ions */
457  limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
458  for( nu=0; nu < limit; nu++ )
459  {
460  double TotBremsAllIons = 0., BremsThisIon;
461 
462  /* do hydrogen first, before main loop since want to add on H- brems */
463  nelem = ipHYDROGEN;
464  ion = 1;
465  BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
466  ASSERT( BremsThisIon >= 0. );
467 
468  /* for case of hydrogen, add H- brems - OpacStack contains the ratio
469  * of the H- to H brems cross section - multiply by this and the fraction in ground */
470  BremsThisIon *= (1.+opac.OpacStack[nu-1+opac.iphmra]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop);
471  TotBremsAllIons = BremsThisIon;
472 
473  /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
474  * using gaunt factors from rfield.gff. */
475  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
476  {
477  /* MAX2 occurs because we want to start at first ion (or above)
478  * and do not want atom */
479  for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
480  {
481  /* eff. charge is ion, so first rfield.gff argument must be "ion". */
482  BremsThisIon = POW2( (realnum)ion )*dense.xIonDense[nelem][ion]*rfield.gff[ion][nu];
483  /*ASSERT( BremsThisIon >= 0. );*/
484 
485  TotBremsAllIons += BremsThisIon;
486  }
487  }
488 
490  /* >>chng 06 apr 05, no free free also turns off emission */
491  TotBremsAllIons *= dense.eden*1.032e-11*rfield.widflx[nu]*rfield.ContBoltz[nu]/rfield.anu[nu]/phycon.sqrte *
493  ASSERT( TotBremsAllIons >= 0.);
494 
495  /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
496  * to outward beam - ConLocNoInter array removed as result
497  * if problems develop with very dense blr clouds, this may be reason */
498  /*rfield.ConLocNoInter[nu] += (realnum)fac;*/
499  /*rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;*/
500 
501  if( nu >= rfield.ipEnergyBremsThin )
502  {
503  /* >>chng 05 feb 20, move into this test on brems opacity - should not be
504  * needed since would use expmtau to limit outward beam */
505  /* >>chng 01 jul 01, move thick brems back to ConEmitLocal but do not add
506  * to outward beam - ConLocNoInter array removed as result
507  * if problems develop with very dense BLR clouds, this may be reason */
508  /*rfield.ConLocNoInter[nu] += (realnum)fac;*/
509  rfield.ConEmitLocal[nzone][nu] += (realnum)TotBremsAllIons;
510 
511  /* do not add optically thick part to outward beam since self absorbed
512  * >>chng 96 feb 27, put back into outward beam since do not integrate
513  * over it anyway. */
514  /* >>chng 99 may 28, take back out of beam since DO integrate over it
515  * in very dense BLR clouds */
516  /* >>chng 01 jul 10, add here, in only one loop, where optically thin */
517  DiffuseEscape[nu] += (realnum)TotBremsAllIons;
518  }
519  }
520 
521  /* grain dust emission */
522  /* >>chng 01 nov 22, moved calculation of grain flux to qheat.c, PvH */
523  if( gv.lgDustOn && gv.lgGrainPhysicsOn )
524  {
525  /* this calculates diffuse emission from grains,
526  * and stores the result in gv.GrainEmission */
528 
529  for( nu=0; nu < rfield.nflux; nu++ )
530  {
532  DiffuseEscape[nu] += gv.GrainEmission[nu];
533  }
534  }
535 
536  /* hminus emission */
537  fac = dense.eden*(double)dense.xIonDense[ipHYDROGEN][0];
538  gn = 1.;
539  gion = 2.;
540  gamma = 0.5*MILNE_CONST*gn/gion/phycon.te/phycon.sqrte;
541  /* >>chng 00 dec 15 change limit to -1 of H edge */
543 
544  if( rfield.ContBoltz[hmi.iphmin-1] > 0. )
545  {
546  for( nu=hmi.iphmin-1; nu < limit; nu++ )
547  {
548  /* H- flux photons cm-3 s-1
549  * ContBoltz is ratio of Boltzmann factor for each freq */
550  factor = gamma*rfield.ContBoltz[nu]/rfield.ContBoltz[hmi.iphmin-1]*rfield.widflx[nu]*
552  rfield.anu2[nu]*fac;
553  rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
554  DiffuseEscape[nu] += (realnum)factor;
555  }
556  }
557  else
558  {
559  for( nu=hmi.iphmin-1; nu < limit; nu++ )
560  {
561  arg = MAX2(0.,TE1RYD*(rfield.anu[nu]-0.05544)/phycon.te);
562  /* this is the limit sexp normally uses */
563  if( arg > SEXP_LIMIT )
564  break;
565  /* H- flux photons cm-3 s-1
566  * flux is in photons per sec per ryd */
567  factor = gamma*exp(-arg)*rfield.widflx[nu]*
569  rfield.anu2[nu]*fac;
570  rfield.ConEmitLocal[nzone][nu] += (realnum)factor;
571  DiffuseEscape[nu] += (realnum)factor;
572  }
573  }
574 
575  /* this dilution is needed to conserve volume in spherical models. tests such
576  * as parispn.in will fault if this is removed */
577  Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );
578 
579  /* this is a unit of energy that will be passed through the code as a test
580  * that all integrations are carried out. A similar test is set in lineset1
581  * and verified in PrtFinal. The opacity at this cell is zero so only
582  * geometrical dilution will affect the integral
583  * Radius is currently outer edge of zone, so radius-drad/2 is radius
584  * of center of zone */
585  rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution;
586  DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;
587 
588  /* opacity should be zero at this energy */
589  ASSERT( opac.opacity_abs[rfield.nflux] == 0. );
590 
591  /* many tmn added to conserve energy in very high Z models
592  * rerun highZ qso model if any tmn ever deleted here
593  *
594  * tmn set in ZoneStart and includes both opacity and dilution
595  *
596  * use diffuse lines and continuum to add flux to outward beam
597  *
598  * NB!!! this routine adds flux to the outward beam
599  * it can only be called once per zone
600  */
601 
602  /* >>chng 96 nov 19, do not consider energies below plasma freq
603  * ipPlasmaFreq points to lowest energy thin to brems and plasma freq */
604 
605  /* add DiffuseEscape continuum to ConInterOut,
606  * lower limit of rfield.ipPlasma-1 since continuum is zero below rfield.ipPlasma-1
607  * due to plasma frequency
608  * note that upper range of sum is <= nflux, which contains the unit
609  * verification token in the highest cell*/
610  for( nu=rfield.ipPlasma-1; nu <= rfield.nflux; nu++ )
611  {
612  /* ConInterOut is the interactive continuum
613  * DiffuseEscape contains all radiation thrown into outward beam */
614  /* radius.dVolOutwrd has factor or 1 or 0.5 for outward part */
615  /* >>chng 05 feb 23, activate this statement */
616  rfield.ConInterOut[nu] += DiffuseEscape[nu]*(realnum)radius.dVolOutwrd;
617  /*rfield.ConInterOut[nu] += rfield.ConEmitLocal[nzone][nu]*opac.ExpmTau[nu]*(realnum)radius.dVolOutwrd;*/
618  ASSERT( rfield.ConInterOut[nu] >= 0.);
619  }
620 
621  {
622  /*@-redef@*/
623  enum {DEBUG_LOC=false};
624  /*@+redef@*/
625  if( DEBUG_LOC )
626  {
627  fprintf(ioQQQ,"rtdiffusedebugg %li\t%.2e\t%.2e\t%.2e\n",
628  nzone, rfield.ConInterOut[1158] , DiffuseEscape[1158],radius.dVolOutwrd);
629  }
630  }
631 
632  /* outward level 1 line photons, 0 is dummy line */
633  for( i=1; i <= nLevel1; i++ )
634  {
635  outline( &TauLines[i] );
636  }
637 
638  for( ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
639  {
640  for( nelem = ipISO; nelem < LIMELM; nelem++ )
641  {
642  if( dense.lgElmtOn[nelem] )
643  {
644  for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ )
645  {
646 
647  SatelliteLines[ipISO][nelem][i].Emis->phots =
648  SatelliteLines[ipISO][nelem][i].Emis->Aul*
649  SatelliteLines[ipISO][nelem][i].Hi->Pop*
650  SatelliteLines[ipISO][nelem][i].Emis->Pesc*
651  dense.xIonDense[nelem][nelem+1-ipISO];
652 
653  outline( &SatelliteLines[ipISO][nelem][i] );
654  }
655  }
656  }
657  }
658 
659  /* outward level 2 line photons */
660  for( i=0; i < nWindLine; i++ )
661  {
662  /* must not also do lines that were already done as part
663  * of the isoelectronic sequences */
664  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
665  {
666  {
667  /*@-redef@*/
668  enum {DEBUG_LOC=false};
669  /*@+redef@*/
670  if( DEBUG_LOC /*&& nzone > 10*/ && i==4821 )
671  {
672  /* set up to dump the Fe 9 169A line */
673  fprintf(ioQQQ,"DEBUG dump lev2 line %li\n", i );
674  DumpLine( &TauLine2[i] );
675  fprintf(ioQQQ,"DEBUG dump %.3e %.3e %.3e\n",
679  }
680  }
681  outline( &TauLine2[i] );
682  /*if( i==2576 ) fprintf(ioQQQ,"DEBUG dump %.3e %.3e \n",
683  rfield.outlin[TauLine2[i].ipCont-1] , rfield.outlin_noplot[TauLine2[i].ipCont-1]);*/
684  }
685  }
686 
687  /* outward hyperfine structure line photons */
688  for( i=0; i < nHFLines; i++ )
689  {
690  outline( &HFLines[i] );
691  }
692 
693  /*Atoms & Molecules*/
694  for( i=0; i < linesAdded2; i++ )
695  {
696  outline( atmolEmis[i].tran );
697  }
698 
699  /* carbon monoxide */
700  for( i=0; i < nCORotate; i++ )
701  {
702  outline( &C12O16Rotate[i] );
703  outline( &C13O16Rotate[i] );
704  }
705 
706  /* H2 emission */
707  H2_RT_diffuse();
708 
709  /* do outward parts of FeII lines, if large atom is turned on */
710  FeII_RT_Out();
714  if( trace.lgTrace )
715  fprintf( ioQQQ, " RT_diffuse returns.\n" );
716 
717  /* >>chng 02 jul 25, zero out all light below plasma freq */
718  for( nu=0; nu<rfield.ipPlasma-1; nu++ )
719  {
720  rfield.flux_beam_const[nu] = 0.;
721  rfield.flux_beam_time[nu] = 0.;
722  rfield.flux_isotropic[nu] = 0.;
723  rfield.flux[nu] = 0.;
724  rfield.ConEmitLocal[nzone][nu] = 0.;
725  rfield.otscon[nu] =0.;
726  rfield.otslin[nu] =0.;
727  rfield.outlin[nu] =0.;
728  rfield.outlin_noplot[nu] =0.;
729  rfield.reflin[nu] = 0.;
730  rfield.TotDiff2Pht[nu] = 0.;
731  rfield.ConOTS_local_photons[nu] = 0.;
732  rfield.ConInterOut[nu] = 0.;
733  }
734 
735  /* find occupation number, also assert that no continua are negative */
736  for( nu=0; nu < rfield.nflux; nu++ )
737  {
738  /* >>chng 00 oct 03, add diffuse continua */
739  /* local diffuse continua */
741 
742  /* confirm that all are non-negative */
743  ASSERT( rfield.flux_beam_const[nu] >= 0.);
744  ASSERT( rfield.flux_beam_time[nu] >= 0.);
745  ASSERT( rfield.flux_isotropic[nu] >= 0.);
746  ASSERT( rfield.flux[nu] >=0.);
747  ASSERT( rfield.ConEmitLocal[nzone][nu] >= 0.);
748  ASSERT( rfield.otscon[nu] >=0.);
749  ASSERT( rfield.otslin[nu] >=0.);
750  ASSERT( rfield.outlin[nu] >=0.);
751  ASSERT( rfield.outlin_noplot[nu] >=0.);
752  ASSERT( rfield.reflin[nu] >= 0.);
753  ASSERT( rfield.TotDiff2Pht[nu] >= 0.);
754  ASSERT( rfield.ConOTS_local_photons[nu] >= 0.);
755  ASSERT( rfield.ConInterOut[nu] >=0.);
756  }
757 
758  /* option to kill outward lines with no outward lines command*/
759  if( rfield.lgKillOutLine )
760  {
761  for( nu=0; nu < rfield.nflux; nu++ )
762  {
763  rfield.outlin[nu] = 0.;
764  rfield.outlin_noplot[nu] = 0.;
765  }
766  }
767 
768  /* option to kill outward continua with no outward continua command*/
769  if( rfield.lgKillOutCont )
770  {
771  for( nu=0; nu < rfield.nflux; nu++ )
772  {
773  rfield.ConInterOut[nu] = 0.;
774  }
775  }
776  return;
777 }

Generated for cloudy by doxygen 1.8.4