cloudy trunk
Loading...
Searching...
No Matches
rt_continuum_shield_fcn.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_continuum_shield_fcn computing continuum shielding due to single line */
4/*conpmp local continuum pumping rate radiative transfer for all lines */
5
6#include "cddefines.h"
7#include "rt.h"
8#include "transition.h"
9#include "thirdparty.h"
10
11/*conpmp local continuum pumping rate radiative transfer for all lines */
12STATIC double conpmp(const TransitionProxy& t);
13STATIC inline double FITTED( double t );
14
15namespace {
16 class my_Integrand
17 {
18 public:
19 double PumpDamp, PumpTau;
20
21 double operator() (double x)
22 {
23 realnum v, rx = realnum(x);
24 VoigtH(realnum(PumpDamp),&rx,&v,1);
25 double opfun_v = sexp(PumpTau*v)*v;
26 return opfun_v;
27 }
28};
29}
30
31
32/*rt_continuum_shield_fcn computing continuum shielding due to single line */
34{
35 double value;
36
37 DEBUG_ENTRY( "rt_continuum_shield_fcn()" );
38
39 value = -1.;
40
41 ASSERT( t.Emis().damp() > 0. );
42
43 if( rt.nLineContShield == LINE_CONT_SHIELD_PESC )
44 {
45 /* set continuum shielding pesc - shieding based on escape probability */
46 if( t.Emis().iRedisFun() == ipPRD )
47 {
48 value = esc_PRD_1side(t.Emis().TauCon(),t.Emis().damp());
49 }
50 else if( t.Emis().iRedisFun() == ipCRD )
51 {
52 value = esca0k2(t.Emis().TauCon());
53 }
54 else if( t.Emis().iRedisFun() == ipCRDW )
55 {
56 value = esc_CRDwing_1side(t.Emis().TauCon(),t.Emis().damp());
57 }
58 else if( t.Emis().iRedisFun() == ipLY_A )
59 {
60 value = esc_PRD_1side(t.Emis().TauCon(),t.Emis().damp());
61 }
62 else
64 }
65 else if( rt.nLineContShield == LINE_CONT_SHIELD_FEDERMAN )
66 {
67 /* set continuum shielding Federman - this is the default */
68 double core, wings;
69
70 /* these expressions implement the appendix of
71 * >>refer line shielding Federman, S.R., Glassgold, A.E., &
72 * >>refercon Kwan, J. 1979, ApJ, 227, 466 */
73 /* doppler core - equation A8 */
74 if( t.Emis().TauCon() < 2. )
75 {
76 core = sexp( t.Emis().TauCon() * 0.66666 );
77 }
78 else if( t.Emis().TauCon() < 10. )
79 {
80 core = 0.638 * pow(t.Emis().TauCon(),(realnum)-1.25f );
81 }
82 else if( t.Emis().TauCon() < 100. )
83 {
84 core = 0.505 * pow(t.Emis().TauCon(),(realnum)-1.15f );
85 }
86 else
87 {
88 core = 0.344 * pow(t.Emis().TauCon(),(realnum)-1.0667f );
89 }
90
91 /* do we add damping wings? */
92 wings = 0.;
93 if( t.Emis().TauCon() > 1.f && t.Emis().damp()>0. )
94 {
95 /* equation A6 */
96 double t1 = 3.02*pow(t.Emis().damp()*1e3,-0.064 );
97 double u1 = sqrt(t.Emis().TauCon()*t.Emis().damp() )/SDIV(t1);
98 wings = t.Emis().damp()/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );
99 /* add very large optical depth tail to converge this with respect
100 * to escape probabilities - if this function falls off more slowly
101 * than escape probability then upper level will become overpopulated.
102 * original paper was not intended for this regime */
103 if( t.Emis().TauCon()>1e7 )
104 wings *= pow( t.Emis().TauCon()/1e7,-1.1 );
105 }
106 value = core + wings;
107 /* some x-ray lines have vastly large damping constants, greater than 1.
108 * in these cases the previous wings value does not work - approximation
109 * is for small damping constant - do not let pump efficiency exceed unity
110 * in this case */
111 if( t.Emis().TauCon()>0. )
112 value = MIN2(1., value );
113 }
114 else if( rt.nLineContShield == LINE_CONT_SHIELD_FERLAND )
115 {
116 /* set continuum shielding ferland */
117 value = conpmp( t );
118 }
119 else if( rt.nLineContShield == 0 )
120 {
121 /* set continuum shielding none */
122 value = 1.;
123 }
124 else
125 {
127 }
128
129 /* the returned pump shield function must be greater than zero,
130 * and less than 1 if a maser did not occur */
131 ASSERT( value>=0 && (value<=1.||t.Emis().TauCon()<0.) );
132
133 return value;
134}
135
137
138static const double BREAK = 3.;
139/* fit to results for tau less than 10 */
140STATIC inline double FITTED( double t )
141{
142 return (0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
143}
144
145/*conpmp local continuum pumping rate radiative transfer for all lines */
147{
148 double a0,
149 conpmp_v,
150 tau,
151 yinc1,
152 yinc2;
153
154 DEBUG_ENTRY( "conpmp()" );
155
156 /* tau used will be optical depth in center of next zone
157 * >>chng 96 july 6, had been ipLnTauIn, did not work when sphere static set */
158 tau = t.Emis().TauCon();
159 /* compute pumping probability */
160 if( tau <= 10. )
161 {
162 /* for tau<10 a does not matter, and one fit for all */
163 conpmp_v = FITTED(tau);
164 }
165 else if( tau > 1e6 )
166 {
167 /* this far in winds line opacity well below electron scattering
168 * so ignore for this problem */
169 conpmp_v = 0.;
170 }
171 else
172 {
173 my_Integrand func;
174 func.PumpDamp = t.Emis().damp();
175 func.PumpTau = tau;
177
178 yinc1 = opfun.sum( 0., BREAK, func );
179 yinc2 = opfun.sum( BREAK, 100., func );
180
181 a0 = 0.886227*(1. + func.PumpDamp);
182 conpmp_v = (yinc1 + yinc2)/a0;
183 }
184
185 /* EscProb is escape probability, will not allow conpmp to be greater than it
186 * on second iteration with thick lines, pump prob=1 and esc=0.5
187 * conpmp = MIN( conpmp , t.t(ipLnEscP) )
188 * */
189 return conpmp_v;
190}
static double a0[83]
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define POW2
Definition cddefines.h:929
float realnum
Definition cddefines.h:103
const int ipCRDW
Definition cddefines.h:294
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipCRD
Definition cddefines.h:292
const int ipLY_A
Definition cddefines.h:296
const int ipPRD
Definition cddefines.h:290
realnum & damp() const
Definition emission.h:563
realnum & TauCon() const
Definition emission.h:453
int & iRedisFun() const
Definition emission.h:403
double sum(double min, double max, Integrand func)
Definition cddefines.h:1531
EmissionList::reference Emis() const
Definition transition.h:408
t_rt rt
Definition rt.cpp:5
double esc_PRD_1side(double tau, double a)
#define LINE_CONT_SHIELD_FERLAND
Definition rt.h:292
#define LINE_CONT_SHIELD_PESC
Definition rt.h:290
double esca0k2(double taume)
#define LINE_CONT_SHIELD_FEDERMAN
Definition rt.h:291
double esc_CRDwing_1side(double tau, double a)
STATIC double FITTED(double t)
double RT_continuum_shield_fcn(const TransitionProxy &t)
STATIC double conpmp(const TransitionProxy &t)
static const double BREAK
void VoigtH(realnum a, const realnum v[], realnum y[], int n)
Definition thirdparty.h:350