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 : ! In general, this routine follows the notation in Schwab et al. (2015)
21 :
22 : module eval_psi
23 :
24 : use const_def, only: dp
25 : use math_lib
26 :
27 : implicit none
28 :
29 : contains
30 :
31 120 : subroutine fd_integral(dk, y, F)
32 : ! wrap eos_fermi_dirac_integral for auto_diff
33 : use auto_diff
34 : use eos_lib, only: eos_fermi_dirac_integral
35 : real(dp), intent(in) :: dk
36 : type(auto_diff_real_2var_order1), intent(in) :: y
37 : type(auto_diff_real_2var_order1), intent(out) :: F
38 : real(dp) :: Fval, dF_dy, dF_dz
39 112 : call eos_fermi_dirac_integral(dk, y% val, 0.0_dp, Fval, dF_dy, dF_dz)
40 112 : F = Fval
41 112 : F% d1val1 = dF_dy * y% d1val1
42 112 : F% d1val2 = dF_dy * y% d1val2
43 112 : end subroutine fd_integral
44 :
45 :
46 8 : subroutine do_psi_Iec_and_Jec(beta, zeta, eta, I, J)
47 :
48 112 : use auto_diff
49 : ! calculate the phase space integral for electron emission (beta-decay)
50 :
51 :
52 : ! auto_diff variables have
53 : ! var1: lnT
54 : ! var2: lnRho
55 :
56 : type(auto_diff_real_2var_order1), intent(in) :: beta ! mec2 / kT
57 : type(auto_diff_real_2var_order1), intent(in) :: zeta ! Q_n / kT
58 : type(auto_diff_real_2var_order1), intent(in) :: eta ! chemical potential / kT
59 :
60 : type(auto_diff_real_2var_order1), intent(out) :: I, J ! phase space integral
61 :
62 : type(auto_diff_real_2var_order1) :: y
63 :
64 : type(auto_diff_real_2var_order1) :: F2, F3, F4, F5
65 : type(auto_diff_real_2var_order1) :: c2, c3
66 :
67 : ! check that assumptions are met
68 8 : if (zeta>-beta) stop "ECAPTURE: zeta > -beta"
69 :
70 8 : y = zeta+eta
71 :
72 8 : call fd_integral(5.0_dp, y, F5)
73 8 : call fd_integral(4.0_dp, y, F4)
74 8 : call fd_integral(3.0_dp, y, F3)
75 8 : call fd_integral(2.0_dp, y, F2)
76 :
77 8 : c3 = -2.0_dp*zeta
78 8 : c2 = zeta*zeta
79 :
80 8 : I = F4 + c3*F3 + c2*F2
81 8 : J = F5 + c3*F4 + c2*F3
82 :
83 8 : I = I / pow5(beta)
84 8 : J = J / pow6(beta)
85 :
86 8 : return
87 :
88 8 : end subroutine do_psi_Iec_and_Jec
89 :
90 :
91 8 : subroutine do_psi_Iee_and_Jee(beta, zeta, eta, I, J)
92 :
93 8 : use auto_diff
94 :
95 : ! calculate the phase space integral for electron emission (beta-decay)
96 :
97 :
98 : ! auto_diff variables have
99 : ! var1: lnT
100 : ! var2: lnRho
101 :
102 : type(auto_diff_real_2var_order1), intent(in) :: beta ! mec2 / kT
103 : type(auto_diff_real_2var_order1), intent(in) :: zeta ! Q_n / kT
104 : type(auto_diff_real_2var_order1), intent(in) :: eta ! chemical potential / kT
105 :
106 : type(auto_diff_real_2var_order1), intent(out) :: I, J ! phase space integral
107 :
108 : type(auto_diff_real_2var_order1) :: y
109 :
110 : type(auto_diff_real_2var_order1) :: F0, F1, F2, F3, F4, F5
111 : type(auto_diff_real_2var_order1) :: c0, c1, c2, c3, c4
112 :
113 : ! check that assumptions are met
114 8 : if (zeta<beta) stop "ECAPTURE: zeta < beta"
115 :
116 :
117 8 : y = zeta-eta
118 :
119 8 : call fd_integral(5.0_dp, y, F5)
120 8 : call fd_integral(4.0_dp, y, F4)
121 8 : call fd_integral(3.0_dp, y, F3)
122 8 : call fd_integral(2.0_dp, y, F2)
123 :
124 8 : c3 = -2.0_dp*zeta
125 8 : c2 = zeta*zeta
126 :
127 8 : I = F4 + c3*F3 + c2*F2
128 8 : J = F5 + c3*F4 + c2*F3
129 :
130 :
131 8 : y = beta-eta
132 :
133 8 : call fd_integral(5.0_dp, y, F5)
134 8 : call fd_integral(4.0_dp, y, F4)
135 8 : call fd_integral(3.0_dp, y, F3)
136 8 : call fd_integral(2.0_dp, y, F2)
137 8 : call fd_integral(1.0_dp, y, F1)
138 8 : call fd_integral(0.0_dp, y, F0)
139 :
140 8 : c3 = (2.0_dp*zeta-4.0_dp*beta)
141 8 : c2 = (6.0_dp*beta*beta - 6.0_dp*beta*zeta + zeta*zeta)
142 8 : c1 = -2.0_dp*beta*(2.0_dp*Beta*beta - 3.0_dp*beta*zeta + zeta*zeta)
143 8 : c0 = beta*beta*(beta-zeta)*(beta-zeta)
144 :
145 8 : I = I - (F4 + c3*F3 + c2*F2 + c1*F1 + c0*F0)
146 :
147 8 : c4 = (3.0_dp*zeta-5.0_dp*beta)
148 8 : c3 = (10.0_dp*beta*beta-12.0_dp*beta*zeta+3.0_dp*zeta*zeta)
149 8 : c2 = (zeta*zeta*zeta - 9.0_dp*beta*zeta*zeta + 18.0_dp*beta*beta*zeta - 10.0_dp*beta*beta*beta)
150 8 : c1 = beta*(5.0_dp*beta - 2.0_dp*zeta)*(beta-zeta)*(beta-zeta)
151 8 : c0 = -beta*beta*pow3(beta-zeta)
152 :
153 8 : J = J - (F5 + c4*F4 + c3*F3 + c2*F2 + c1*F1 + c0*F0)
154 :
155 :
156 8 : I = I / pow5(beta)
157 8 : J = J / pow6(beta)
158 :
159 8 : return
160 :
161 8 : end subroutine do_psi_Iee_and_Jee
162 :
163 : end module eval_psi
|
