cloudy trunk
Loading...
Searching...
No Matches
rt_escprob.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/*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
4/*esc_PRD_1side fundamental escape probability radiative transfer routine for incomplete redistribution */
5/*RTesc_lya escape prob for hydrogen atom Lya, using Hummer and Kunasz results,
6 * called by hydropesc */
7/*esc_PRD escape probability radiative transfer for incomplete redistribution */
8/*esc_CRDwing escape probability for CRD with wings */
9/*esc_CRDcore escape probability for CRD with no wings */
10/*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
11/*RT_DestProb returns line destruction probability due to continuum opacity */
12/*RT_LineWidth determine half width of any line with known optical depths */
13#include "cddefines.h"
14#include "physconst.h"
15#include "dense.h"
16#include "conv.h"
17#include "rfield.h"
18#include "opacity.h"
19#include "lines_service.h"
20#include "taulines.h"
21#include "doppvel.h"
22#include "pressure.h"
23#include "wind.h"
24#include "rt.h"
25#include "iso.h"
26#include "thirdparty.h"
27
28inline double tau_from_here(double tau_tot, double tau_in)
29{
30 const double SCALE = 2.;
31 double tt = tau_tot - tau_in;
32 if (1)
33 {
34 /* help convergence by not letting tau go to zero at back edge of
35 * when there was a bad guess for the total optical depth
36 * note that this test is seldom hit since RTMakeStat does check
37 * for overrun */
38 if( tt < 0. )
39 {
40 tt = tau_in/SCALE;
41 }
42 }
43 else
44 {
45 // Alternatively, allow tau_from_here to go to zero, and then
46 // increase again. This will give more or less the right
47 // distribution of tau values, though with tau_out estimated to
48 // be zero somewhere inside the layer rather than at the edge.
49 //
50 // The iteration 1 inward-only treatment is equivalent to this,
51 // if the estimated tau_out is set to be zero.
52 //
53 // I think you'll find, I'm playing all of the right
54 // notes... but not necessarily in the right order.
55 //
56 // -- Eric Morecambe
57 tt = fabs(tt);
58 }
59 return tt;
60}
61
62/*escConE2 one of the forms of the continuum escape probability */
64{
65public:
67
68 double operator() (double x)
69 {
70 return exp(-chnukt_ContTkt*(x-1.))/x*e2(chnukt_ctau/POW3(x));
71 }
72};
73
74/* a continuum escape probability */
76{
77public:
79
80 double operator() (double x)
81 {
82 return exp(-chnukt_ContTkt*(x-1.))/x;
83 }
84};
85
86
87/*escmase escape probability for negative (masing) optical depths,*/
88STATIC double escmase(double tau);
89/*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
90STATIC void RTesc_lya_1side(double taume,
91 double beta,
92 realnum *esc,
93 realnum *dest,
94 /* position of line in frequency array on c scale */
95 long ipLine );
96
97double esc_PRD_1side(double tau,
98 double a)
99{
100 double atau,
101 b,
102 escinc_v;
103
104 DEBUG_ENTRY( "esc_PRD_1side()" );
105 ASSERT( a>0.0 );
106
107 /* this is one of the three fundamental escape probability routines
108 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
109 * it computes esc prob for incomplete redistribution
110 * */
111# if 0
112 if( strcmp(rt.chEscFunSubord,"SIMP") == 0 )
113 {
114 /* this set with "escape simple" command, used for debugging */
115 escinc_v = 1./(1. + tau);
116 return escinc_v;
117 }
118# endif
119
120 if( tau < 0. )
121 {
122 /* line mased */
123 escinc_v = escmase(tau);
124 }
125 else
126 {
127 /* first find coefficient b(tau) */
128 atau = a*tau;
129 if( atau > 1. )
130 {
131 b = 1.6 + (3.*pow(2.*a,-0.12))/(1. + atau);
132 }
133 else
134 {
135 double sqrtatau = sqrt(atau);
136 b = 1.6 + (3.*pow(2.*a,-0.12))*sqrtatau/(1. + sqrtatau);
137 }
138 b = MIN2(6.,b);
139
140 escinc_v = 1./(1. + b*tau);
141 }
142 return escinc_v;
143}
144
145// Implement equation (157) of Avrett & Loeser 1966
146// 1966SAOSR.201.....A for comparison with escape probability with
147// damping. Uses transformation x -> y/(1-y) to map domain of
148// integral to [0,1)
150{
152public:
154
156 {
157 if (y >= 1.)
158 return 0;
159 realnum x = y/(1.-y);
160 realnum phi;
161 VoigtU(damp,&x,&phi,1);
162 if (phi <= 0.)
163 {
164 return 0.;
165 }
166 else
167 {
168 return phi*expn(2,phi*tau)/POW2(1.-y);
169 }
170 }
171};
172
173template<class F>
174double trapezium(realnum xmin, realnum xmax, const F func)
175{
176 double tol = 1e-6;
177 vector<double> vxmin, vxmax, fxmin, fxmax, step;
178 vxmin.push_back(xmin);
179 vxmax.push_back(xmax);
180 fxmin.push_back(func(xmin));
181 fxmax.push_back(func(xmax));
182 step.push_back(0.5*(xmax-xmin)*(fxmin[0]+fxmax[0]));
183 for (long level = 0; level < 25; ++level)
184 {
185 size_t nstep = step.size();
186 double current = 0.0;
187 for (size_t i=0; i<nstep; ++i)
188 current += step[i];
189 for (size_t i=0; i<nstep; ++i)
190 {
191 if (vxmax[i] > vxmin[i])
192 {
193 if (level <= 3 || step[i] > tol*current
194 || (i == nstep-1 && current <= 0.0) )
195 {
196 double xbar=0.5*(vxmin[i]+vxmax[i]);
197 vxmin.push_back(xbar);
198 fxmin.push_back(func(xbar));
199 vxmax.push_back(vxmax[i]);
200 fxmax.push_back(fxmax[i]);
201 long ilast = vxmax.size()-1;
202 vxmax[i] = xbar;
203 fxmax[i] = fxmin[ilast];
204 step[i] = 0.5*(vxmax[i]-vxmin[i])*(fxmin[i]+fxmax[i]);
205 //fprintf(ioQQQ,"%g %g\n",vxmin[ilast],fxmin[ilast],
206 // step[i]);
207 step.push_back(0.5*(vxmax[ilast]-vxmin[ilast])*
208 (fxmin[ilast]+fxmax[ilast]));
209 }
210 }
211 }
212 }
213 double current = 0.0;
214 for (size_t i=0; i<step.size(); ++i)
215 current += step[i];
216 return current;
217}
218
219
220/*esc_CRDwing_1side fundamental escape probability radiative transfer routine, for complete redistribution */
221double esc_CRDwing_1side(double tau,
222 double a )
223{
224 DEBUG_ENTRY( "esc_CRDwing_1side()" );
225
226 /* this is one of the three fundamental escape probability routines
227 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
228 * it computes esc prob for complete redistribution with wings
229 * computes escape prob for complete redistribution in one direction
230 * */
231
232 /* this is the only case that this routine computes,
233 * and is the usual case for subordinate lines,
234 * complete redistribution with damping wings */
235
236 if (0)
237 {
238 // try to compare the formulae from this function with A&L
239 // exact expression
240 a = 1000.;
241 double taustep = 2., taumin=1e-8,taumax=1e14;
243 for (tau = taumin; tau < taumax; tau *= taustep)
244 {
245 double esccom_v = esca0k2(tau);
246 double sqrta = sqrt(a);
247 double pwing = tau > 0.0 ? sqrta/(sqrta+0.5*3.0*sqrt(SQRTPI*tau)) :
248 1.0;
249 double esctot = esccom_v*(1.0-pwing)+pwing;
250 double intgral = IntDamp.sum(0.,1.,k2DampArg(a,SQRTPI*tau));
251 double intgralt = trapezium(0.,1.,k2DampArg(a,SQRTPI*tau));
252 fprintf(ioQQQ,"cfesc %15.8g %15.8g %15.8g %15.8g %15.8g %15.8g\n",
253 tau,a,esccom_v,esctot,2.0*intgral,2.0*intgralt);
254 }
255 exit(0);
256 }
257 double esccom_v = esca0k2(tau);
258
259 // Escape probability correction for finite damping
260 // Results agree to +/- 20% from a=1e-3->1e3, no change for a->0
261
262 double sqrta = sqrt(a);
263 double scal = a*(1.0+a+tau)/(POW2(1.0+a)+a*tau);
264 double pwing = scal*((tau > 0.0) ? sqrta/sqrt(a+2.25*SQRTPI*tau) : 1.0);
265 return esccom_v*(1.0-pwing)+pwing;
266}
267
268/*RTesc_lya escape prob for hydrogen atom Lya, using
269 >>refer La escp Hummer, D.G., & Kunasz, P.B., 1980, ApJ, 236, 609
270 * called by hydropesc, return value is escape probability */
272 /* the inward escape probability */
273 double *esin,
274 /* the destruction probility */
275 double *dest,
276 /* abundance of the species */
277 double abund,
278 const TransitionProxy& t,
279 realnum DopplerWidth)
280{
281 double beta,
282 conopc,
283 escla_v;
284 realnum dstin,
285 dstout;
286
287 DEBUG_ENTRY( "RTesc_lya()" );
288
289 /*
290 * this is one of the three fundamental escape probability functions
291 * the three are esc_CRDwing_1side, esc_PRD_1side, and RTesc_lya
292 * evaluate esc prob for LA
293 * optical depth in outer direction always defined
294 */
295
296 if( t.Emis().TauTot() - t.Emis().TauIn() < 0. )
297 {
298 /* this is the case if we overrun the optical depth scale
299 * just leave things as they are */
300 escla_v = t.Emis().Pesc();
301 rt.fracin = t.Emis().FracInwd();
302 *esin = rt.fracin;
303 *dest = t.Emis().Pdest();
304 return escla_v;
305 }
306
307 /* incomplete redistribution */
308 conopc = opac.opacity_abs[ t.ipCont()-1 ];
309 if( abund > 0. )
310 {
311 /* the continuous opacity is positive, we have a valid soln */
312 beta = conopc/(abund/SQRTPI*t.Emis().opacity()/
313 DopplerWidth + conopc);
314 }
315 else
316 {
317 /* abundance is zero, set miniumum dest prob */
318 beta = 1e-10;
319 }
320
321 /* find rt.wayin, the escape prob in inward direction */
323 t.Emis().TauIn(),
324 beta,
325 &rt.wayin,
326 &dstin,
327 /* position of line in energy array on C scale */
328 t.ipCont()-1);
329
330 ASSERT( (rt.wayin <= 1.) && (rt.wayin >= 0.) && (dstin <= 1.) && (dstin >= 0.) );
331
332 /* find rt.wayout, the escape prob in outward direction */
333 RTesc_lya_1side(MAX2(t.Emis().TauTot()/100.,
334 t.Emis().TauTot()-t.Emis().TauIn()),
335 beta,
336 &rt.wayout,
337 &dstout,
338 t.ipCont()-1);
339
340 ASSERT( (rt.wayout <= 1.) && (rt.wayout >= 0.) && (dstout <= 1.) && (dstout >= 0.) );
341
342 /* esc prob is mean of in and out */
343 escla_v = (rt.wayin + rt.wayout)/2.;
344 /* the inward escaping part of the line */
345 *esin = rt.wayin;
346
347 /* dest prob is mean of in and out */
348 *dest = (dstin + dstout)/2.f;
349 /* >>chng 02 oct 02, sum of escape and dest prob must be less then unity,
350 * for very thin models this forces dest prob to go to zero,
351 * rather than the value of DEST0, caught by Jon Slavin */
352 *dest = (realnum)MIN2( *dest , 1.-escla_v );
353 /* but dest prob can't be negative */
354 *dest = (realnum)MAX2(0., *dest );
355
356 /* fraction of line emitted in inward direction */
357 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
358 ASSERT( escla_v >=0. && *dest>=0. && *esin>=0. );
359 return escla_v;
360}
361
362/*esc_PRD escape probability radiative transfer for incomplete redistribution */
363double esc_PRD(double tau,
364 double tau_out,
365 double damp )
366{
367 double escgrd_v,
368 tt;
369
370 DEBUG_ENTRY( "esc_PRD()" );
371
372 ASSERT( damp > 0. );
373
374 /* find escape prob for incomp redis, average of two 1-sided probs*/
375
376 if( iteration > 1 )
377 {
378 tt = tau_from_here(tau_out, tau);
379
380 rt.wayin = (realnum)esc_PRD_1side(tau,damp);
381 rt.wayout = (realnum)esc_PRD_1side(tt,damp);
382 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
383 escgrd_v = 0.5*(rt.wayin + rt.wayout);
384 }
385 else
386 {
387 /* outward optical depth not defined, dont estimate fraction out */
388 rt.fracin = 0.;
389 rt.wayout = 1.;
390 escgrd_v = esc_PRD_1side(tau,damp);
391 rt.wayin = (realnum)escgrd_v;
392 }
393
394 ASSERT( escgrd_v > 0. );
395 return escgrd_v;
396}
397
398/*esc_CRDwing escape probability radiative transfer for CRDS in core only */
399double esc_CRDwing(double tau_in,
400 double tau_out,
401 double damp)
402{
403 double escgrd_v,
404 tt;
405
406 DEBUG_ENTRY( "esc_CRDwing()" );
407
408 /* find escape prob for CRD with damping wings, average of two 1-sided probs*/
409
410 /* crd with wings */
411 if( iteration > 1 )
412 {
413 /* outward optical depth if defined */
414 /* >>chng 03 jun 07, add test for masers here */
415 if( tau_out <0 || tau_in < 0. )
416 {
417 /* we have a maser, use smallest optical depth to damp it out */
418 tt = MIN2( tau_out , tau_in );
419 tau_in = tt;
420 }
421 else
422 {
423 tt = tau_from_here(tau_out, tau_in);
424 }
425
426 rt.wayin = (realnum)esc_CRDwing_1side(tau_in,damp);
427 rt.wayout = (realnum)esc_CRDwing_1side(tt,damp);
428 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
429 escgrd_v = 0.5*(rt.wayin + rt.wayout);
430 }
431 else
432 {
433 /* outward optical depth not defined, dont estimate fraction out */
434 rt.fracin = 0.;
435 rt.wayout = 1.;
436 escgrd_v = esc_CRDwing_1side(tau_in,damp);
437 rt.wayin = (realnum)escgrd_v;
438 }
439
440 ASSERT( escgrd_v > 0. );
441 return escgrd_v;
442}
443
444/*esc_CRDwing escape probability radiative transfer for incomplete redistribution */
445double esc_CRDcore(double tau_in,
446 double tau_out)
447{
448 double escgrd_v,
449 tt;
450
451 DEBUG_ENTRY( "esc_CRDcore()" );
452
453 /* find escape prob for CRD with damping wings, average of two 1-sided probs*/
454
455 /* crd with wings */
456 if( iteration > 1 )
457 {
458 /* outward optical depth if defined */
459 /* >>chng 03 jun 07, add test for masers here */
460 if( tau_out <0 || tau_in < 0. )
461 {
462 /* we have a maser, use smallest optical depth to damp it out */
463 tt = MIN2( tau_out , tau_in );
464 tau_in = tt;
465 }
466 else
467 {
468 tt = tau_from_here(tau_out, tau_in);
469 }
470
471 rt.wayin = (realnum)esca0k2(tau_in);
472 rt.wayout = (realnum)esca0k2(tt);
473 rt.fracin = rt.wayin/(rt.wayin + rt.wayout);
474 escgrd_v = 0.5*(rt.wayin + rt.wayout);
475 }
476 else
477 {
478 /* outward optical depth not defined, dont estimate fraction out */
479 rt.fracin = 0.;
480 rt.wayout = 1.;
481 escgrd_v = esca0k2(tau_in);
482 rt.wayin = (realnum)escgrd_v;
483 }
484
485 ASSERT( escgrd_v > 0. );
486 return escgrd_v;
487}
488
489/*esca0k2 derive Hummer's K2 escape probability for Doppler core only */
490double esca0k2(double taume)
491{
492 double arg,
493 esca0k2_v,
494 suma,
495 sumb,
496 sumc,
497 sumd,
498 tau;
499 static const double a[5]={1.00,-0.1117897,-0.1249099917,-9.136358767e-3,
500 -3.370280896e-4};
501 static const double b[6]={1.00,0.1566124168,9.013261660e-3,1.908481163e-4,
502 -1.547417750e-7,-6.657439727e-9};
503 static const double c[5]={1.000,19.15049608,100.7986843,129.5307533,-31.43372468};
504 static const double d[6]={1.00,19.68910391,110.2576321,169.4911399,-16.69969409,
505 -36.664480000};
506
507 DEBUG_ENTRY( "esca0k2()" );
508
509 /* compute Hummer's K2 escape probability function for a=0
510 * using approx from
511 * >>refer line escp Hummer, D.G., xxxx, JQRST, 26, 187.
512 *
513 * convert to David's opacity */
514 tau = taume*SQRTPI;
515
516 if( tau < 0. )
517 {
518 /* the line mased */
519 esca0k2_v = escmase(taume);
520
521 }
522 else if( tau < 0.01 )
523 {
524 esca0k2_v = 1. - 2.*tau;
525
526 }
527 else if( tau <= 11. )
528 {
529 suma = a[0] + tau*(a[1] + tau*(a[2] + tau*(a[3] + a[4]*tau)));
530 sumb = b[0] + tau*(b[1] + tau*(b[2] + tau*(b[3] + tau*(b[4] +
531 b[5]*tau))));
532 esca0k2_v = tau/2.5066283*log(tau/SQRTPI) + suma/sumb;
533
534 }
535 else
536 {
537 /* large optical depth limit */
538 arg = 1./log(tau/SQRTPI);
539 sumc = c[0] + arg*(c[1] + arg*(c[2] + arg*(c[3] + c[4]*arg)));
540 sumd = d[0] + arg*(d[1] + arg*(d[2] + arg*(d[3] + arg*(d[4] +
541 d[5]*arg))));
542 esca0k2_v = (sumc/sumd)/(2.*tau*sqrt(log(tau/SQRTPI)));
543 }
544 return esca0k2_v;
545}
546
547/*escmase escape probability for negative (masing) optical depths */
548STATIC void FindNeg( void )
549{
550 long int i;
551
552 DEBUG_ENTRY( "FindNeg()" );
553
554 /* do the level 1 lines */
555 for( i=1; i <= nLevel1; i++ )
556 {
557 /* check if a line was a strong maser */
558 if( TauLines[i].Emis().TauIn() < -1. )
559 DumpLine(TauLines[i]);
560 }
561
562 /* Generic atoms & molecules from databases
563 * added by Humeshkar Nemala*/
564 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
565 {
566 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
567 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
568 {
569 if((*em).TauIn() < -1. )
570 DumpLine((*em).Tran());
571 }
572 }
573
574 /* now do the level 2 lines */
575 for( i=0; i < nWindLine; i++ )
576 {
577 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
578 {
579 /* check if a line was a strong maser */
580 if( TauLine2[i].Emis().TauIn() < -1. )
581 DumpLine(TauLine2[i]);
582 }
583 }
584
585 /* now do the hyperfine structure lines */
586 for( i=0; i < nHFLines; i++ )
587 {
588 /* check if a line was a strong maser */
589 if( HFLines[i].Emis().TauIn() < -1. )
590 DumpLine(HFLines[i]);
591 }
592
593 return;
594}
595
596STATIC double escmase(double tau)
597{
598 double escmase_v;
599
600 DEBUG_ENTRY( "escmase()" );
601
602 /* this is the only routine that computes maser escape probabilities */
603 ASSERT( tau <= 0. );
604
605 if( tau > -0.1 )
606 {
607 escmase_v = 1. - tau*(0.5 + tau/6.);
608 }
609 else if( tau > -30. )
610 {
611 escmase_v = (1. - exp(-tau))/tau;
612 }
613 else
614 {
615 fprintf( ioQQQ, " DISASTER escmase called with 2big tau%10.2e\n",
616 tau );
617 fprintf( ioQQQ, " This is zone number%4ld\n", nzone );
618 FindNeg();
619 ShowMe();
621 }
622
623 ASSERT( escmase_v >= 1. );
624 return escmase_v;
625}
626
627/*esccon continuum escape probability */
628double esccon(double tau,
629 double hnukt)
630{
631 double dinc,
632 escpcn_v,
633 sumesc,
634 sumrec;
635
636 DEBUG_ENTRY( "esccon()" );
637
638 /* computes continuum escape probabilities */
639 if( tau < 0.01 )
640 {
641 escpcn_v = 1.;
642 return escpcn_v;
643 }
644
645 else if( hnukt > 1. && tau > 100. )
646 {
647 escpcn_v = 1e-20;
648 return escpcn_v;
649 }
650
651 my_Integrand_conrec func_conrec;
652 func_conrec.chnukt_ContTkt = hnukt;
654
655 my_Integrand_escConE2 func_escConE2;
656 func_escConE2.chnukt_ContTkt = hnukt;
657 func_escConE2.chnukt_ctau = tau;
659
660 dinc = 10./hnukt;
661 sumrec = conrec.sum(1.,1.+dinc,func_conrec);
662 sumesc = escConE2.sum(1.,1.+dinc,func_escConE2);
663
664 if( sumrec > 0. )
665 {
666 escpcn_v = sumesc/sumrec;
667 }
668 else
669 {
670 escpcn_v = 0.;
671 }
672 return escpcn_v;
673}
674
675/*RTesc_lya_1side fit Hummer and Kunasz escape probability for hydrogen atom Lya */
676STATIC void RTesc_lya_1side(double taume,
677 double beta,
678 realnum *esc,
679 realnum *dest,
680 /* position of line in frequency array on c scale */
681 long ipLine )
682{
683 double esc0,
684 fac,
685 fac1,
686 fac2,
687 tau,
688 taucon,
689 taulog;
690
691 /* DEST0 is the smallest destruction probability to return
692 * in high metallicity models, in rt.h
693 const double DEST0=1e-8;*/
694
695 DEBUG_ENTRY( "RTesc_lya_1side()" );
696
697 /* fits to numerical results of Hummer and Kunasz Ap.J. 80 */
698 tau = taume*SQRTPI;
699
700 /* this is the real escape probability */
701 esc0 = 1./((0.6451 + tau)*(0.47 + 1.08/(1. + 7.3e-6*tau)));
702
703 esc0 = MAX2(0.,esc0);
704 esc0 = MIN2(1.,esc0);
705
706 if( tau > 0. )
707 {
708 taulog = log10(MIN2(1e8,tau));
709 }
710 else
711 {
712 /* the line mased
713 *>>chng 06 sep 08, kill xLaMase
714 hydro.xLaMase = MIN2(hydro.xLaMase,(realnum)tau);*/
715 taulog = 0.;
716 *dest = 0.;
717 *esc = (realnum)esc0;
718 }
719
720 if( beta > 0. )
721 {
722 taucon = MIN2(2.,beta*tau);
723
724 if( taucon > 1e-3 )
725 {
726 fac1 = -1.25 + 0.475*taulog;
727 fac2 = -0.485 + 0.1615*taulog;
728 fac = -fac1*taucon + fac2*POW2(taucon);
729 fac = pow(10.,fac);
730 fac = MIN2(1.,fac);
731 }
732 else
733 {
734 fac = 1.;
735 }
736
737 *esc = (realnum)(esc0*fac);
738 /* MIN puts cat at 50 */
739 *dest = (realnum)(beta/(0.30972 - MIN2(.28972,0.03541667*taulog)));
740 }
741
742 else
743 {
744 *dest = 0.;
745 *esc = (realnum)esc0;
746 }
747
748 *dest = MIN2(*dest,1.f-*esc);
749 *dest = MAX2(0.f,*dest);
750
751 /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
752 * in this case scattering is much more likely than absorption on this event */
753 *dest = (realnum)( (1. - opac.albedo[ipLine]) * *dest + opac.albedo[ipLine]*DEST0);
754 /* this is for debugging H Lya */
755 {
756 /*@-redef@*/
757 enum {BUG=false};
758 /*@+redef@*/
759 if( BUG )
760 {
761 fprintf(ioQQQ,"scatdest tau %.2e beta%.2e 1-al%.2e al%.2e dest%.2e \n",
762 taume,
763 beta,
764 (1. - opac.albedo[ipLine]),
765 opac.albedo[ipLine] ,
766 *dest
767 );
768 }
769 }
770 return;
771}
772
773/*RT_DestProb returns line destruction probability due to continuum opacity */
775 /* abundance of species */
776 double abund,
777 /* its line absorption cross section */
778 double crsec,
779 /* pointer to energy within continuum array, to get background opacity,
780 * this is on the f not c scale */
781 long int ipanu,
782 /* line width */
783 double widl,
784 /* escape probability */
785 double escp,
786 /* type of redistribution function */
787 int nCore)
788{
789 /* this will be the value we shall return */
790 double eovrlp_v;
791
792 double conopc,
793 beta;
794
795 /* DEST0 is the smallest destruction probability to return
796 * in high metallicity models
797 * this was set to 1e-8 until 99nov18,
798 * in cooling flow model the actual Lya ots dest prob was 1e-16,
799 * and this lower limit of 1e-8 caused energy balance problems,
800 * since dest prob was 8 orders of magnitude too great.
801 * >>chng 99 nov 18, to 1e-20, but beware that comments indicate that
802 * this will cause problems with high metallicity clouds(?) */
803 /* >>chng 00 jun 04, to 0 since very feeble ionization clouds, with almost zero opacity,
804 * this was a LARGE number */
805 /*const double DEST0=1e-20;
806 const double DEST0=0.;*/
807
808 DEBUG_ENTRY( "RT_DestProb()" );
809
810 /* computes "escape probability" due to continuum destruction of
811 *
812 * if esc prob gt 1 then line is masing - return small number for dest prob */
813 /* >>>chng 99 apr 10, return min dest when scattering greater than abs */
814 /* no idea of opacity whatsoever, on very first soln for this model */
815 /* >>chng 05 mar 20, add test on line being above upper bound of frequency
816 * do not want to evaluate opacity in this case since off scale */
817 if( escp >= 1.0 || !conv.nTotalIoniz || ipanu >= rfield.nflux )
818 {
819 eovrlp_v = 0.;
820 return eovrlp_v;
821 }
822
823 /* find continuum opacity */
824 conopc = opac.opacity_abs[ipanu-1];
825
826 ASSERT( crsec > 0. );
827
828 /* may be no population, cannot use above test since return 0 not DEST0 */
829 if( abund <= 0. || conopc <= 0. )
830 {
831 /* do not set this to DEST0 since energy not then conserved */
832 eovrlp_v = 0.;
833 return eovrlp_v;
834 }
835
836 /* fac of 1.7 convert to Hummer convention for line opacity */
837 beta = conopc/(abund*SQRTPI*crsec/widl + conopc);
838 /* >>chng 04 may 10, rm * 1-pesc)
839 beta = MIN2(beta,(1.-escp)); */
840
841 if( nCore == ipDEST_INCOM )
842 {
843 /* fits to
844 * >>>refer la esc Hummer and Kunasz 1980 Ap.J. 236,609.
845 * the max value of 1e-3 is so that we do not go too far
846 * beyond what Hummer and Kunasz did, discussed in
847 * >>refer rt esc proc Ferland, G.J., 1999, ApJ, 512, 247 */
850 eovrlp_v = MIN2(1e-3,8.5*beta);/**/
851 }
852 else if( nCore == ipDEST_K2 )
853 {
854 /* Doppler core only; a=0., Hummer 68
855 eovrlp_v = RT_DestHummer(beta);*/
856 eovrlp_v = MIN2(1e-3,8.5*beta);/**/
857 }
858 else if( nCore == ipDEST_SIMPL )
859 {
860 /* this for debugging only
861 eovrlp_v = 8.5*beta;*/
862 /* >>chng 04 may 13, use same min function */
863 eovrlp_v = MIN2(1e-3,8.5*beta);/**/
864 }
865 else
866 {
867 fprintf( ioQQQ, " chCore of %i not understood by RT_DestProb.\n",
868 nCore );
870 }
871
872 /* renorm to unity */
873 eovrlp_v /= 1. + eovrlp_v;
874
875 /* multiply by 1-escape prob, since no destruction when optically thin */
876 eovrlp_v *= 1. - escp;
877
878 /*check results in bounds */
879 ASSERT( eovrlp_v >= 0. );
880 ASSERT( eovrlp_v <= 1. );
881
882 {
883 /* debugging code for Lya problems */
884 /*@-redef@*/
885 enum {DEBUG_LOC=false};
886 /*@+redef@*/
887 if( DEBUG_LOC )
888 {
889 if( rfield.anu[ipanu-1]>0.73 && rfield.anu[ipanu-1]<0.76 &&
890 fp_equal( abund, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().PopOpc() ) )
891 {
892 fprintf(ioQQQ,"%li RT_DestProb\t%g\n",
893 nzone, eovrlp_v );
894 }
895 }
896 }
897
898 /* >>chng 04 may 10, rm min */
899 /* this hack removed since no fundamental reason for it to be here,
900 * this should be added to scattering escape, if included at all */
901# if 0
902 /* >>chng 99 apr 12, limit destruction prob in case where gas dominated by scattering.
903 * in this case scattering is much more likely than absorption on this event
904 eovrlp_v = (1. - opac.albedo[ipanu-1]) * eovrlp_v +
905 opac.albedo[ipanu-1]*DEST0; */
906 /* >>chng 01 aug 11, add factor of 3 for increase in mean free path, and min on 0 */
907 /*eovrlp_v = MAX2(DEST0,1. - 3.*opac.albedo[ipanu-1]) * eovrlp_v +
908 opac.albedo[ipanu-1]*DEST0;*/
909 eovrlp_v = POW2(1. - opac.albedo[ipanu-1]) * eovrlp_v +
910 opac.albedo[ipanu-1]*0.;
911# endif
912
913 return eovrlp_v;
914}
915
916/*RT_LineWidth compute line width (cm/sec), using optical depth array information
917 * this is where effects of wind are done */
918double RT_LineWidth(const TransitionProxy& t, realnum DopplerWidth)
919{
920 double RT_LineWidth_v,
921 aa,
922 atau,
923 b,
924 r,
925 vth;
926 realnum tau;
927
928 DEBUG_ENTRY( "RT_LineWidth()" );
929
930 /* uses line width from
931 * >>refer esc prob Bonilha et al. Ap.J. (1979) 233 649
932 * return value is half velocity width*(1-ESC PROB) [cm s-1]
933 * this assumes incomplete redistribution, damp.tau^1/3 width */
934
935 /* thermal width */
936 vth = DopplerWidth;
937
938 /* optical depth in outer direction is defined
939 * on second and later iterations.
940 * smaller of inner and outer optical depths is chosen for esc prob */
941 if( iteration > 1 )
942 {
943 /* optical depth to shielded face */
944 realnum tauout = t.Emis().TauTot() - t.Emis().TauIn();
945
946 /* >>chng 99 apr 22 use smaller optical depth */
947 tau = MIN2( t.Emis().TauIn() , tauout );
948 }
949 else
950 {
951 tau = t.Emis().TauIn();
952 }
953 /* do not evaluate line width if quite optically thin - will be dominated
954 * by noise in this case */
955 if( tau <1e-3 )
956 return 0;
957
958 t.Emis().damp() = t.Emis().dampXvel() / DopplerWidth;
959 ASSERT( t.Emis().damp() > 0. );
960
961 double Pesc = esc_PRD_1side( tau , t.Emis().damp());
962
963 /* max optical depth is thermalization length */
964 realnum therm = (realnum)(5.3e16/MAX2(1e-15,dense.eden));
965 if( tau > therm )
966 {
967 /* \todo 2 this seems to create an inconsistency as it changes tau
968 * for the purposes of this routine (to return the line-width),
969 * but this leaves the actual optical depth unchanged. */
970 pressure.lgPradDen = true;
971 tau = therm;
972 }
973
974 /* >>chng 01 jun 23, use wind vel instead of rt since rt deleted */
975 /* >>chng 04 may 13, use thermal for subsonic cases */
980 if( ! wind.lgBallistic() )
981 {
982 /* static geometry */
983 /* esc prob has noise if smaller than FLT_EPSILON, or is masing */
984 if( (tau-opac.taumin)/100. < FLT_EPSILON )
985 {
986 RT_LineWidth_v = 0.;
987 }
988 else if( tau <= 20. )
989 {
990 atau = -6.907755;
991 if( tau > 1e-3 )
992 atau = log(tau);
993 aa = 4.8 + 5.2*tau + (4.*tau - 1.)*atau;
994 b = 6.5*tau - atau;
995 double escProb = Pesc + t.Emis().Pelec_esc() + t.Emis().Pdest();
996 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1., escProb ) );
997 /* small number roundoff can dominate this process */
998 if( escProb >= 1. - 100. * FLT_EPSILON )
999 RT_LineWidth_v = 0.;
1000 }
1001 else
1002 {
1003 ASSERT( t.Emis().damp()*tau >= 0.);
1004 atau = log(MAX2(0.0001,tau));
1005 aa = 1. + 2.*atau/pow(1. + 0.3*t.Emis().damp()*tau,0.6667) + pow(6.5*
1006 t.Emis().damp()*tau,0.333);
1007 b = 1.6 + 1.5/(1. + 0.20*t.Emis().damp()*tau);
1008 RT_LineWidth_v = vth*0.8862*aa/b*(1. - MIN2( 1. ,
1009 (Pesc+ t.Emis().Pelec_esc() + t.Emis().Pdest())) );
1010 }
1011
1012 /* we want full width, not half width */
1013 RT_LineWidth_v *= 2.;
1014
1015 }
1016 else
1017 {
1018 /* ballistic wind */
1019 r = t.Emis().damp()*tau/PI;
1020 if( r <= 1. )
1021 {
1022 RT_LineWidth_v = vth*sqrt(log(MAX2(1.,tau))*.2821);
1023 }
1024 else
1025 {
1026 RT_LineWidth_v = 2.*fabs(wind.windv0);
1027 if( r*vth <= RT_LineWidth_v )
1028 {
1029 RT_LineWidth_v = vth*r*log(RT_LineWidth_v/(r*vth));
1030 }
1031 }
1032 }
1033
1034 ASSERT( RT_LineWidth_v >= 0. );
1035 return RT_LineWidth_v;
1036}
1037
1038/*RT_DestHummer evaluate Hummer's betaF(beta) function */
1039double RT_DestHummer(double beta) /* beta is ratio of continuum to mean line opacity,
1040 * returns dest prob = beta F(beta) */
1041{
1042 double fhummr_v,
1043 x;
1044
1045 DEBUG_ENTRY( "RT_DestHummer()" );
1046
1047 /* evaluates Hummer's F(beta) function for case where damping
1048 * constant is zero, are returns beta*F(beta)
1049 * fit to Table 1, page 80, of Hummer MNRAS 138, 73-108.
1050 * beta is ratio of continuum to line opacity; FUMMER is
1051 * product of his F() times beta; the total destruction prob
1052 * this beta is Hummer's normalization of the Voigt function */
1053
1054 ASSERT( beta >= 0.);/* non-positive is unphysical */
1055 if( beta <= 0. )
1056 {
1057 fhummr_v = 0.;
1058 }
1059 else
1060 {
1061 x = log10(beta);
1062 if( x < -5.5 )
1063 {
1064 fhummr_v = 3.8363 - 0.56329*x;
1065 }
1066 else if( x < -3.5 )
1067 {
1068 fhummr_v = 2.79153 - 0.75325*x;
1069 }
1070 else if( x < -2. )
1071 {
1072 fhummr_v = 1.8446 - 1.0238*x;
1073 }
1074 else
1075 {
1076 fhummr_v = 0.72500 - 1.5836*x;
1077 }
1078 fhummr_v *= beta;
1079 }
1080 return fhummr_v;
1081}
t_abund abund
Definition abund.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
const int ipDEST_INCOM
Definition cddefines.h:300
#define POW3
Definition cddefines.h:936
#define MIN2
Definition cddefines.h:761
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipDEST_K2
Definition cddefines.h:298
const int ipDEST_SIMPL
Definition cddefines.h:302
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
double e2(double t)
Definition service.cpp:299
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
long nWindLine
Definition cdinit.cpp:19
EmissionProxy::iterator iterator
Definition emission.h:317
realnum & damp() const
Definition emission.h:563
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & TauIn() const
Definition emission.h:423
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
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
long & ipCont() const
Definition transition.h:450
EmissionList::reference Emis() const
Definition transition.h:408
realnum tau
realnum operator()(realnum y) const
k2DampArg(realnum damp, realnum tau)
realnum damp
double operator()(double x)
double operator()(double x)
t_conv conv
Definition conv.cpp:5
t_dense dense
Definition dense.cpp:24
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
UNUSED const double PI
Definition physconst.h:29
UNUSED const double SQRTPI
Definition physconst.h:44
t_pressure pressure
Definition pressure.cpp:5
static long int * ipLine
t_rfield rfield
Definition rfield.cpp:8
t_rt rt
Definition rt.cpp:5
#define DEST0
Definition rt.h:227
double RT_DestHummer(double beta)
double tau_from_here(double tau_tot, double tau_in)
double esc_PRD_1side(double tau, double a)
double esc_CRDcore(double tau_in, double tau_out)
double trapezium(realnum xmin, realnum xmax, const F func)
double esca0k2(double taume)
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
STATIC void FindNeg(void)
double RT_DestProb(double abund, double crsec, long int ipanu, double widl, double escp, int nCore)
STATIC double escmase(double tau)
double esc_CRDwing_1side(double tau, double a)
double esc_CRDwing(double tau_in, double tau_out, double damp)
double RTesc_lya(double *esin, double *dest, double abund, const TransitionProxy &t, realnum DopplerWidth)
double esc_PRD(double tau, double tau_out, double damp)
double esccon(double tau, double hnukt)
STATIC void RTesc_lya_1side(double taume, double beta, realnum *esc, realnum *dest, long ipLine)
long int nSpecies
Definition taulines.cpp:21
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
double expn(int n, double x)
void VoigtU(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:364
void DumpLine(const TransitionProxy &t)
Wind wind
Definition wind.cpp:5