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