Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2021 Josiah Schwab & 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_ecapture
21 :
22 : use rates_def
23 : use const_def, only: dp, ln2
24 : use math_lib
25 : use auto_diff
26 :
27 : implicit none
28 :
29 :
30 : contains
31 :
32 4 : subroutine do_eval_ecapture_reaction_info( &
33 4 : n, ids, cc, T9, YeRho, &
34 : etak, d_etak_dlnT, d_etak_dlnRho, & ! these are kinetic chemical potentials
35 4 : lambda, dlambda_dlnT, dlambda_dlnRho, &
36 4 : Q, dQ_dlnT, dQ_dlnRho, &
37 4 : Qneu, dQneu_dlnT, dQneu_dlnRho, &
38 : ierr)
39 : use rates_def, only: Coulomb_Info
40 : use eval_psi
41 : use eval_coulomb
42 : use const_def, only: ln10, kerg, mev_to_ergs, keV, me, clight, ev2erg, fine, pi
43 : use chem_def
44 : use utils_lib, only: is_bad, &
45 : integer_dict_lookup, integer_dict_size
46 :
47 : integer, intent(in) :: n, ids(:)
48 : type(Coulomb_Info), intent(in) :: cc
49 : real(dp), intent(in) :: T9, YeRho, etak, d_etak_dlnT, d_etak_dlnRho
50 : real(dp), dimension(:), intent(inout) :: &
51 : lambda, dlambda_dlnT, dlambda_dlnRho, &
52 : Q, dQ_dlnT, dQ_dlnRho, &
53 : Qneu, dQneu_dlnT, dQneu_dlnRho
54 : integer, intent(out) :: ierr
55 :
56 : logical, parameter :: dbg = .false.
57 :
58 : integer :: i, ir, j, lhs, rhs
59 : integer :: offset, ntrans, lo, hi
60 : integer :: offset_lhs, nstates_lhs, lo_lhs, hi_lhs
61 : integer :: offset_rhs, nstates_rhs, lo_rhs, hi_rhs
62 :
63 : type(auto_diff_real_2var_order1) :: mue, eta, beta, kT
64 :
65 : character(len=iso_name_length) :: ecapture_lhs, ecapture_rhs
66 :
67 4 : real(dp) :: Qx, mec2, Ei, Ef, G_GM, ln2ft
68 :
69 308 : real(dp), dimension(max_num_ecapture_states) :: E_lhs, E_rhs, J_lhs, J_rhs
70 : type(auto_diff_real_2var_order1), dimension(max_num_ecapture_states) :: bf, PPi
71 : type(auto_diff_real_2var_order1) :: Z, Eavg
72 :
73 : integer, dimension(max_num_ecapture_transitions) :: Si, Sf
74 104 : real(dp), dimension(max_num_ecapture_transitions) :: logft
75 :
76 : type(auto_diff_real_2var_order1), dimension(max_num_ecapture_transitions) :: zeta
77 : type(auto_diff_real_2var_order1) :: Ie_ad, Je_ad, lambda_ad, neutrino_ad, Qneu_ad
78 : type(auto_diff_real_2var_order1), dimension(max_num_ecapture_transitions) :: Pj
79 :
80 : ! for coulomb corrections
81 4 : real(dp) :: Z_lhs, Z_rhs
82 : type(auto_diff_real_2var_order1) :: muI_lhs, muI_rhs, delta_muI, delta_Q, Vs, f_muE
83 :
84 : character(len=2*iso_name_length+1) :: key
85 :
86 : integer :: rxn_type, ii
87 : integer, parameter :: rxn_ecapture = 1, rxn_betadecay = -1
88 :
89 : include 'formats'
90 :
91 4 : ierr = 0
92 :
93 : ! auto_diff variables have
94 : ! var1: lnT
95 : ! var2: lnRho
96 :
97 : ! first, translate the density, temperature, etc into appropriate units
98 :
99 4 : kT = 1d3 * keV * T9 ! in MeV
100 4 : kT% d1val1 = kT% val
101 4 : kT% d1val2 = 0d0
102 :
103 4 : mec2 = me * clight*clight / mev_to_ergs ! in MeV
104 4 : beta = mec2/kT ! dimesionless
105 :
106 : ! the chemical potentials from the equation of state are kinetic
107 : ! so add in the rest mass terms
108 :
109 4 : eta% val = etak + beta% val
110 4 : eta% d1val1 = d_etak_dlnT + beta % d1val1
111 4 : eta% d1val2 = d_etak_dlnRho + beta % d1val2
112 :
113 : ! only evaluate this for really degenerate stuff
114 4 : if (eta < 2*beta) return
115 :
116 : ! also need chemical potential in MeV
117 4 : mue = eta * kT
118 :
119 934 : do i = 1, n ! loop over reactions
120 :
121 : ! if there's not a weak reaction index, don't bother
122 930 : ir = ids(i)
123 930 : if (ir <= 0) then
124 : if (dbg) write(*,'(a,i3)') "No weak reaction for ", ir
125 : cycle
126 : end if
127 :
128 : ! get reactant names & ids
129 930 : ecapture_lhs = weak_lhs_nuclide_name(ir)
130 930 : ecapture_rhs = weak_rhs_nuclide_name(ir)
131 :
132 : ! generate reaction dictionary key (e.g. "mg24 na24")
133 930 : call create_ecapture_dict_key(ecapture_lhs, ecapture_rhs, key)
134 :
135 : ! get the transition data
136 930 : call integer_dict_lookup(ecapture_transitions_number_dict, key, ntrans, ierr)
137 930 : if (ierr /=0) then
138 : if (dbg) write(*,*) key, "is not a reaction included in ecapture module"
139 922 : ierr = 0
140 922 : cycle
141 : end if
142 : if (dbg) write(*,*) key, "is a reaction included in ecapture module"
143 :
144 8 : call integer_dict_lookup(ecapture_transitions_offset_dict, key, offset, ierr)
145 8 : if (ierr /=0) stop "ERROR: ecapture (transitions)"
146 : if (dbg) write(*,*) ntrans, offset
147 :
148 : ! get nuclide indices
149 8 : lhs = weak_lhs_nuclide_id(ir)
150 8 : rhs = weak_rhs_nuclide_id(ir)
151 :
152 : ! get species charges (which will be one different)
153 8 : Z_lhs = chem_isos% Z(lhs)
154 8 : Z_rhs = chem_isos% Z(rhs)
155 :
156 : ! determine whether this is a beta-decay or an electron capture
157 8 : rxn_type = int(Z_lhs - Z_rhs)
158 :
159 8 : Qx = 0
160 8 : G_GM = 1
161 4 : select case (rxn_type)
162 : case(rxn_ecapture)
163 : ! use the *atomic* isotope data to get the *nuclear* mass excess
164 4 : Qx = chem_isos% mass_excess(lhs) - chem_isos% mass_excess(rhs) - mec2
165 4 : G_GM = exp(pi* fine * Z_lhs)
166 4 : if (dbg) write(*,*) "ecapture ", key
167 : case(rxn_betadecay)
168 : ! use the *atomic* isotope data to get the *nuclear* mass excess
169 4 : Qx = chem_isos% mass_excess(lhs) - chem_isos% mass_excess(rhs) + mec2
170 4 : G_GM = exp(pi * fine * Z_rhs)
171 8 : if (dbg) write(*,*) "betadecay ", key
172 : end select
173 :
174 8 : lo = offset + 1
175 8 : hi = offset + ntrans
176 :
177 24 : Si(1:ntrans) = ecapture_transitions_data(lo:hi, i_Si)
178 24 : Sf(1:ntrans) = ecapture_transitions_data(lo:hi, i_Sf)
179 24 : logft(1:ntrans) = ecapture_logft_data(lo:hi)
180 :
181 : ! get the left state info (energies and spins)
182 :
183 8 : call integer_dict_lookup(ecapture_states_number_dict, ecapture_lhs, nstates_lhs, ierr)
184 8 : if (ierr /=0) then
185 0 : write(*,*) 'ecapture_lhs ' // trim(ecapture_lhs)
186 0 : write(*,*) 'ecapture_rhs ' // trim(ecapture_rhs)
187 0 : write(*,*) 'size of dict', integer_dict_size(ecapture_states_number_dict)
188 0 : stop "ERROR: ecapture_states_number_dict (lhs states)"
189 : end if
190 8 : call integer_dict_lookup(ecapture_states_offset_dict, ecapture_lhs, offset_lhs, ierr)
191 8 : if (ierr /=0) stop "ERROR: ecapture_states_offset_dict (lhs states)"
192 :
193 8 : lo_lhs = offset_lhs + 1
194 8 : hi_lhs = offset_lhs + nstates_lhs
195 :
196 40 : E_lhs(1:nstates_lhs) = ecapture_states_data(lo_lhs:hi_lhs, i_E)
197 40 : J_lhs(1:nstates_lhs) = ecapture_states_data(lo_lhs:hi_lhs, i_J)
198 :
199 : ! get the right state info (energies and spins)
200 :
201 8 : call integer_dict_lookup(ecapture_states_number_dict, ecapture_rhs, nstates_rhs, ierr)
202 8 : if (ierr /=0) stop "ERROR: ecapture_states_number_dict (rhs states)"
203 8 : call integer_dict_lookup(ecapture_states_offset_dict, ecapture_rhs, offset_rhs, ierr)
204 8 : if (ierr /=0) stop "ERROR: ecapture_states_offset_dict (rhs states)"
205 :
206 8 : lo_rhs = offset_rhs + 1
207 8 : hi_rhs = offset_rhs + nstates_rhs
208 :
209 40 : E_rhs(1:nstates_rhs) = ecapture_states_data(lo_rhs:hi_rhs, i_E)
210 : J_rhs(1:nstates_rhs) = ecapture_states_data(lo_rhs:hi_rhs, i_J)
211 :
212 : ! we assume the left hand side states have thermal occupation fractions
213 :
214 : ! calculate boltzmann factor (bf), partition function (Z), and <E>
215 8 : Z = 0d0
216 8 : Eavg = 0d0
217 40 : do ii=1,nstates_lhs
218 32 : bf(ii) = (2 * J_lhs(ii) + 1) * exp(-E_lhs(ii)/ kT)
219 32 : Z = Z + bf(ii)
220 40 : Eavg = Eavg + bf(ii) * E_lhs(ii)
221 : end do
222 8 : Eavg = Eavg / Z
223 :
224 : ! occupation probability
225 40 : do ii=1,nstates_lhs
226 40 : PPi(ii) = bf(ii)/Z
227 : end do
228 :
229 : ! now rearrange these rates so they apply to the transitions
230 24 : do j = 1, ntrans
231 24 : Pj(j) = PPi(Si(j))
232 : end do
233 :
234 : if (dbg) then
235 : do j = 1, ntrans
236 : write(*,"(2(F5.2, I2, F5.2))") E_lhs(Si(j)), int(J_lhs(Si(j))), Pj(j), &
237 : E_rhs(Sf(j)), int(J_rhs(Sf(j))), logft(j)
238 : end do
239 : end if
240 :
241 : ! get chemical potential difference, related to coulomb correction
242 :
243 : ! get ion chemical potentials (already in units of kT)
244 8 : muI_lhs = do_muI_coulomb(cc, Z_lhs)
245 8 : muI_rhs = do_muI_coulomb(cc, Z_rhs)
246 :
247 : ! shift should be negative (for e-capture) given our definitions
248 8 : delta_muI = muI_lhs - muI_rhs
249 8 : delta_Q = delta_muI * kT
250 :
251 : ! get screening potential (in fraction of fermi energy)
252 8 : Vs = do_Vs_coulomb(cc, Z_lhs)
253 8 : f_muE = 1d0 - Vs
254 :
255 : if (dbg) write(*,*) eta, muI_lhs, muI_rhs, delta_muI, Vs
256 :
257 8 : lambda_ad = 0d0
258 8 : neutrino_ad = 0d0
259 :
260 : ! do the phase space integrals
261 24 : do j = 1, ntrans
262 :
263 16 : Ei = E_lhs(Si(j))
264 16 : Ef = E_rhs(Sf(j))
265 :
266 : ! calculate energy difference, including electron rest mass
267 16 : zeta(j) = (Qx + delta_Q - Ef + Ei)/ kT
268 :
269 8 : select case(rxn_type)
270 : case(rxn_ecapture)
271 8 : call do_psi_Iec_and_Jec(beta, zeta(j), f_muE*eta, Ie_ad, Je_ad)
272 : case(rxn_betadecay)
273 16 : call do_psi_Iee_and_Jee(beta, zeta(j), f_muE*eta, Ie_ad, Je_ad)
274 : end select
275 :
276 : ! apply fermi-function correction factor
277 16 : Ie_ad = G_GM * Ie_ad
278 16 : Je_ad = G_GM * Je_ad
279 :
280 : ! convert to rates
281 16 : ln2ft = ln2 * exp10(-logft(j))
282 : ! protect against 0s
283 24 : if ((Ie_ad > 0) .and. (Je_ad > 0)) then
284 16 : lambda_ad = lambda_ad + Ie_ad * ln2ft * Pj(j)
285 16 : neutrino_ad = neutrino_ad + mec2 * Je_ad * ln2ft * Pj(j)
286 : end if
287 :
288 : end do
289 :
290 8 : if (lambda_ad > 1d-30) then
291 4 : Qneu_ad = neutrino_ad / lambda_ad
292 : else
293 4 : Qneu_ad = 0d0
294 : end if
295 :
296 : ! unpack from auto_diff
297 8 : lambda(i) = lambda_ad % val
298 8 : dlambda_dlnT(i) = lambda_ad % d1val1
299 8 : dlambda_dlnRho(i) = lambda_ad % d1val2
300 :
301 8 : Qneu(i) = Qneu_ad% val
302 8 : dQneu_dlnT(i) = Qneu_ad% d1val1
303 8 : dQneu_dlnRho(i) = Qneu_ad% d1val2
304 :
305 : ! this is the *total* energy per decay (neu losses are subtracted later)
306 :
307 : ! in the past, these Q values used to include terms associated
308 : ! with the electron and ion chemical potentials.
309 : ! these terms are now handled elsewhere, so Q is just the change in rest mass.
310 : ! note that Qx is defined differently here than in eval_weak,
311 : ! which is why we explicitly add/subtract mec2.
312 :
313 974 : select case(rxn_type)
314 : case(rxn_ecapture)
315 :
316 4 : Q(i) = Qx + mec2
317 4 : dQ_dlnT(i) = 0
318 4 : dQ_dlnRho(i) = 0
319 :
320 : case(rxn_betadecay)
321 :
322 4 : Q(i) = Qx - mec2
323 4 : dQ_dlnT(i) = 0
324 8 : dQ_dlnRho(i) = 0
325 :
326 : end select
327 :
328 : end do
329 :
330 4 : ierr = 0
331 :
332 4 : end subroutine do_eval_ecapture_reaction_info
333 :
334 : end module eval_ecapture
|
