cloudy trunk
Loading...
Searching...
No Matches
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
34void 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 */
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();
477
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
549void 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
void atmdat_outer_shell(long int iz, long int in, long int *imax, long int *ig0, long int *ig1)
void FeII_RT_Out(void)
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define ASSERT(exp)
Definition cddefines.h:578
const int ipRecEsc
Definition cddefines.h:279
const double SEXP_LIMIT
Definition cddefines.h:1476
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition cddefines.h:961
const int NISO
Definition cddefines.h:261
const int ipRecNetEsc
Definition cddefines.h:281
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
long nWindLine
Definition cdinit.cpp:19
realnum & Pesc() const
Definition emission.h:523
realnum & Aul() const
Definition emission.h:613
double & phots() const
Definition emission.h:503
TransitionProxy::iterator iterator
Definition transition.h:280
void outline(double nonScatteredFraction, bool lgDoChecks) const
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
vector< freeBound > fb
Definition iso.h:452
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long int numLevels_local
Definition iso.h:498
double RadRec_caseB
Definition iso.h:513
vector< two_photon > TwoNu
Definition iso.h:586
qList st
Definition iso.h:453
realnum CaseBCheck
Definition iso.h:510
t_conv conv
Definition conv.cpp:5
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
void GrainMakeDiffuse(void)
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_Heavy Heavy
Definition heavy.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
double emergent_line(double emissivity_in, double emissivity_out, long int ipCont)
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double MILNE_CONST
Definition physconst.h:233
UNUSED const double TE1RYD
Definition physconst.h:183
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void RT_iso_integrate_RRC(const long ipISO, const long nelem, const bool lgUpdateContinuum)
void RT_diffuse(void)
long int nSpecies
Definition taulines.cpp:21
vector< vector< TransitionList > > SatelliteLines
Definition taulines.cpp:38
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
multi_arr< int, 3 > ipSatelliteLines
Definition taulines.cpp:37
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
species * dBaseSpecies
Definition taulines.cpp:14
t_trace trace
Definition trace.cpp:5
void DumpLine(const TransitionProxy &t)
void CalcTwoPhotonEmission(two_photon &tnu, bool lgDoInduced)
void PrtTwoPhotonEmissCoef(const two_photon &tnu, const double &densityProduct)