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 eval_weak
21 :
22 : use const_def, only: dp
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 : use rates_def
26 : use suzuki_tables
27 :
28 : implicit none
29 :
30 : contains
31 :
32 98224 : subroutine do_eval_weak_reaction_info( &
33 98224 : n, ids, reaction_ids, T9, YeRho, &
34 : eta, d_eta_dlnT, d_eta_dlnRho, &
35 98224 : lambda, dlambda_dlnT, dlambda_dlnRho, &
36 98224 : Q, dQ_dlnT, dQ_dlnRho, &
37 98224 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
38 : ierr)
39 : use const_def, only: ln10, kerg, mev_to_ergs
40 : use chem_def
41 : use interp_1d_def
42 : use utils_lib, only: is_bad
43 : integer, intent(in) :: n, ids(:), reaction_ids(:)
44 : real(dp), intent(in) :: T9, YeRho, eta, d_eta_dlnT, d_eta_dlnRho
45 : real(dp), dimension(:), intent(inout) :: &
46 : lambda, dlambda_dlnT, dlambda_dlnRho, &
47 : Q, dQ_dlnT, dQ_dlnRho, &
48 : Qneu, dQneu_dlnT, dQneu_dlnRho
49 : integer, intent(out) :: ierr
50 :
51 98224 : real(dp) :: alfa, beta, d_alfa_dlnT, alfa_hi_Z, beta_hi_Z, d_alfa_hi_Z_dlnT
52 : integer :: i, ir, cid
53 :
54 : include 'formats'
55 :
56 : call do_eval_weaklib_reaction_info( &
57 : n, ids, T9, YeRho, &
58 : eta, d_eta_dlnT, d_eta_dlnRho, &
59 : lambda, dlambda_dlnT, dlambda_dlnRho, &
60 : Q, dQ_dlnT, dQ_dlnRho, &
61 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
62 98224 : ierr)
63 98224 : if (ierr /= 0) then
64 : !write(*,*) 'failed in do_eval_weaklib_reaction_info'
65 : return
66 : end if
67 :
68 : !if (T9 >= max(T9_weaklib_full_on, T9_weaklib_full_on_hi_Z)) return
69 :
70 : ! revise lambda using rate for low T
71 : ! alfa is fraction from weaklib
72 98224 : if (T9 >= T9_weaklib_full_on) then
73 : alfa = 1d0
74 : d_alfa_dlnT = 0d0
75 72600 : else if (T9 > T9_weaklib_full_off) then
76 : alfa = (T9 - T9_weaklib_full_off)/ &
77 19284 : (T9_weaklib_full_on - T9_weaklib_full_off)
78 : d_alfa_dlnT = T9 / &
79 19284 : (T9_weaklib_full_on - T9_weaklib_full_off)
80 : else
81 : alfa = 0d0
82 : d_alfa_dlnT = 0d0
83 : end if
84 98224 : beta = 1d0 - alfa ! beta is fraction for low T
85 :
86 98224 : if (T9 >= T9_weaklib_full_on_hi_Z) then
87 : alfa_hi_Z = 1d0
88 : d_alfa_hi_Z_dlnT = 0d0
89 80008 : else if (T9 > T9_weaklib_full_off_hi_Z) then
90 : alfa_hi_Z = (T9 - T9_weaklib_full_off_hi_Z)/ &
91 8 : (T9_weaklib_full_on_hi_Z - T9_weaklib_full_off_hi_Z)
92 : d_alfa_hi_Z_dlnT = T9 / &
93 8 : (T9_weaklib_full_on_hi_Z - T9_weaklib_full_off_hi_Z)
94 : else
95 : alfa_hi_Z = 0d0
96 : d_alfa_hi_Z_dlnT = 0d0
97 : end if
98 98224 : beta_hi_Z = 1d0 - alfa_hi_Z
99 :
100 855618 : do i = 1, n
101 757394 : ir = reaction_ids(i)
102 757394 : if (ir == 0) cycle
103 757386 : cid = weak_reaction_info(1,ir)
104 757386 : if (weak_lowT_rate(ir) <= 0d0) cycle
105 27350 : if (cid <= 0) cycle
106 125574 : if (ids(i) <= 0) then
107 18216 : lambda(i) = weak_lowT_rate(ir)
108 18216 : dlambda_dlnT(i) = 0d0
109 18216 : dlambda_dlnRho(i) = 0d0
110 9134 : else if (chem_isos% Z(cid) >= weaklib_blend_hi_Z) then
111 9124 : lambda(i) = alfa_hi_Z*lambda(i) + beta_hi_Z*weak_lowT_rate(ir)
112 : dlambda_dlnT(i) = alfa_hi_Z*dlambda_dlnT(i) + &
113 9124 : d_alfa_hi_Z_dlnT * (lambda(i) - weak_lowT_rate(ir))
114 9124 : dlambda_dlnRho(i) = alfa_hi_Z*dlambda_dlnRho(i)
115 : else
116 10 : lambda(i) = alfa*lambda(i) + beta*weak_lowT_rate(ir)
117 : dlambda_dlnT(i) = alfa*dlambda_dlnT(i) + &
118 10 : d_alfa_dlnT * (lambda(i) - weak_lowT_rate(ir))
119 10 : dlambda_dlnRho(i) = alfa*dlambda_dlnRho(i)
120 : end if
121 : end do
122 :
123 98224 : end subroutine do_eval_weak_reaction_info
124 :
125 :
126 98224 : subroutine do_eval_weaklib_reaction_info( &
127 196448 : n, ids, T9_in, YeRho_in, &
128 : eta, d_eta_dlnT, d_eta_dlnRho, &
129 98224 : lambda, dlambda_dlnT, dlambda_dlnRho, &
130 98224 : Q, dQ_dlnT, dQ_dlnRho, &
131 98224 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
132 : ierr)
133 98224 : use const_def, only: ln10, kerg, mev_to_ergs
134 : use chem_def
135 : use interp_1d_def
136 : use utils_lib, only: is_bad, integer_dict_lookup
137 : integer, intent(in) :: n, ids(:)
138 : real(dp), intent(in) :: T9_in, YeRho_in, eta, d_eta_dlnT, d_eta_dlnRho
139 : real(dp), dimension(:), intent(inout) :: &
140 : lambda, dlambda_dlnT, dlambda_dlnRho, &
141 : Q, dQ_dlnT, dQ_dlnRho, &
142 : Qneu, dQneu_dlnT, dQneu_dlnRho
143 : integer, intent(out) :: ierr
144 :
145 : logical, parameter :: dbg = .false.
146 :
147 98224 : real(dp) :: T, T9, YeRho, lYeRho
148 : integer :: i, ir, in, out, rxn_idx
149 : logical :: neg
150 98224 : real(dp) :: Qx, conv, mue, d_mue_dlnRho, d_mue_dlnT
151 : character(len=iso_name_length) :: weak_lhs, weak_rhs
152 : integer, parameter :: nwork = pm_work_size
153 :
154 : real(dp) :: &
155 98224 : s_lambda, s_dlambda_dlnT, s_dlambda_dlnRho, &
156 98224 : s_Qneu, s_dQneu_dlnT, s_dQneu_dlnRho
157 :
158 : character(len=2*iso_name_length+1) :: key
159 :
160 : class(weak_rate_table), pointer :: table
161 :
162 : include 'formats'
163 :
164 98224 : ierr = 0
165 :
166 98224 : T9 = T9_in
167 98224 : YeRho = YeRho_in
168 98224 : lYeRho = log10(YeRho_in)
169 98224 : if (is_bad(lYeRho)) then
170 0 : ierr = -1
171 0 : return
172 :
173 : write(*,1) 'lYeRho', lYeRho
174 : write(*,1) 'YeRho_in', YeRho_in
175 : write(*,1) 'log10(YeRho_in)', log10(YeRho_in)
176 : !call mesa_error(__FILE__,__LINE__,'weak lYeRho')
177 : end if
178 :
179 98224 : if (n == 0) then
180 0 : write(*,*) 'problem in eval_weak_reaction_info: n == 0'
181 0 : write(*,2) 'n', n
182 0 : write(*,1) 'T9', T9
183 0 : return
184 : end if
185 :
186 855618 : do i = 1, n
187 :
188 757394 : lambda(i) = 0d0
189 757394 : dlambda_dlnT(i) = 0d0
190 757394 : dlambda_dlnRho(i) = 0d0
191 757394 : Q(i) = 0d0
192 757394 : dQ_dlnT(i) = 0d0
193 757394 : dQ_dlnRho(i) = 0d0
194 757394 : Qneu(i) = 0d0
195 757394 : dQneu_dlnT(i) = 0d0
196 757394 : dQneu_dlnRho(i) = 0d0
197 :
198 757394 : ir = ids(i)
199 757394 : if (ir <= 0) cycle
200 10060 : if (ir > size(weak_reactions_tables)) then
201 0 : call show_stuff
202 0 : write(*,*) 'bad ir', ir, i
203 0 : call mesa_error(__FILE__,__LINE__)
204 0 : ierr = -1
205 0 : return
206 : end if
207 :
208 10060 : table => weak_reactions_tables(ir) % t
209 :
210 10060 : T9 = T9_in
211 10060 : YeRho = YeRho_in
212 10060 : lYeRho = log10(YeRho_in)
213 :
214 : ! convert to MeV
215 10060 : conv = kerg/mev_to_ergs
216 10060 : T = T9*1d9
217 10060 : mue = eta*conv*T
218 10060 : d_mue_dlnRho = d_eta_dlnRho*conv*T
219 10060 : if (d_eta_dlnT == 0) then
220 10060 : d_mue_dlnT = 0
221 : else
222 : d_mue_dlnT = d_eta_dlnT*conv*T + mue
223 : end if
224 :
225 : ! clip small values to edge of table
226 10060 : if (T9 < table % T9s(1)) &
227 8 : T9 = table % T9s(1)
228 10060 : if (lYeRho < table % lYeRhos(1)) &
229 0 : lYeRho = table % lYeRhos(1)
230 :
231 : ! clip large values to edge of table
232 10060 : if (T9 > table % T9s(table % num_T9)) &
233 0 : T9 = table % T9s(table % num_T9)
234 10060 : if (lYeRho > table % lYeRhos(table % num_lYeRho)) &
235 2 : lYeRho = table % lYeRhos(table % num_lYeRho)
236 :
237 : call table% interpolate(T9, lYeRho, &
238 : lambda(i), dlambda_dlnT(i), dlambda_dlnRho(i), &
239 10060 : Qneu(i), dQneu_dlnT(i), dQneu_dlnRho(i), ierr)
240 :
241 10060 : in = weak_lhs_nuclide_id(ir)
242 10060 : out = weak_rhs_nuclide_id(ir)
243 10060 : Qx = chem_isos% mass_excess(in) - chem_isos% mass_excess(out)
244 :
245 10060 : if (use_suzuki_tables) then
246 : ! now, if there's a suzuki reaction, use that one instead
247 930 : call create_weak_dict_key(weak_lhs_nuclide_name(ir), weak_rhs_nuclide_name(ir), key)
248 930 : call integer_dict_lookup(suzuki_reactions_dict, key, rxn_idx, ierr)
249 930 : if (ierr /=0) then
250 : if (dbg) write(*,*) key, "is not a reaction included in the Suzuki tables"
251 788 : ierr = 0
252 : else
253 : if (dbg) write(*,*) key, "is a reaction included in the Suzuki tables"
254 142 : table => suzuki_reactions_tables(rxn_idx) %t
255 : call table % interpolate(T9, lYeRho, &
256 : s_lambda, s_dlambda_dlnT, s_dlambda_dlnRho, &
257 142 : s_Qneu, s_dQneu_dlnT, s_dQneu_dlnRho, ierr)
258 142 : if (ierr == 0) then
259 4 : lambda(i) = s_lambda
260 4 : dlambda_dlnT(i) = s_dlambda_dlnT
261 4 : dlambda_dlnRho(i) = s_dlambda_dlnRho
262 4 : Qneu(i) = s_Qneu
263 4 : dQneu_dlnT(i) = s_dQneu_dlnT
264 4 : dQneu_dlnRho(i) = s_dQneu_dlnRho
265 : !write(*,*) lYeRho, T9, s_lambda, s_Qneu
266 : end if
267 142 : ierr = 0
268 : end if
269 : end if
270 :
271 : ! neg is true for electron capture and positron emission
272 10060 : neg = ((chem_isos% Z(in) - chem_isos% Z(out)) == 1.0d0)
273 :
274 : ! in the past, these Q values used to include terms
275 : ! associated with the electron and ion chemical potentials
276 : ! these terms are now handled elsewhere, so Q is just the change in rest mass.
277 : ! since Qx is made from atomic mass excesses, it includes the electron rest mass.
278 :
279 10060 : if (neg) then ! electron capture and positron emission
280 9594 : Q(i) = Qx
281 9594 : dQ_dlnT(i) = 0
282 9594 : dQ_dlnRho(i) = 0
283 : else ! positron capture and electron emission
284 466 : Q(i) = Qx
285 466 : dQ_dlnT(i) = 0
286 466 : dQ_dlnRho(i) = 0
287 : end if
288 :
289 10060 : if (lambda(i) < 1d-30) then
290 268 : Qneu(i) = 0
291 268 : dQneu_dlnT(i) = 0
292 268 : dQneu_dlnRho(i) = 0
293 : end if
294 :
295 108284 : if (is_bad(Qneu(i))) then
296 0 : ierr = -1
297 0 : return
298 :
299 : write(*,2) 'lambda', i, lambda(i)
300 : write(*,2) 'Qneu', i, Qneu(i)
301 : call show_stuff
302 : end if
303 :
304 : end do
305 :
306 98224 : if (is_bad(lYeRho)) then
307 0 : ierr = -1
308 0 : return
309 : call show_stuff
310 : end if
311 :
312 :
313 : contains
314 :
315 0 : subroutine show_stuff
316 : integer :: i
317 : include 'formats'
318 0 : write(*,1) 'T9', T9
319 0 : write(*,1) 'lYeRho', lYeRho
320 0 : write(*,1) 'eta', eta
321 0 : write(*,1) 'mue', mue
322 0 : write(*,'(A)')
323 0 : do i = 1, n
324 0 : ir = ids(i)
325 0 : in = weak_lhs_nuclide_id(i)
326 0 : out = weak_rhs_nuclide_id(i)
327 0 : if (.true.) then
328 0 : write(*,4) 'ir, i, n', ir, i, n
329 0 : if (ir <= 0 .or. ir > ubound(weak_lhs_nuclide_id,dim=1)) cycle
330 0 : weak_lhs = chem_isos% name(weak_lhs_nuclide_id(ir))
331 0 : weak_rhs = chem_isos% name(weak_rhs_nuclide_id(ir))
332 0 : write(*,'(a30,3i5)') weak_lhs // weak_rhs, ir, i, n
333 : !write(*,1) 'chem_isos% mass_excess(in)', chem_isos% mass_excess(in)
334 : !write(*,1) 'chem_isos% mass_excess(out)', chem_isos% mass_excess(out)
335 0 : write(*,2) 'Qx', i, chem_isos% mass_excess(in) - chem_isos% mass_excess(out)
336 0 : write(*,2) 'Q', i, Q(i)
337 0 : write(*,2) 'Qneu', i, Qneu(i)
338 0 : write(*,'(A)')
339 : end if
340 : end do
341 0 : call mesa_error(__FILE__,__LINE__)
342 98224 : end subroutine show_stuff
343 :
344 : end subroutine do_eval_weaklib_reaction_info
345 :
346 : end module eval_weak
|
