cloudy trunk
Loading...
Searching...
No Matches
rt_line_one.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_line_escape do line radiative transfer,
4 * evaluates escape and destruction probability */
5/* NB no wind distinction - that is done where optical depths are incremented
6 * with line opacity - rt_line_one_tauinc and in line width for rad pressure */
7/*RT_line_fine_opacity do fine opacities for one line */
8/*RT_line_electron_scatter evaluate electron scattering escape probability */
9/*RT_line_pumping pumping by external and locally emitted radiation fields */
10#include "cddefines.h"
11#include "rfield.h"
12#include "doppvel.h"
13#include "dense.h"
14#include "opacity.h"
15#include "transition.h"
16#include "conv.h"
17#include "radius.h"
18#include "rt.h"
19#include "physconst.h"
20#include "cosmology.h"
21#include "thirdparty.h"
22#include "hydrogenic.h"
23
24/*RT_line_pumping pumping by external and locally emitted radiation fields */
26 /* the em line we will work on */
27 const TransitionProxy& t ,
28 /* this is option to not include line self shielding across this zone.
29 * this can cause pump to depend on zone thickness, and leads to unstable
30 * feedback in some models with the large H2 molecule, due to Solomon
31 * process depending on zone thickness and level populations. */
32 bool lgShield_this_zone,
33 realnum DopplerWidth)
34{
35 DEBUG_ENTRY( "RT_line_pumping()" );
36
37 ASSERT( t.ipCont() >= 1 );
38
39 /* pumping by incident and diffuse continua */
40 /* option to kill induced processes */
41 if( !rfield.lgInducProcess )
42 {
43 t.Emis().pump() = 0.;
44 }
45 else if( conv.lgFirstSweepThisZone || t.Emis().iRedisFun() == ipLY_A )
46 {
47 double dTau;
48 double shield_continuum;
49
50 /* continuum shielding into this function */
51 shield_continuum = RT_continuum_shield_fcn( t );
52
53 /* continuum upward pumping rate, A gu/gl abs prob occnum
54 * the "no induced" command causes continuum pumping to be set to 0
55 * this includes pumping by diffuse continuum */
56 t.Emis().pump() = t.Emis().Aul() * (*t.Hi()).g() / (*t.Lo()).g() * shield_continuum *(
57 rfield.OccNumbIncidCont[t.ipCont()-1] + rfield.OccNumbContEmitOut[t.ipCont()-1] );
58
59 dTau = (t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth +
60 opac.opacity_abs[t.ipCont()-1])* radius.drad_x_fillfac;
61
62 if( lgShield_this_zone && dTau > 1e-3 )
63 {
64 double dTauRel = dTau / (1.+t.Emis().TauCon());
65 /* correction for line optical depth across zone */
66 t.Emis().pump() *= log(1. + dTauRel ) / dTauRel;
67 }
68
69 /*
70 * This is an option to account for intrinsic absorption or emission line the lyman
71 * lines. This is done without changing the incident coarse continuum.
72 * Check if continuum pumping of H Lyman lines to be multiplied by scale factor
73 * hydro.xLymanPumpingScaleFactor is set with atom h-like Lyman pumping scale command
74 * Lya pump rate is always updated - this is simplest way to finese intrinsic absorption
75 * or emission
76 */
77 if ( t.ipLo() == 0 && t.systemIs(iso_sp[ipH_LIKE][ipHYDROGEN].tr) )
78 {
79 t.Emis().pump() *= hydro.xLymanPumpingScaleFactor;
80 }
81 /* NB NB line pumping by local diffuse emission is not included. */
82 }
83
84 return;
85}
86
87/*RT_line_electron_scatter evaluate electron scattering escape probability */
89 /* the em line we will work on */
90 const TransitionProxy& t ,
91 realnum DopplerWidth)
92{
93
94 DEBUG_ENTRY( "RT_line_electron_scatter()" );
95
96 /* escape following scattering off an electron */
97 /* this is turned off with the no scattering escape command */
98 if( !rt.lgElecScatEscape )
99 {
100 t.Emis().Pelec_esc() = 0.;
101 return;
102 }
103
104 // in the transition structure PopOpc is the pop relative to the ion,
105 // but has been converted to a physical population before this routine
106 // was called. No need to convert here */
107 double opac_line = (*t.Lo()).Pop() * t.Emis().opacity()/DopplerWidth;
108
109 /* the opacity in electron scattering */
110 double opac_electron = dense.eden*6.65e-25;
111 /* this is equation 5 of
112 *>>refer line desp Netzer, H., Elitzur, M., & Ferland, G. J. 1985, ApJ, 299, 752*/
113 double opacity_ratio = opac_electron/(opac_electron+opac_line);
114 /* keep total probability less than 0.1 */
115 t.Emis().Pelec_esc() = (realnum)opacity_ratio * MAX2(0.f,1.f-t.Emis().Pesc()-t.Emis().Pdest());
116
117 return;
118}
119
120/*RT_line_escape do line radiative transfer escape and destruction probabilities
121 * this routine sets */
123 /* the em line we will work on */
124 const TransitionProxy &t,
125 /* Stark escape probability to be added to Pesc */
126 realnum pestrk,
127 realnum DopplerWidth,
128 bool lgGoodTau)
129{
130 int nRedis = -1;
131
132 DEBUG_ENTRY( "RT_line_escape()" );
133
134 /* a runaway maser */
135 if( t.Emis().TauIn() < -30. )
136 {
137 fprintf( ioQQQ, "PROBLEM RT_line_escape called with large negative "
138 "optical depth, zone %.2f, setting lgAbort true.\n",
139 fnzone );
140 DumpLine(t);
141 /* return busted true instead of exit here, so that code will
142 * not hang under MPI */
143 lgAbort = true;
144 return;
145 }
146
147 if( cosmology.lgDo )
148 {
149 /* Sobolev escape */
150 if( conv.lgFirstSweepThisZone && lgGoodTau )
151 {
152 realnum tau_Sobolev = t.Emis().TauIn();
153
154 if( tau_Sobolev < 1E-5 )
155 {
156 t.Emis().Pesc() = 1.;
157 }
158 else
159 {
160 t.Emis().Pesc() = ( 1.f - exp( -1.f * tau_Sobolev ) )/ tau_Sobolev;
161 }
162
163 /* inward escaping fraction */
164 t.Emis().FracInwd() = rt.fracin;
165 }
166 fixit(); // is this correct?
167 nRedis = ipDEST_K2;
168 }
169 /* static solution - which type of line will determine
170 * which redistribution function */
171 /* iRedisFun() == 1 - alpha resonance line, partial redistribution,
172 * ipPRD == 1 */
173 else if( t.Emis().iRedisFun() == ipPRD )
174 {
175 /* incomplete redistribution with wings */
176 if( conv.lgFirstSweepThisZone && lgGoodTau )
177 {
178 t.Emis().Pesc() = (realnum)esc_PRD( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );
179
180 /* >>chng 03 jun 07, do not clobber esp prob when line is masing -
181 * this had effect of preventing total escape prob from getting larger than 1 */
182 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
183 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
184
185 /* inward escaping fraction */
186 t.Emis().FracInwd() = rt.fracin;
187 }
188 nRedis = ipDEST_INCOM;
189 }
190
191 /* complete redistribution without wings - t.ipLnRedis is ipCRD == -1 */
192 else if( t.Emis().iRedisFun() == ipCRD )
193 {
194 if( conv.lgFirstSweepThisZone && lgGoodTau )
195 {
196 /* >>chng 01 mar -6, escsub will call any of several esc prob routines,
197 * depending of how core is set. We always want core-only for this option,
198 * so call esca0k2(tau) directly */
199 t.Emis().Pesc() = (realnum)esc_CRDcore( t.Emis().TauIn(), t.Emis().TauTot() );
200
201 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
202 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
203
204 /* inward escaping fraction */
205 t.Emis().FracInwd() = rt.fracin;
206 }
207 nRedis = ipDEST_K2;
208 }
209
210 /* CRD with wings, = 2 */
211 else if( t.Emis().iRedisFun() == ipCRDW )
212 {
213 /* complete redistribution with damping wings */
214 if( conv.lgFirstSweepThisZone && lgGoodTau )
215 {
216 t.Emis().Pesc() = (realnum)esc_CRDwing( t.Emis().TauIn(), t.Emis().TauTot(), t.Emis().damp() );
217
218 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
219 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
220
221 /* inward escaping fraction */
222 t.Emis().FracInwd() = rt.fracin;
223 }
224 nRedis = ipDEST_K2;
225 }
226
227 /* Lya is special case */
228 else if( t.Emis().iRedisFun() == ipLY_A )
229 {
230 /* incomplete redistribution with wings, for special case of Lya
231 * uses fits to Hummer & Kunasz numerical results
232 * this routine is different because escape and dest probs
233 * are evaluated together, so no test of lgDoEsc */
234 if( lgGoodTau )
235 {
236 double dest , esin;
237
238 /* this will always evaluate escape prob, no matter what lgDoEsc is.
239 * The destruction prob comes back as dest */
240 t.Emis().Pesc() = (realnum)RTesc_lya( &esin, &dest, t.Emis().PopOpc(), t, DopplerWidth );
241
242 if( pestrk > 0.f && t.Emis().Pesc() < 1.f )
243 t.Emis().Pesc() = min( 1.f, t.Emis().Pesc() + pestrk );
244
245 /* this is current destruction rate */
246 t.Emis().Pdest() = (realnum)dest;
247
248 /* this is fraction of line which is inward escaping */
249 t.Emis().FracInwd() = rt.fracin;
250 }
251 }
252 else
253 {
254 fprintf( ioQQQ, " RT_line_escape called with impossible redistribution function %d\n",
255 t.Emis().iRedisFun());
256 ShowMe();
258 }
259
260 /* only do this if not Lya special case, since dest prob already done */
261 if( lgGoodTau && t.Emis().iRedisFun() != ipLY_A && t.Emis().opacity() > 0. )
262 {
263 t.Emis().Pdest() = (realnum)RT_DestProb(
264 /* the abundance of the species */
265 t.Emis().PopOpc(),
266 /* line center opacity in funny units (needs a vel) */
267 t.Emis().opacity(),
268 /* array index for line in continuum array,
269 * on f not c scale */
270 t.ipCont(),
271 /* line width in velocity units */
272 DopplerWidth,
273 /* escape probability */
274 t.Emis().Pesc(),
275 /* redistribution function */
276 nRedis);
277 }
278
279 return;
280}
281
282/*RT_line_fine_opacity do fine opacities for one line */
284 /* the em line we will work on */
285 const TransitionProxy& t ,
286 realnum DopplerWidth)
287
288{
289 DEBUG_ENTRY( "RT_line_fine_opacity()" );
290
291 /* this is line center frequency, including bulk motion of gas */
292 long int ipLineCenter = t.Emis().ipFine() + rfield.ipFineConVelShift;
293
294 /* define fine opacity fine grid fine mesh */
295 /* rfield.lgOpacityFine flag set false with no fine opacities command */
296 if( !conv.lgLastSweepThisZone || ipLineCenter < 0 || t.Emis().PopOpc() < SMALLFLOAT ||
297 ipLineCenter>rfield.nfine || !rfield.lgOpacityFine )
298 {
299 return;
300 }
301
302 /* number of fine opacity cells corresponding to one doppler width for current
303 * temperature, velocity field, and nuclear mass,
304 * rfield.fine_opac_velocity_width is width per cell, cm/s */
305 realnum cells_wide_1x = DopplerWidth/rfield.fine_opac_velocity_width;
306
307 /* line center opacity - type realnum since will add to fine opacity array,
308 * which is realnum */
309 realnum opac_line = (realnum)t.Emis().PopOpc() * t.Emis().opacity() / DopplerWidth;
310
311 // this is effective optical depth to this point. Do not do line if
312 // this product is less than SMALLFLOAT
313 double dTauEffec = opac_line*radius.depth_x_fillfac;
314 if( dTauEffec < SMALLFLOAT )
315 return;
316
317 /* core width of optically thick line, do 4x with exponential Doppler core,
318 * must be at least one cell, but profile is symmetric */
319 const bool doDamp = dTauEffec*t.Emis().damp()/9. > 0.1;
320 long int nCells_core = (long)(cells_wide_1x*4.f + 1.5f);
321 /* core is symmetric - make sure both upper and lower bounds
322 * are within continuum energy bin */
323 if( ipLineCenter - nCells_core < 1 )
324 nCells_core = ipLineCenter - 1;
325 if( ipLineCenter + nCells_core > rfield.nfine )
326 nCells_core = ipLineCenter -rfield.nfine - 1;
327
328 /* want core to be at least one cell wide */
329 nCells_core = MAX2( 1 , nCells_core );
330
331 long int nCells_damp;
332 /* include damping wings if optical depth is large */
333 if( doDamp )
334 {
335 // find number of cells to extend the damping wings - cells_wide_1x is one dop width
336 // tests with th85orion, stopping at half -h2 point, showed that 0.01 was
337 // needed for outer edge, given the definition of dTauEffec
338 realnum x = (realnum)sqrt( dTauEffec * t.Emis().damp()*100./SQRTPI ) * cells_wide_1x;
339 // test on size of x, which can exceed
340 // limits of long in extreme optical depths */
341 if( x<LONG_MAX )
342 {
343 nCells_damp = (long)x;
344
345 if( ipLineCenter-nCells_damp < 1 )
346 nCells_damp = ipLineCenter-1;
347
348 if( ipLineCenter+nCells_damp > rfield.nfine )
349 nCells_damp = rfield.nfine - ipLineCenter-1;
350 }
351 else
352 {
353 /* x was too big, just set to extreme range, which is
354 * number of cells to nearest boundary */
355 nCells_damp = MIN2( rfield.nfine-ipLineCenter , ipLineCenter )-1;
356 }
357 }
358 else
359 {
360 nCells_damp = nCells_core;
361 }
362
363 static vector<realnum> xprofile, profile;
364 xprofile.resize(nCells_damp);
365 profile.resize(nCells_damp);
366
367 for( long int i=0; i<nCells_damp; ++i )
368 {
369 /* distance from line center in units of doppler width */
370 xprofile[i] = (realnum) i/cells_wide_1x;
371 }
372
373 VoigtH(t.Emis().damp(), &xprofile[0], &profile[0], nCells_damp);
374
375 /* line center itself, must not double count here */
376 rfield.fine_opac_zone[ipLineCenter] += profile[0]*opac_line;
377 for( long int i=1; i<nCells_damp; ++i )
378 {
379 rfield.fine_opac_zone[ipLineCenter+i] += profile[i]*opac_line;
380 rfield.fine_opac_zone[ipLineCenter-i] += profile[i]*opac_line;
381 }
382
383 return;
384}
385
386/*RT_line_one do rt for emission line structure - calls RT_line_escape or RT_line_wind */
388 /* the em line we will work on */
389 const TransitionProxy &t,
390 /* this is option to not include line self shielding across this zone.
391 * this can cause pump to depend on zone thickness, and leads to unstable
392 * feedback in some models with the large H2 molecule, due to Solomon
393 * process depending on zone thickness and level populations. */
394 bool lgShield_this_zone,
395 /* Stark escape probability to be added to Pesc */
396 realnum pestrk,
397 realnum DopplerWidth )
398{
399 DEBUG_ENTRY( "RT_line_one()" );
400
401 // do nothing is population and this is not the very first call
402 // skip line transfer if requested with 'no line transfer' command, but never skip Lya
403 if( !rfield.lgDoLineTrans && (t.Emis().iRedisFun() != ipLY_A) )
404 {
405 return;
406 }
407
408 /* line damping constant at current temperature */
409 t.Emis().damp() = t.Emis().dampXvel() / DopplerWidth;
410 ASSERT( t.Emis().damp() > 0. );
411
412 // do not evaluate if no population
413 if( (*t.Lo()).Pop()<=SMALLFLOAT )
414 {
415 /* zero population, return after setting everything with side effects */
416 t.Emis().Pesc() = 1.f;
417
418 /* inward escaping fraction */
419 t.Emis().FracInwd() = 0.5;
420
421 /* pumping rate */
422 t.Emis().pump() = 0.;
423
424 /* destruction probability */
425 t.Emis().Pdest() = 0.;
426 t.Emis().Pelec_esc() = 0.;
427
428 return;
429 }
430
431 /* option to keep track of population values during calls,
432 * print out data to make histogram */
433 enum {DEBUG_LOC=false};
434 if( DEBUG_LOC )
435 {
436 static long int nTau[100];
437 long n;
438
439 if( nzone==0 )
440 {
441 for(n=0; n<100; ++n )
442 nTau[n] = 0;
443 }
444 if( (*t.Lo()).Pop()<=SMALLFLOAT )
445 n = 0;
446 else
447 n = (long)log10( (*t.Lo()).Pop() )+37;
448 n = MIN2( n , 99 );
449 n = MAX2( n , 0 );
450 ++nTau[n];
451 if( nzone > 183 )
452 {
453 for(n=0; n<100; ++n )
454 fprintf(ioQQQ,"%li\t%li\n", n , nTau[n] );
456 }
457 }
458
459 // transition is below plasma frequency - photons not emitted
460 if( t.EnergyErg() / EN1RYD <= rfield.plsfrq )
461 {
462 t.Emis().Pesc() = SMALLFLOAT;
463 t.Emis().Pdest() = SMALLFLOAT;
464 t.Emis().Pelec_esc() = SMALLFLOAT;
465 t.Emis().pump() = SMALLFLOAT;
466 }
467 else
468 {
469
470 /* this checks if we have overrun the optical depth scale,
471 * in which case the inward optical depth is greater than the
472 * previous iteration's total optical depth.
473 * We do not reevaluate escape probabilities if the optical depth
474 * scale has been overrun due to huge bogus change in solution
475 * that would result */
476 bool lgGoodTau = lgTauGood( t );
477
478 // the last sweep through this zone is to do the fine opacities
479 // the populations are not updated on the last sweep so the
480 // line transfer details also should not be updated
481 if( conv.lgLastSweepThisZone )
482 RT_line_fine_opacity( t , DopplerWidth );
483 else
484 {
485 RT_line_escape( t, pestrk, DopplerWidth , lgGoodTau);
486 RT_line_electron_scatter( t , DopplerWidth );
487 RT_line_pumping( t , lgShield_this_zone , DopplerWidth );
488 }
489 }
490
491 return;
492}
493
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipDEST_INCOM
Definition cddefines.h:300
#define MIN2
Definition cddefines.h:761
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipDEST_K2
Definition cddefines.h:298
const int ipCRDW
Definition cddefines.h:294
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipCRD
Definition cddefines.h:292
long min(int a, long b)
Definition cddefines.h:723
const int ipLY_A
Definition cddefines.h:296
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
const int ipPRD
Definition cddefines.h:290
realnum & damp() const
Definition emission.h:563
long int & ipFine() const
Definition emission.h:413
double & PopOpc() const
Definition emission.h:603
int & iRedisFun() const
Definition emission.h:403
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
realnum & TauIn() const
Definition emission.h:423
realnum & opacity() const
Definition emission.h:593
double & pump() const
Definition emission.h:473
realnum & FracInwd() const
Definition emission.h:463
realnum & Pdest() const
Definition emission.h:543
realnum & TauTot() const
Definition emission.h:433
realnum & dampXvel() const
Definition emission.h:553
int & ipLo() const
Definition transition.h:458
realnum EnergyErg() const
Definition transition.h:78
qList::iterator Lo() const
Definition transition.h:392
long & ipCont() const
Definition transition.h:450
bool systemIs(const TransitionList *query) const
Definition transition.h:343
qList::iterator Hi() const
Definition transition.h:396
EmissionList::reference Emis() const
Definition transition.h:408
t_conv conv
Definition conv.cpp:5
t_cosmology cosmology
Definition cosmology.cpp:11
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_hydro hydro
Definition hydrogenic.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double SQRTPI
Definition physconst.h:44
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
double esc_CRDwing(double tau, double tout, double damp)
double RT_continuum_shield_fcn(const TransitionProxy &t)
double esc_CRDcore(double tau, double tout)
double RT_DestProb(double abund, double crsec, long int ipanu, double widl, double escp, int nCore)
double RTesc_lya(double *esin, double *dest, double abund, const TransitionProxy &t, realnum DopplerWidth)
double esc_PRD(double tau, double tout, double damp)
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
STATIC void RT_line_pumping(const TransitionProxy &t, bool lgShield_this_zone, realnum DopplerWidth)
STATIC void RT_line_fine_opacity(const TransitionProxy &t, realnum DopplerWidth)
STATIC void RT_line_electron_scatter(const TransitionProxy &t, realnum DopplerWidth)
STATIC void RT_line_escape(const TransitionProxy &t, realnum pestrk, realnum DopplerWidth, bool lgGoodTau)
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:350
void DumpLine(const TransitionProxy &t)
bool lgTauGood(const TransitionProxy &t)
Definition transition.h:570