cloudy
trunk
Main Page
Related Pages
Data Structures
Files
File List
Globals
All
Data Structures
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
source
rt_continuum_shield_fcn.cpp
Go to the documentation of this file.
1
/* This file is part of Cloudy and is copyright (C)1978-2008 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
/*conpmp local continuum pumping rate radiative transfer for all lines */
6
#include "
cddefines.h
"
7
#include "
rt.h
"
8
/*conpmp local continuum pumping rate radiative transfer for all lines */
9
STATIC
double
conpmp
(
transition
* t);
10
STATIC
inline
double
FITTED
(
double
t );
11
static
double
PumpDamp
,
PumpTau
;
12
13
/*rt_continuum_shield_fcn computing continuum shielding due to single line */
14
double
RT_continuum_shield_fcn
(
transition
*t )
15
{
16
double
value;
17
18
DEBUG_ENTRY
(
"rt_continuum_shield_fcn()"
);
19
20
value = -1.;
21
22
if
(
rt
.
nLineContShield
==
LINE_CONT_SHIELD_PESC
)
23
{
24
/* set continuum shielding pesc - shieding based on escape probability */
25
if
( t->
Emis
->
iRedisFun
==
ipPRD
)
26
{
27
value =
esc_PRD_1side
(t->
Emis
->
TauCon
,t->
Emis
->
damp
);
28
}
29
else
if
( t->
Emis
->
iRedisFun
==
ipCRD
)
30
{
31
value =
esca0k2
(t->
Emis
->
TauCon
);
32
}
33
else
if
( t->
Emis
->
iRedisFun
==
ipCRDW
)
34
{
35
value =
esc_CRDwing_1side
(t->
Emis
->
TauCon
,t->
Emis
->
damp
);
36
}
37
else
if
( t->
Emis
->
iRedisFun
==
ipLY_A
)
38
{
39
value =
esc_PRD_1side
(t->
Emis
->
TauCon
,t->
Emis
->
damp
);
40
}
41
else
42
TotalInsanity
();
43
}
44
else
if
(
rt
.
nLineContShield
==
LINE_CONT_SHIELD_FEDERMAN
)
45
{
46
/* set continuum shielding Federman - this is the default */
47
double
core, wings;
48
49
/* these expressions implement the appendix of
50
* >>refer line shielding Federman, S.R., Glassgold, A.E., &
51
* >>refercon Kwan, J. 1979, ApJ, 227, 466 */
52
/* doppler core - equation A8 */
53
if
( t->
Emis
->
TauCon
< 2. )
54
{
55
core =
sexp
( t->
Emis
->
TauCon
* 0.66666 );
56
}
57
else
if
( t->
Emis
->
TauCon
< 10. )
58
{
59
core = 0.638 * pow(t->
Emis
->
TauCon
,(
realnum
)-1.25f );
60
}
61
else
if
( t->
Emis
->
TauCon
< 100. )
62
{
63
core = 0.505 * pow(t->
Emis
->
TauCon
,(
realnum
)-1.15f );
64
}
65
else
66
{
67
core = 0.344 * pow(t->
Emis
->
TauCon
,(
realnum
)-1.0667f );
68
}
69
70
/* do we add damping wings? */
71
wings = 0.;
72
if
( t->
Emis
->
TauCon
> 1.f && t->
Emis
->
damp
>0. )
73
{
74
/* equation A6 */
75
double
t1 = 3.02*pow(t->
Emis
->
damp
*1e3,-0.064 );
76
double
u1 = sqrt(t->
Emis
->
TauCon
*t->
Emis
->
damp
)/
SDIV
(t1);
77
wings = t->
Emis
->
damp
/
SDIV
(t1)/sqrt( 0.78540 +
POW2
(u1) );
78
/* add very large optical depth tail to converge this with respect
79
* to escape probabilities - if this function falls off more slowly
80
* than escape probability then upper level will become overpopulated.
81
* original paper was not intended for this regime */
82
if
( t->
Emis
->
TauCon
>1e7 )
83
wings *= pow( t->
Emis
->
TauCon
/1e7,-1.1 );
84
}
85
value = core + wings;
86
/* some x-ray lines have vastly large damping constants, greater than 1.
87
* in these cases the previous wings value does not work - approximation
88
* is for small damping constant - do not let pump efficiency exceed unity
89
* in this case */
90
if
( t->
Emis
->
TauCon
>0. )
91
value =
MIN2
(1., value );
92
}
93
else
if
(
rt
.
nLineContShield
==
LINE_CONT_SHIELD_FERLAND
)
94
{
95
/* set continuum shielding ferland */
96
value =
conpmp
( t );
97
}
98
else
if
(
rt
.
nLineContShield
== 0 )
99
{
100
/* set continuum shielding none */
101
value = 1.;
102
}
103
else
104
{
105
TotalInsanity
();
106
}
107
108
/* the returned pump shield function must be greater than zero,
109
* and less than 1 if a maser did not occur */
110
ASSERT
( value>=0 && (value<=1.||t->Emis->TauCon<0.) );
111
112
return
value;
113
}
114
115
/*opfun routine used to get continuum pumping of lines
116
* used in conpmp in call to qg32 */
117
STATIC
double
opfun
(
double
x)
118
{
119
double
opfun_v,
120
v;
121
122
DEBUG_ENTRY
(
"opfun()"
);
123
124
v =
vfun
(
PumpDamp
,x);
125
opfun_v =
sexp
(
PumpTau
*v)*v;
126
return
( opfun_v );
127
}
128
129
static
const
double
BREAK
= 3.;
130
/* fit to results for tau less than 10 */
131
// #define FITTED(t) ((0.98925439 + 0.084594094*(t))/(1. + (t)*(0.64794212 + (t)*0.44743976)))
132
STATIC
inline
double
FITTED
(
double
t )
133
{
134
return
(0.98925439 + 0.084594094*t)/(1. + t*(0.64794212 + t*0.44743976));
135
}
136
137
/*conpmp local continuum pumping rate radiative transfer for all lines */
138
STATIC
double
conpmp
(
transition
* t)
139
{
140
double
a0
,
141
conpmp_v,
142
tau,
143
yinc1,
144
yinc2;
145
146
DEBUG_ENTRY
(
"conpmp()"
);
147
148
/* tau used will be optical depth in center of next zone
149
* >>chng 96 july 6, had been ipLnTauIn, did not work when sphere static set */
150
tau = t->
Emis
->
TauCon
;
151
/* compute pumping probability */
152
if
( tau <= 10. )
153
{
154
/* for tau<10 a does not matter, and one fit for all */
155
conpmp_v =
FITTED
(tau);
156
}
157
else
if
( tau > 1e6 )
158
{
159
/* this far in winds line opacity well below electron scattering
160
* so ignore for this problem */
161
conpmp_v = 0.;
162
}
163
else
164
{
165
/* following two are passed on to later subs */
166
PumpDamp
= t->
Emis
->
damp
;
167
PumpTau
= tau;
168
a0 = 0.886227*(1. +
PumpDamp
);
169
yinc1 =
qg32
(0.,
BREAK
,
opfun
);
170
yinc2 =
qg32
(
BREAK
,100.,
opfun
);
171
conpmp_v = (yinc1 + yinc2)/a0;
172
}
173
174
/* EscProb is escape probability, will not allow conpmp to be greater than it
175
* on second iteration with thick lines, pump prob=1 and esc=0.5
176
* conpmp = MIN( conpmp , t->t(ipLnEscP) )
177
* */
178
return
( conpmp_v );
179
}
Generated for cloudy by
1.8.1.2