Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module pycno
21 :
22 : use rates_def
23 : use utils_lib
24 : use const_def, only: dp
25 : use math_lib
26 :
27 : implicit none
28 :
29 : contains
30 :
31 2 : subroutine FL_epsnuc_3alf(T, Rho, Y, UE, r, drdT, drdRho)
32 : real(dp), intent(in) :: T ! temperature
33 : real(dp), intent(in) :: Rho ! density
34 : real(dp), intent(in) :: Y ! helium mass fraction
35 : real(dp), intent(in) :: UE ! electron molecular weight
36 : real(dp), intent(out) :: r ! rate in ergs/g/sec
37 : real(dp), intent(out) :: drdT ! partial wrt temperature
38 : real(dp), intent(out) :: drdRho ! partial wrt density
39 :
40 2 : real(dp) :: T6, R6, R6T, R6T13, R6T16, T62, T612, T623, T653, T632, T613, U, AF
41 2 : real(dp) :: G1, dG1dRho, dG1dT, G2, dG2dRho, dG2dT
42 2 : real(dp) :: B1, dB1dRho, dB1dT, B2, dB2dRho, dB2dT
43 2 : real(dp) :: dUdT, dUdRho, U32, U52, dAFdT, dAFdRho
44 2 : real(dp) :: E1, dE1dT, dE1dRho, E2, dE2dT, dE2dRho
45 2 : real(dp) :: F1, dF1dT, dF1dRho, F2, dF2dT, dF2dRho
46 2 : real(dp) :: dR6dRho, dR6TdRho, dR6T13dRho, dR6T16dRho
47 2 : real(dp) :: dT6dT, dT612dT, dT62dT, dT613dT, dT623dT, dT632dT, dT653dT
48 :
49 : ! DEBUG
50 : real(dp), parameter :: AF_0 = 1.9005324047511074D+00
51 : real(dp), parameter :: B1_denom_0 = 2.9602238143383192D-01
52 : real(dp), parameter :: B1_0 = 1.2227955158250397D-08
53 : real(dp), parameter :: B2_denom_0 = 1.7563773044362474D+00
54 : real(dp), parameter :: B2_0 = 1.0173166567483392D-14
55 : real(dp), parameter :: E1_0 = -2.2308014220480969D+00
56 : real(dp), parameter :: F1_0 = 1.5176626709750911D-04
57 : real(dp), parameter :: E2_0 = -2.2350904778008243D+01
58 : real(dp), parameter :: F2_0 = 2.7741209323605414D-13
59 : real(dp), parameter :: T_0 = 7.9432823472428218D+07
60 : real(dp), parameter :: RHO_0 = 3.1622776917911558D+09
61 : real(dp), parameter :: r_0 = 2.2348420508311778D+20
62 : real(dp), parameter :: G1_0 = 1.5177849505266735D-04
63 : real(dp), parameter :: G2_0 = 2.8758525980353755D-13
64 : real(dp), parameter :: U_0 = 1.0723431204522564D+00
65 :
66 2 : real(dp) :: tmp, &
67 2 : B1_denom, dB1_denom_dRho, dB1_denom_dT, &
68 2 : B2_denom, dB2_denom_dRho, dB2_denom_dT
69 2 : real(dp) :: A1, dA1dT, B1_numerator, dB1_numerator_dT
70 2 : real(dp) :: A2, dA2dT, B2_numerator, dB2_numerator_dT
71 :
72 : include 'formats'
73 :
74 2 : R6=RHO*1d-6
75 2 : dR6dRho = 1d-6
76 2 : R6T=2d0*R6/UE
77 2 : dR6TdRho = 2d0*dR6dRho/UE
78 :
79 2 : R6T16=pow(R6T,1d0/6d0)
80 2 : dR6T16dRho = (1d0/6d0)*dR6TdRho*R6T16/R6T
81 2 : R6T13=R6T16*R6T16
82 2 : dR6T13dRho = 2*R6T16*dR6T16dRho
83 2 : T6=T*1d-6
84 2 : dT6dT=1d-6
85 2 : dT62dT=2d0*T6*dT6dT
86 2 : T613=pow(T6,1d0/3d0)
87 2 : dT613dT=(1d0/3d0)*dT6dT*T613/T6
88 2 : T623=T613*T613
89 2 : dT623dT=2*T613*dT613dT
90 2 : T653=T623*T6
91 2 : dT653dT = dT623dT*T6 + T623*dT6dT
92 :
93 2 : T62=T6*T6
94 2 : T612=sqrt(T6)
95 2 : dT612dT=0.5d0*dT6dT/T612
96 2 : T632=T6*T612
97 2 : dT632dT=1.5d0*T612*dT6dT
98 :
99 2 : U=1.35D0*R6T13/T623
100 2 : dUdT = -U * dT623dT / T623
101 2 : dUdRho = U * dR6T13dRho / R6T13
102 :
103 2 : U32 = U*sqrt(U)
104 2 : U52 = U*U32
105 :
106 2 : if (U < 1) then ! strong screening regime, eqn 4.8a in F&L
107 :
108 0 : A1 = pow2(1d0-4.222D-2*T623) + 2.643D-5*T653
109 0 : dA1dT = -2d0*4.222d-2*dT623dT*(1d0 - 4.222D-2*T623) + 2.643D-5*dT653dT
110 :
111 0 : B1_denom=A1*T623
112 0 : dB1_denom_dT = dA1dT*T623 + A1*dT623dT
113 :
114 0 : B1_numerator = 16.16D0*exp(-134.92d0/T613)
115 0 : dB1_numerator_dT = B1_numerator*134.92d0*dT613dT/(T613*T613)
116 :
117 0 : B1=B1_numerator/B1_denom
118 0 : dB1dT = dB1_numerator_dT/B1_denom - B1*dB1_denom_dT/B1_denom
119 :
120 0 : A2 = pow2(1d0-2.807D-2*T623) + 2.704D-6*T653
121 0 : dA2dT = -2*2.807D-2*dT623dT*(1d0-2.807D-2*T623) + 2.704D-6*dT653dT
122 :
123 0 : B2_denom=A2*T623
124 0 : dB2_denom_dT = dA2dT*T623 + A2*dT623dT
125 :
126 0 : B2_numerator = 244.6D0*pow5(1D0+3.528D-3*T623) * exp(-235.72D0/T613)
127 : dB2_numerator_dT = B2_numerator* &
128 0 : (5D0*3.528D-3*dT623dT/(1D0+3.528D-3*T623) + 235.72D0*dT613dT/T623)
129 :
130 0 : B2=B2_numerator/B2_denom
131 0 : dB2dT = dB2_numerator_dT/B2_denom - B2*dB2_denom_dT/B2_denom
132 :
133 0 : if (5.458D3 > R6T) then
134 :
135 0 : E1 = -1065.1D0/T6
136 0 : dE1dT = -E1 * dT6dT / T6
137 :
138 0 : F1 = exp(E1)/T632
139 0 : dF1dT = F1 * (dE1dT - dT632dT/T632)
140 :
141 0 : B1=B1+F1
142 0 : dB1dT = dB1dT + dF1dT
143 :
144 : end if
145 :
146 0 : if (1.836D4 > R6T) then
147 :
148 0 : E2 = -3336.4D0/T6
149 0 : dE2dT = -E2 * dT6dT / T6
150 :
151 0 : F2 = exp(E2)/T632
152 0 : dF2dT = F2 * (dE2dT - dT632dT/T632)
153 :
154 0 : B2=B2+F2
155 0 : dB2dT = dB2dT + dF2dT
156 :
157 : end if
158 :
159 0 : G1=B1*exp(60.492D0*R6T13/T6)
160 0 : dG1dT = G1*(dB1dT/B1 - 60.492D0*R6T13*dT6dT/(T6*T6))
161 0 : dG1dRho=0
162 :
163 0 : G2=B2*exp(106.35D0*R6T13/T6)
164 0 : dG2dT = G2*(dB2dT/B2 - 106.35D0*R6T13*dT6dT/(T6*T6))
165 0 : dG2dRho=0
166 :
167 : else ! pycnonuclear regime, eqn 4.8b in F&L
168 :
169 2 : AF=1d0/U32 + 1d0
170 2 : dAFdT = -1.5d0 * dUdT/U52
171 2 : dAFdRho = -1.5d0 * dUdRho/U52
172 :
173 2 : B1_denom=T612*(pow2(1.0d0-5.680D-2*R6T13)+8.815D-7*T62)
174 2 : dB1_denom_dT = B1_denom*dT612dT/T612 + T612*8.815D-7*dT62dT
175 2 : dB1_denom_dRho = -2*5.680D-2*(1.0d0-5.680D-2*R6T13)*T612*dR6T13dRho
176 :
177 2 : B1=1.178D0*AF*exp(-77.554d0/R6T16)/B1_denom
178 2 : dB1dT = B1 * (dAFdT/AF - dB1_denom_dT/B1_denom)
179 2 : dB1dRho = B1 * (dAFdRho/AF + 77.554d0*dR6T16dRho/(R6T16*R6T16) - dB1_denom_dRho/B1_denom)
180 :
181 2 : B2_denom=T612*(pow2(1.0d0-3.791D-2*R6T13)+5.162D-8*T62)
182 2 : dB2_denom_dT = B2_denom*dT612dT/T612 + T612*5.162D-8*dT62dT
183 :
184 2 : tmp = pow(Rho/UE,1d0/3d0)
185 2 : dB2_denom_dRho = T612*(-0.000252733d0 + 9.58112d-8*tmp)*tmp/Rho
186 :
187 2 : B2=13.48D0*AF*pow5(1.0d0+5.070D-3*R6T13)*exp(-135.08D0/R6T16)/B2_denom
188 2 : dB2dT = B2 * (dAFdT/AF - dB2_denom_dT/B2_denom)
189 2 : dB2dRho = B2 * (dAFdRho/AF + 135.08D0*dR6T16dRho/(R6T16*R6T16) - dB2_denom_dRho/B2_denom)
190 :
191 2 : if (5.458D3 > R6T) then
192 :
193 0 : E1 = (60.492d0*R6T13 - 1065.1D0)/T6
194 0 : dE1dT = -E1 * dT6dT / T6
195 0 : dE1dRho = 60.492d0*dR6T13dRho/T6
196 :
197 0 : F1 = exp(E1)/T632
198 0 : dF1dT = F1 * (dE1dT - dT632dT/T632)
199 0 : dF1dRho = F1 * dE1dRho
200 :
201 : !write(*,1) 'E1', E1
202 : !write(*,1) 'F1', F1
203 :
204 0 : G1=B1+F1
205 0 : dG1dT = dB1dT + dF1dT
206 0 : dG1dRho = dB1dRho + dF1dRho
207 :
208 : else
209 :
210 : G1=B1; dG1dRho = dB1dRho; dG1dT = dB1dT
211 :
212 : end if
213 :
214 2 : if (1.836D4 > R6T) then
215 :
216 2 : E2 = (106.35D0*R6T13 - 3336.4D0)/T6
217 2 : dE2dT = -E2 * dT6dT / T6
218 2 : dE2dRho = 106.35D0*dR6T13dRho/T6
219 :
220 2 : F2 = exp(E2)/T632
221 2 : dF2dT = F2 * (dE2dT - dT632dT/T632)
222 2 : dF2dRho = F2 * dE2dRho
223 :
224 : !write(*,1) 'E2', E2
225 : !write(*,1) 'F2', F2
226 :
227 2 : G2=B2+F2
228 2 : dG2dT = dB2dT + dF2dT
229 2 : dG2dRho = dB2dRho + dF2dRho
230 :
231 : else
232 :
233 : G2=B2; dG2dRho = dB2dRho; dG2dT = dB2dT
234 :
235 : end if
236 :
237 : end if
238 :
239 2 : r=5.120D29*G1*G2*Y*Y*Y*R6*R6 ! ergs/g/sec, eqn 4.7 in F&L
240 :
241 2 : if (r < 1d-99 .or. G1 < 1d-99 .or. G2 < 1d-99) then
242 0 : drdT = 0
243 0 : drdRho = 0
244 : else
245 2 : drdT = r * (dG1dT/G1 + dG2dT/G2)
246 2 : drdRho = r * (dG1dRho/G1 + dG2dRho/G2 + 2*dR6dRho/R6)
247 :
248 2 : return
249 :
250 : write(*,1) 'T', T
251 : write(*,1) 'RHO', RHO
252 : write(*,1) 'r', r
253 : write(*,1) 'G1', G1
254 : write(*,1) 'G2', G2
255 : write(*,1) 'U', U
256 : write(*,'(A)')
257 :
258 : write(*,1) 'abs(Rho_0 - Rho)', abs(Rho_0 - Rho)
259 :
260 : if (.true. .and. abs(Rho_0 - Rho) > 1d-2) then
261 : write(*,'(A)')
262 : write(*,1) 'analytic drdRho', drdRho
263 : write(*,1) 'numeric drdRho', (r_0 - r) / (Rho_0 - Rho)
264 : write(*,'(A)')
265 : write(*,1) 'analytic dG1dRho', dG1dRho
266 : write(*,1) 'numeric dG1dRho', (G1_0 - G1) / (Rho_0 - Rho)
267 : write(*,'(A)')
268 : write(*,1) 'analytic dG2dRho', dG2dRho
269 : write(*,1) 'numeric dG2dRho', (G2_0 - G2) / (Rho_0 - Rho)
270 : write(*,'(A)')
271 : write(*,1) 'analytic dUdRho', dUdRho
272 : write(*,1) 'numeric dUdRho', (U_0 - U) / (Rho_0 - Rho)
273 : write(*,'(A)')
274 : write(*,1) 'analytic AF', dAFdRho
275 : write(*,1) 'numeric AF', (AF_0 - AF) / (Rho_0 - Rho)
276 : write(*,'(A)')
277 : write(*,1) 'analytic B1_denom', dB1_denom_dRho
278 : write(*,1) 'numeric B1_denom', (B1_denom_0 - B1_denom) / (Rho_0 - Rho)
279 : write(*,'(A)')
280 : write(*,1) 'analytic B1', dB1dRho
281 : write(*,1) 'numeric B1', (B1_0 - B1) / (Rho_0 - Rho)
282 : write(*,'(A)')
283 : write(*,1) 'analytic B2_denom', dB2_denom_dRho
284 : write(*,1) 'numeric B2_denom', (B2_denom_0 - B2_denom) / (Rho_0 - Rho)
285 : write(*,'(A)')
286 : write(*,1) 'analytic B2', dB2dRho
287 : write(*,1) 'numeric B2', (B2_0 - B2) / (Rho_0 - Rho)
288 : write(*,'(A)')
289 : write(*,1) 'analytic E1', dE1dRho
290 : write(*,1) 'numeric E1', (E1_0 - E1) / (Rho_0 - Rho)
291 : write(*,'(A)')
292 : write(*,1) 'analytic F1', dF1dRho
293 : write(*,1) 'numeric F1', (F1_0 - F1) / (Rho_0 - Rho)
294 : write(*,'(A)')
295 : write(*,1) 'analytic E2', dE2dRho
296 : write(*,1) 'numeric E2', (E2_0 - E2) / (Rho_0 - Rho)
297 : write(*,'(A)')
298 : write(*,1) 'analytic F2', dF2dRho
299 : write(*,1) 'numeric F2', (F2_0 - F2) / (Rho_0 - Rho)
300 : write(*,'(A)')
301 : call mesa_error(__FILE__,__LINE__,'FL_epsnuc_3alf')
302 : end if
303 :
304 : end if
305 :
306 : end subroutine FL_epsnuc_3alf
307 :
308 : end module pycno
|
