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
highen.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
/*highen do high energy radiation field - gas interaction, Compton scattering, etc */
4
#include "
cddefines.h
"
5
#include "
physconst.h
"
6
#include "
trace.h
"
7
#include "
heavy.h
"
8
#include "
radius.h
"
9
#include "
magnetic.h
"
10
#include "
hextra.h
"
11
#include "
thermal.h
"
12
#include "
dense.h
"
13
#include "
doppvel.h
"
14
#include "
ionbal.h
"
15
#include "
rfield.h
"
16
#include "
opacity.h
"
17
#include "
gammas.h
"
18
#include "
highen.h
"
19
20
void
highen
(
void
)
21
{
22
long
int
i,
23
ion,
24
nelem;
25
26
double
CosRayDen,
27
crsphi,
28
heatin,
29
sqthot;
30
31
double
hin;
32
33
DEBUG_ENTRY
(
"highen()"
);
34
35
36
/**********************************************************************
37
*
38
* COMPTON RECOIL IONIZATION
39
*
40
* bound electron scattering of >2.3 kev photons if neutral
41
* comoff usually 1, 0 if "NO COMPTON EFFECT" command given
42
* lgCompRecoil usually t, f if "NO RECOIL IONIZATION" command given
43
*
44
**********************************************************************/
45
46
/* comoff turned off with no compton command,
47
* lgCompRecoil - no recoil ionization */
48
if
(
rfield
.
comoff
> .0 &&
ionbal
.
lgCompRecoil
)
49
{
50
for
( nelem=0; nelem<
LIMELM
; ++nelem )
51
{
52
for
( ion=0; ion<nelem+1; ++ion )
53
{
54
ionbal
.
CompRecoilIonRate
[nelem][ion] = 0.;
55
ionbal
.
CompRecoilHeatRate
[nelem][ion] = 0.;
56
if
(
dense
.
xIonDense
[nelem][ion] > 0. )
57
{
58
/* recoil ionization starts at 194 Ryd = 2.6 keV */
59
/* this is the ionization potential of the valence shell */
60
/* >>chng 02 may 27, lower limit is now 1 beyond actual threshold
61
* since recoil energy at threshold was very small, sometimes negative */
62
for
( i=
ionbal
.
ipCompRecoil
[nelem][ion]; i <
rfield
.
nflux
; ++i)
63
{
64
double
recoil_energy;
65
crsphi =
opac
.
OpacStack
[i-1+
opac
.
iopcom
] *
rfield
.
SummedCon
[i] *
66
/* this is number of electrons in valence shell of this ion */
67
ionbal
.
nCompRecoilElec
[nelem-ion];
68
69
/* direct hydrogen ionization due to compton scattering,
70
* does not yet include secondaries,
71
* last term accounts for number of valence electrons that contribute */
72
ionbal
.
CompRecoilIonRate
[nelem][ion] += crsphi;
73
74
/* recoil energy in Rydbergs
75
* heating modified for suprathermal secondaries below; ANU2=ANU**2 */
76
/* >>chng 02 mar 27, use general formula for recoil energy */
77
/*energy = 2.66e-5*rfield.anu2[i] - 1.;*/
78
recoil_energy =
rfield
.
anu2
[i] /
79
(
EN1RYD
*
ELECTRON_MASS
*
SPEEDLIGHT
*
SPEEDLIGHT
) *
EN1RYD
*
EN1RYD
-
80
Heavy
.
Valence_IP_Ryd
[nelem][ion];
81
82
/* heating is in rydbergs because SecIon2PrimaryErg, SecExcitLya2PrimaryErg, HeatEfficPrimary in ryd */
83
ionbal
.
CompRecoilHeatRate
[nelem][ion] += crsphi*recoil_energy;
84
}
85
/* net heating rate, per atom, convert ryd/sec/cm3 to ergs/sec/atom */
86
ionbal
.
CompRecoilHeatRate
[nelem][ion] *=
EN1RYD
;
87
88
ASSERT
(
ionbal
.
CompRecoilHeatRate
[nelem][ion] >= 0.);
89
ASSERT
(
ionbal
.
CompRecoilIonRate
[nelem][ion] >= 0.);
90
}
91
}
92
}
93
}
94
else
95
{
96
for
( nelem=0; nelem<
LIMELM
; ++nelem )
97
{
98
for
( ion=0; ion<nelem+1; ++ion )
99
{
100
ionbal
.
CompRecoilIonRate
[nelem][ion] = 0.;
101
ionbal
.
CompRecoilHeatRate
[nelem][ion] = 0.;
102
}
103
}
104
}
105
106
/**********************************************************************
107
*
108
* COSMIC RAYS
109
*
110
* heating and ionization by cosmic rays, other relativistic particles
111
* CRYDEN=density (1/CM**3), neutral rate assumes 15ev total
112
* energy loss, 13.6 into ionization, 1.4 into heating
113
* units erg/sec/cm**3
114
* iff not specified, CRTEMP is 2.6E9K
115
*
116
**********************************************************************/
117
118
if
(
hextra
.
cryden
> 0. )
119
{
120
ASSERT
(
hextra
.
crtemp
> 0. );
121
/* this is current cosmic ray density, as determined from original density times
122
* possible dependence on radius */
123
if
(
hextra
.
lg_CR_B_equipartition
)
124
{
125
/* >>chng 06 jun 02, add this option
126
* this is case where cr are in equipartition with magnetic field -
127
* set with COSMIC RAY EQUIPARTITION command */
128
CosRayDen =
hextra
.
background_density
*
129
/* ratio of energy density in current B to typical galactic
130
* galactic background energy density of 1.8 eV cm-3 is from
131
*>>refer cr background Webber, W.R. 1998, ApJ, 506, 329 */
132
magnetic
.
energydensity
/
133
(
CR_EDEN_GAL_BACK_EV_CMM3
/*1.8eV cm-3*/
*
EN1EV
/*erg eV-1*/
);
134
}
135
else
136
{
137
/* this is usual case, CR density may depend on radius, usually does not */
138
CosRayDen =
hextra
.
cryden
*pow(
radius
.
Radius
/
radius
.
rinner
,(
double
)
hextra
.
crpowr
);
139
}
140
141
/* cosmic ray energy density rescaled by ratio to background ion rate and B field */
142
hextra
.
cr_energydensity
= CosRayDen/
hextra
.
background_density
*
143
(
CR_EDEN_GAL_BACK_EV_CMM3
/*1.8eV cm-3*/
*
EN1EV
/*erg eV-1*/
);
144
145
/* related to current temperature, when thermal */
146
sqthot = sqrt(
hextra
.
crtemp
);
147
148
/* rate hot electrons heat cold ones, Balbus and McKee 1982
149
* units erg sec^-1 cm^-3,
150
* in sumheat we will multipy this rate by sum of neturals, but for this
151
* term we actually want eden, so mult by eden over sum of neut */
152
ionbal
.
CosRayHeatThermalElectrons
= 5.5e-14/sqthot*CosRayDen;
153
154
/* ionization rate; Balbus and McKee */
155
ionbal
.
CosRayIonRate
= 1.22e-4/sqthot*
156
log(2.72*pow(
hextra
.
crtemp
/1e8,0.097))*CosRayDen;
157
158
/* option for thermal CRs, first is the usual (and default) relativistic case */
159
if
(
hextra
.
crtemp
> 2e9 )
160
{
161
/* usual circumstance; relativistic cosmic rays,
162
* cosmic ray ionization rate s-1 cm-3; ext rel limit */
163
ionbal
.
CosRayIonRate
*= 3.;
164
165
}
166
else
167
{
168
/* option for thermal cosmic rays */
169
ionbal
.
CosRayIonRate
*= 10.;
170
}
171
/* >>chng 04 jan 27, from 0.093 to 2.574 as per following */
172
/* cr heating from Table 1 of
173
*>>refer cr heating Wolfire et al.1995, ApJ, 443, 152
174
* For every ionization due to cosmic rays, ~35eV of heat is added
175
* to the system. This manifests itself in the ionbal.CosRayHeatNeutralParticles term
176
* by the 2.574*EN1RYD term, which is just the energy in ergs in 35 eV.
177
* Change made by Nick Abel and Gargi Shaw, 04 Jan 27. In heatsum
178
* we multiply by the number of secondaries that occur */
179
ionbal
.
CosRayHeatNeutralParticles
=
ionbal
.
CosRayIonRate
*2.574*
EN1RYD
;
180
181
if
(
trace
.
lgTrace
)
182
{
183
fprintf(
ioQQQ
,
" highen: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n"
,
184
CosRayDen,
ionbal
.
CosRayIonRate
,
ionbal
.
CosRayHeatNeutralParticles
,
hextra
.
crtemp
);
185
}
186
}
187
else
188
{
189
ionbal
.
CosRayIonRate
= 0.;
190
ionbal
.
CosRayHeatNeutralParticles
= 0.;
191
}
192
/* >>chng 06 may 23, Penning ionization
193
ionbal.CosRayIonRate += 1e-9 *
194
StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*dense.xIonDense[ipHELIUM][1]; */
195
196
/*fprintf(ioQQQ,"DEBUG cr %.2f %.3e %.3e %.3e\n",
197
fnzone,
198
hextra.cryden ,
199
ionbal.CosRayIonRate ,
200
ionbal.CosRayHeatNeutralParticles );*/
201
202
/**********************************************************************
203
*
204
* add on extra heating due to turbulence, goes into [1] of [x][0][11][0]
205
*
206
**********************************************************************/
207
208
/* TurbHeat added with hextra command, DispScale added with turbulence dissipative */
209
if
( (
hextra
.
TurbHeat
+
DoppVel
.
DispScale
) != 0. )
210
{
211
/* turbulent heating only goes into the low-energy heat of this element */
212
/* >>>>chng 00 apr 28, functional form of radius dependence had bee turrad/depth
213
* and so went to infinity at the illuminated face. Changed to current form as
214
* per Ivan Hubeny comment */
215
if
(
hextra
.
lgHextraDepth
)
216
{
217
/* if turrad is >0 then vary heating with depth */
218
ionbal
.
ExtraHeatRate
=
219
hextra
.
TurbHeat
*
sexp
(
radius
.
depth
/
hextra
.
turrad
);
220
221
/* >>chng 00 aug 16, added option for heating from back side */
222
if
(
hextra
.
turback
!= 0. )
223
{
224
ionbal
.
ExtraHeatRate
+=
225
hextra
.
TurbHeat
*
sexp
((
hextra
.
turback
-
radius
.
depth
) /
hextra
.
turrad
);
226
}
227
}
228
else
if
(
hextra
.
lgHextraDensity
)
229
{
230
/* depends on density */
231
ionbal
.
ExtraHeatRate
=
232
hextra
.
TurbHeat
*
dense
.
gas_phase
[
ipHYDROGEN
]/
hextra
.
HextraScaleDensity
;
233
}
234
/* this is turbulence dissipate command */
235
else
if
(
DoppVel
.
DispScale
> 0. )
236
{
237
double
turb =
DoppVel
.
TurbVel
*
sexp
(
radius
.
depth
/
DoppVel
.
DispScale
);
238
/* if turrad is >0 then vary heating with depth */
239
/* >>chng 01 may 10, add extra factor of length over 1e13 cm */
240
ionbal
.
ExtraHeatRate
= 3.45e-28 / 2.82824 * turb * turb * turb
241
/ (
DoppVel
.
DispScale
/1e13);
242
}
243
else
244
{
245
/* constant extra heating */
246
ionbal
.
ExtraHeatRate
=
hextra
.
TurbHeat
;
247
}
248
}
249
250
else
251
{
252
ionbal
.
ExtraHeatRate
= 0.;
253
}
254
255
/**********************************************************************
256
*
257
* option to add on fast neutron heating, goes into [0] & [2] of [x][0][11][0]
258
*
259
**********************************************************************/
260
if
(
hextra
.
lgNeutrnHeatOn
)
261
{
262
/* hextra.totneu is energy flux erg cm-2 s-1
263
* CrsSecNeutron is 4E-26 cm^-2, cross sec for stopping rel neutrons
264
* this is heating erg s-1 due to fast neutrons, assumed to secondary ionize */
265
/* neutrons assumed to only secondary ionize */
266
ionbal
.
xNeutronHeatRate
=
hextra
.
totneu
*
hextra
.
CrsSecNeutron
;
267
}
268
else
269
{
270
ionbal
.
xNeutronHeatRate
= 0.;
271
}
272
273
274
/**********************************************************************
275
*
276
* pair production in elec field of nucleus
277
*
278
**********************************************************************/
279
ionbal
.
PairProducPhotoRate
[0] =
GammaK
(
opac
.
ippr
,
rfield
.
nflux
,
opac
.
ioppr
,1.);
280
ionbal
.
PairProducPhotoRate
[1] =
thermal
.
HeatLowEnr
;
281
ionbal
.
PairProducPhotoRate
[2] =
thermal
.
HeatHiEnr
;
282
283
/**********************************************************************
284
*
285
* Compton energy exchange
286
*
287
**********************************************************************/
288
rfield
.
cmcool
= 0.;
289
rfield
.
cmheat
= 0.;
290
heatin = 0.;
291
/* comoff usually 1, turns off Compton */
292
if
(
rfield
.
comoff
>0.01 )
293
{
294
for
( i=0; i <
rfield
.
nflux
; i++ )
295
{
296
297
/* Compton cooling
298
* CSIGC is Tarter expression times ANU(I)*3.858E-25
299
* 6.338E-6 is k in inf mass Rydbergs, still needs factor of TE */
300
rfield
.
comup
[i] = (double)(
rfield
.
flux
[i]+
rfield
.
ConInterOut
[i]+
301
rfield
.
outlin
[i]+
rfield
.
outlin_noplot
[i])*
rfield
.
csigc
[i]*(
dense
.
eden
*4.e0*
302
6.338e-6*1e-15);
303
304
rfield
.
cmcool
+=
rfield
.
comup
[i];
305
306
/* Compton heating
307
* CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25
308
* CMHEAT is just spontaneous, HEATIN is just induced */
309
rfield
.
comdn
[i] = (double)(
rfield
.
flux
[i]+
rfield
.
ConInterOut
[i]+
310
rfield
.
outlin
[i]+
rfield
.
outlin_noplot
[i])*
rfield
.
csigh
[i]*
dense
.
eden
*1e-15;
311
312
/* induced Compton heating */
313
hin = (double)(
rfield
.
flux
[i]+
rfield
.
ConInterOut
[i]+
rfield
.
outlin
[i]+
rfield
.
outlin_noplot
[i])*
314
rfield
.
csigh
[i]*
rfield
.
OccNumbIncidCont
[i]*
dense
.
eden
*1e-15;
315
rfield
.
comdn
[i] += hin;
316
heatin += hin;
317
318
/* following is total compton heating */
319
rfield
.
cmheat
+=
rfield
.
comdn
[i];
320
}
321
322
/* remember how important induced compton heating is */
323
if
(
rfield
.
cmheat
> 0. )
324
{
325
rfield
.
cinrat
=
MAX2
(
rfield
.
cinrat
,heatin/
rfield
.
cmheat
);
326
}
327
328
if
(
trace
.
lgTrace
&&
trace
.
lgComBug
)
329
{
330
fprintf(
ioQQQ
,
" HIGHEN: COOL num=%8.2e HEAT num=%8.2e\n"
,
331
rfield
.
cmcool
,
rfield
.
cmheat
);
332
}
333
}
334
335
/* save compton heating rate into main heating save array,
336
* no secondary ionizations from compton heating cooling */
337
thermal
.
heating
[0][19] =
rfield
.
cmheat
;
338
339
if
(
trace
.
lgTrace
&&
trace
.
lgComBug
)
340
{
341
fprintf(
ioQQQ
,
342
" HIGHEN finds heating fracs= frac(compt)=%10.2e "
343
" f(pair)%10.2e totHeat=%10.2e\n"
,
344
thermal
.
heating
[0][19]/
thermal
.
htot
,
345
thermal
.
heating
[0][21]/
thermal
.
htot
,
346
thermal
.
htot
);
347
}
348
return
;
349
}
Generated for cloudy by
1.8.1.2