-
Notifications
You must be signed in to change notification settings - Fork 5
/
cgsconstants.f90
292 lines (229 loc) · 11.5 KB
/
cgsconstants.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
!>
!! \brief This module contains definitions of physical constants
!! and conversion factors, as well as routines to set recombination
!! and collisional ionization coefficients.
!!
!! Module for Capreole / C2-Ray (f90)
!!
!! \b Author: Garrelt Mellema, Martina M. Friedrich
!!
!! \b Date: 2015-06-30 (previously undated)
!!
!! \b Version: cgs units
!!
module cgsconstants
use precision, only: dp
use mathconstants, only: pi
implicit none
! A collection of physical constants and conversion factors
! Units: cgs
!> proton mass
real(kind=dp), parameter :: m_p=1.672661e-24_dp
!> speed of light
real(kind=dp), parameter :: c=2.997925e+10_dp
!> Planck constant
real(kind=dp), parameter :: hplanck=6.6260755e-27_dp
!> Stefan-Boltzmann constant
real(kind=dp), parameter :: sigma_SB=5.670e-5_dp
!> Boltzmann constant
real(kind=dp), parameter :: k_B=1.381e-16_dp
!> Gravitational constant
real(kind=dp), parameter :: G_grav=6.6732e-8_dp
!> ev2k - conversion factor between evs and kelvins
real(kind=dp),parameter :: ev2k=1.0/8.617e-05
!> ev2erg - conversion factor between evs and ergs
real(kind=dp),parameter :: ev2erg=1.602e-12
!> ev2j - conversion factor between ergs and Joules
real(kind=dp),parameter :: erg2j=1e-7
! The following are scaled to frequency scaling
!> Frequency scaling factor,
!! this scaling parameter is independent of any main program scaling
!! (see scaling.f90), and may only be used in the radiation physics
!! subroutines (currently switched off)
real(kind=dp),parameter :: sclfre=1.0e15
!> conversion between evs and frequency
real(kind=dp),parameter :: ev2fr=0.241838e15!/sclfre
! h/k, Planck/Boltzm
! Check this number, seems scaled
!real(kind=dp),parameter :: hoverk=47979.72484
!> Planck constant scaled
real(kind=dp),parameter :: hscl=hplanck!*sclfre
!> tpic2 - 2*pi/c^2 times scaling factors needed for the integral cores
real(kind=dp),parameter :: two_pi_over_c_square=2.0*pi/(c*c)
!> Hydrogen recombination parameter (power law index)
real(kind=dp),parameter :: albpow=-0.7_dp !in the old code -0.79
!> Hydrogen recombination parameter (value at 10^4 K)
real(kind=dp),parameter :: bh00=2.59e-13_dp ! OTS value, alpha_B
!> Helium0 recombination parameter (power law index)
real(kind=dp), parameter :: alcpow=-0.672_dp
!> Helium0 recombination parameter (value at 10^4 K)
real(kind=dp), parameter :: bhe00=4.26e-13_dp !alpha_b+alpha_1
!> Helium1 recombination parameter (value at 10^4 K)
real(kind=dp), parameter :: bhe10=1.53e-12_dp !different in the book!
!here it was 5.91e-12 I replace with book value of 2*1.1e-12
!> Hydrogen ionization energy (in eV)
real(kind=dp), parameter :: eth0=13.598
!> Hydrogen ionization energy (in erg)
real(kind=dp),parameter :: hionen=eth0*ev2erg
!> Hydrogen ionization energy expressed in K
real(kind=dp),parameter :: temph0=eth0*ev2k
!> Hydrogen collisional ionization parameter 1
real(kind=dp),parameter :: xih0=1.0
!> Hydrogen collisional ionization parameter 2
real(kind=dp),parameter :: fh0=0.83
!> Hydrogen collisional ionization parameter
real(kind=dp),parameter :: colh0=1.3e-8*fh0*xih0/(eth0*eth0)
! critical electron density at which colisions become important for y-fraction ionization (Osterbrock)
real(kind=dp),parameter :: n_el_crit=4.0e3
!> Helium ionization energy (in eV)
real(kind=dp), dimension(0:1), parameter :: ethe=(/24.587,54.416/)
!> Helium ionization energy (in erg)
real(kind=dp), dimension(0:1), parameter :: heionen=(/ethe(0)*ev2erg,ethe(1)*ev2erg/)
!> Helium ionization energy expressed in K
real(kind=dp), dimension(0:1), parameter :: temphe=(/ethe(0)*ev2k,ethe(1)*ev2k/)
!> Helium collisional ionization parameter 1
real(kind=dp),dimension(0:1),parameter :: xihe=(/2.0,1.0/)
!> Helium collisional ionization parameter 2
real(kind=dp),dimension(0:1),parameter :: fhe=(/0.63,1.30/)
!> Helium collisional ionization parameter
real(kind=dp),dimension(0:1),parameter :: colhe= &
(/1.3e-8*fhe(0)*xihe(0)/(ethe(0)*ethe(0)), &
1.3e-8*fhe(1)*xihe(1)/(ethe(1)*ethe(1))/)
!> Hydrogen 0 A recombination parameter
real(kind=dp) :: arech0
!> Hydrogen 0 B recombination parameter
real(kind=dp) :: brech0
!> Helium 0 A recombination parameter
real(kind=dp) :: areche0
!> Helium 0 B recombination parameter
real(kind=dp) :: breche0
!> Helium 0 1 recombination parameter
real(kind=dp) :: oreche0
!> Helium + A recombination parameter
real(kind=dp) :: areche1
!> Helium + B recombination parameter
real(kind=dp) :: breche1
!> Helium + 2 recombination parameter
real(kind=dp) :: treche1
!> H0 collisional ionization parameter at T=temp0
real(kind=dp) :: colli_HI
!> He0 collisional ionization parameter at T=temp0
real(kind=dp) :: colli_HeI
!> He1 collisional ionization parameter at T=temp0
real(kind=dp) :: colli_HeII
!> Fraction fo He++ -> He+ recombination photons that goes into 2photon decay
real(kind=dp) :: v
contains
!> This subroutine contains all temperature dependent recombination and
!! collisional ionization rates. It should be called every time after the
!! temperature is recalculated.
subroutine ini_rec_colion_factors (temperature)
real(kind=dp),intent(in):: temperature
call ini_hydrogen_recombination (temperature)
call ini_helium0_recombination (temperature)
call ini_helium1_recombination (temperature)
call ini_hydrogen_helium_collisional_ionization (temperature)
end subroutine ini_rec_colion_factors
! ---------------------------------------------------------------------------
subroutine ini_hydrogen_recombination (temperature)
! Recombination from H^+ to H^0 coefficients.
! MMF: Implement the new fit as mentioned in email from 2008-06-08
!
! fit1=bh00*1.4*(T/1e4)^albpow*0.9
! fit2=bh00*5*(T/1e4)^albpow*1.95
! brech0=1/[(1/fit1)+(1/fit2)]
!
real(kind=dp),intent(in):: temperature
!real(kind=dp) :: fit1,fit2
real(kind=dp) :: lambda
! Hydrogen A & B recombination rates, fits from Hui & Gnedin (1997)
lambda = 2.0_dp*(temph0/temperature)
arech0=1.269e-13*lambda**1.503_dp/(1.0_dp+(lambda/0.522)**0.470)**1.923
brech0=2.753e-14*lambda**1.500_dp/(1.0_dp+(lambda/2.740)**0.407)**2.242
end subroutine ini_hydrogen_recombination
!----------------------------------------------------------------------------
subroutine ini_helium0_recombination (temperature)
! Recombination from He^+ to He^0 coefficients.
! MMF: See email from 2008-08-11
! For now, use old fit.
real(kind=dp),intent(in):: temperature
real(kind=dp) :: lambda
real(kind=dp) :: dielectronic
! He0 recombination rates including the dielectronic recombination
if (temperature < 9.e3_dp) then ! was 1.5e4
lambda=2.0_dp*(temph0/temperature)
areche0=1.269e-13_dp*lambda**1.503_dp/(1.0_dp+(lambda/0.522)**0.470)**1.923
breche0=2.753e-14_dp*lambda**1.500_dp/(1.0_dp+(lambda/2.740)**0.407)**2.242
else
lambda= 2.0_dp*(temphe(0)/temperature)
dielectronic = 1.9e-3_dp*temperature**(-1.5_dp)*exp(-4.7e5_dp/temperature)*(1.0_dp+0.3_dp*exp(-9.4e4_dp/temperature))
areche0= 3.000e-14_dp*lambda**0.654_dp + dielectronic
breche0= 1.260e-14_dp*lambda**0.750_dp + dielectronic
endif
oreche0=areche0-breche0
! above 1.5e4 I follow the fit from Hui & Gnedin (1997) for He0 for radiation recombination,
! this fits the data from Hummer&Storey (1998) good.
! below that, the data from Hummer&Storey is fitted by the fit for H from Hui%Gnedin1997
! comparing the He0 and H data from Hummer&Storey (1998) and Hummer (1994) shows, that
! the two never diviate more than 4% above T=1.5e4.
! without dielectronic recombinations
!areche0 = bhe00*(temperature/1.0e4_dp)**alcpow
!oreche0 =1.54e-13_dp*(temperature/1.0e4_dp)**(-0.5_dp) ! values from Osterbrock interpolated
!breche0 = areche0-oreche0 !B recombination
! NOT USED ANYMORE
end subroutine ini_helium0_recombination
!----------------------------------------------------------------------------
subroutine ini_helium1_recombination (temperature)
! Recombination from He^2+ to He^+ coefficients.
! MMF: Implement the new fit as mentioned in email from 2008-06-08
real(kind=dp),intent(in):: temperature
real(kind=dp) :: lambda
! fits to He1 recombination rates to the data from Hummer (1994)
! using the fits from Hui & Gnedin (1997), they are excellent fits
! to the table data from Hummer (1994)
! The same is true for the A rates
lambda=2.0_dp*(temphe(1)/temperature)
breche1=5.5060e-14_dp*lambda**1.5_dp/(1.0_dp+(lambda/2.740_dp)**0.407_dp)**2.242_dp
areche1=2.538e-13*lambda**1.503_dp/(1.0_dp+(lambda/0.522_dp)**0.470_dp)**1.923_dp
! For the B-C values, that is recombination to n=2, there is no data in
! Hummer (1994), neither a fit in Hui&Gnedin (1997)
treche1 = 3.4e-13_dp*(temperature/1.0e4_dp)**(-0.6_dp) ! extrapolate Osterbrok B value minus C value, p 38
! fraction of He++ recombination photons that go into two-photon decay.
v=0.285_dp*(temperature/1.0e4_dp)**0.119_dp ! fit to Hummer and Seaton, 1964 table five
end subroutine ini_helium1_recombination
!----------------------------------------------------------------------------
subroutine ini_hydrogen_helium_collisional_ionization (temperature)
! Collisional ionization coefficients for H^0 to H^+ (temperature dependent)
real(kind=dp),intent(in):: temperature
real(kind=dp) :: sqrtt0
! Fits from Cox (1970)
sqrtt0 =sqrt(temperature)
colli_HI =colh0*sqrtt0*exp(-temph0/temperature)
colli_HeI=colhe(0)*sqrtt0*exp(-temphe(0)/temperature)
colli_HeII=colhe(1)*sqrtt0*exp(-temphe(1)/temperature)
! I might consider using the Hui&Gnedin fits instead. However, the difference is not so big
!sqrtt0 = temperature**(-1.5_dp)*exp(-temph0/temperature)
!acolh0 =21.11_dp*sqrtt0*exp(-temph0/temperature)* &
!(2.0_dp*temph0/temperature)**-1.089_dp/(1.0_dp+((2.0_dp*temph0/temperature)/0.354_dp)**0.874)**1.101_dp
!acolhe0=32.38_dp*sqrtt0*exp(-temphe(0)/temperature)* &
! (2.0_dp*temphe(0)/temperature)**-1.146_dp/(1.0_dp+((2.0_dp*temphe(0)/temperature)/0.416_dp)**0.987)**1.056;
!acolhe1=19.95_dp*sqrtt0*exp(-temphe(1)/temperature)* &
!(2.0_dp*temphe(1)/temperature)**-1.089_dp/(1.0_dp+((2.0_dp*temphe(1)/temperature)/0.553_dp)**0.735)**1.275;
end subroutine ini_hydrogen_helium_collisional_ionization
!----------------------------------------------------------------------------
subroutine ini_recombinations_collisional_ionizations_constant
! SET THEM HERE EXPLIXITLY FOR T=10 000 K TO THE ONES FROM GABRIEL
brech0 = 2.59182e-13_dp ! H0 b-recombination coefficient
breche0 = 2.61613e-13_dp ! He0 b-recombination coefficient
breche1 = 1.54528e-12_dp ! He1 b-recombination coefficient
areche0 = 4.22471e-13_dp ! He0 a-recombination coefficient
areche1 = 2.22561e-12_dp ! He1 a-recombination coefficient
arech0 = 4.29695e-13_dp ! H0 a-recombination coefficient
oreche0 =areche0-breche0 ! he0 1 recombination coefficient
treche1 = 3.46e-13_dp !Osterbrok B value minus C value, p 38
v=0.285_dp
colli_HI = 8.96396e-16_dp ! H collisional ionization coefficient
colli_HeI = 7.46415e-22_dp ! He0 collisional ionization coefficient
colli_HeII = 2.28059e-37_dp ! He1 collisional ionization coefficient
end subroutine ini_recombinations_collisional_ionizations_constant
end module cgsconstants