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 atm_T_tau_relations
21 :
22 : use const_def, only: dp, iln10, two_thirds
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : private
29 : public :: get_T_tau_base
30 : public :: eval_T_tau
31 : public :: eval_T_tau_dq_dtau
32 :
33 : contains
34 :
35 9 : subroutine get_T_tau_base (id, tau_base, ierr)
36 :
37 : use atm_def, only: &
38 : ATM_T_TAU_EDDINGTON, &
39 : ATM_T_TAU_SOLAR_HOPF, &
40 : ATM_T_TAU_KRISHNA_SWAMY, &
41 : ATM_T_TAU_TRAMPEDACH_SOLAR
42 :
43 : integer, intent(in) :: id
44 : real(dp), intent(out) :: tau_base
45 : integer, intent(out) :: ierr
46 :
47 9 : ierr = 0
48 :
49 : ! Get the base optical depth
50 :
51 12 : select case (id)
52 : case (ATM_T_TAU_EDDINGTON)
53 3 : tau_base = 2._dp/3._dp
54 : case (ATM_T_TAU_SOLAR_HOPF)
55 2 : tau_base = 0.4116433502_dp
56 : case (ATM_T_TAU_KRISHNA_SWAMY)
57 2 : tau_base = 0.3121563_dp
58 : case (ATM_T_TAU_TRAMPEDACH_SOLAR)
59 2 : tau_base = 0.5147929058057147_dp
60 : case default
61 0 : write(*,*) 'Invalid id in get_T_tau_base: ', id
62 9 : call mesa_error(__FILE__,__LINE__)
63 : end select
64 :
65 9 : return
66 :
67 : end subroutine get_T_tau_base
68 :
69 :
70 6266 : subroutine eval_T_tau (id, tau, Teff, lnT, ierr)
71 :
72 : use atm_def, only: &
73 : ATM_T_TAU_EDDINGTON, &
74 : ATM_T_TAU_SOLAR_HOPF, &
75 : ATM_T_TAU_KRISHNA_SWAMY, &
76 : ATM_T_TAU_TRAMPEDACH_SOLAR
77 :
78 : integer, intent(in) :: id
79 : real(dp), intent(in) :: tau
80 : real(dp), intent(in) :: Teff
81 : real(dp), intent(out) :: lnT
82 : integer, intent(out) :: ierr
83 :
84 6266 : ierr = 0
85 :
86 : ! Evaluate the T-tau relation
87 :
88 9262 : select case (id)
89 : case (ATM_T_TAU_EDDINGTON)
90 4094 : call eval_Eddington(tau, Teff, lnT)
91 : case (ATM_T_TAU_SOLAR_HOPF)
92 2136 : call eval_solar_Hopf(tau, Teff, lnT)
93 : case (ATM_T_TAU_KRISHNA_SWAMY)
94 2172 : call eval_Krishna_Swamy(tau, Teff, lnT)
95 : case (ATM_T_TAU_TRAMPEDACH_SOLAR)
96 1134 : call eval_Trampedach_solar(tau, Teff, lnT)
97 : case default
98 0 : write(*,*) 'Invalid id in eval_T_tau: ', id
99 6266 : call mesa_error(__FILE__,__LINE__)
100 : end select
101 :
102 6266 : return
103 :
104 : end subroutine eval_T_tau
105 :
106 :
107 2996 : subroutine eval_Eddington (tau, Teff, lnT)
108 :
109 : real(dp), intent(in) :: tau
110 : real(dp), intent(in) :: Teff
111 : real(dp), intent(out) :: lnT
112 :
113 2996 : real(dp) :: Teff4
114 : real(dp) :: T4
115 :
116 : ! Evaluate the Eddington T-tau relation
117 :
118 2996 : Teff4 = Teff*Teff*Teff*Teff
119 2996 : T4 = 0.75d0*Teff4*(tau + two_thirds)
120 :
121 2996 : lnT = log(T4)*0.25d0
122 :
123 2996 : return
124 :
125 : end subroutine eval_Eddington
126 :
127 :
128 1098 : subroutine eval_solar_Hopf (tau, Teff, lnT)
129 :
130 : real(dp), intent(in) :: tau
131 : real(dp), intent(in) :: Teff
132 : real(dp), intent(out) :: lnT
133 :
134 : real(dp), parameter :: Q1 = 1.0361_dp
135 : real(dp), parameter :: Q2 = -0.3134_dp
136 : real(dp), parameter :: Q3 = 2.44799995_dp
137 : real(dp), parameter :: Q4 = -0.29589999_dp
138 : real(dp), parameter :: Q5 = 30._dp
139 :
140 1098 : real(dp) :: e1
141 1098 : real(dp) :: e2
142 1098 : real(dp) :: Teff4
143 1098 : real(dp) :: T4
144 :
145 : ! Evaluate the T-tau relation for an approximate Hopf function
146 : ! tuned to solar data; see MESA II, Sec. A.5. This is essentially
147 : ! equivalent to the fit given by Sonoi et al. (2019, A&A, 621,
148 : ! 84)
149 :
150 :
151 1098 : e1 = exp(-Q3*tau)
152 1098 : e2 = exp(-Q5*tau)
153 :
154 1098 : Teff4 = Teff*Teff*Teff*Teff
155 1098 : T4 = 0.75d0*Teff4*(tau + Q1 + Q2*e1 + Q4*e2)
156 :
157 1098 : lnT = log(T4)*0.25d0
158 :
159 1098 : return
160 :
161 : end subroutine eval_solar_Hopf
162 :
163 :
164 1038 : subroutine eval_Krishna_Swamy (tau, Teff, lnT)
165 :
166 : real(dp), intent(in) :: tau
167 : real(dp), intent(in) :: Teff
168 : real(dp), intent(out) :: lnT
169 :
170 : real(dp), parameter :: Q1 = 1.39_dp
171 : real(dp), parameter :: Q2 = -0.815_dp
172 : real(dp), parameter :: Q3 = 2.54_dp
173 : real(dp), parameter :: Q4 = -0.025_dp
174 : real(dp), parameter :: Q5 = 30._dp
175 :
176 1038 : real(dp) :: e1
177 1038 : real(dp) :: e2
178 1038 : real(dp) :: Teff4
179 1038 : real(dp) :: T4
180 :
181 : ! Evaluate the T-tau relation from Krishna-Swamy (1966, ApJ 145,
182 : ! 174–194)
183 :
184 1038 : e1 = exp(-Q3*tau)
185 1038 : e2 = exp(-Q5*tau)
186 :
187 1038 : Teff4 = Teff*Teff*Teff*Teff
188 1038 : T4 = 0.75d0*Teff4*(tau + Q1 + Q2*e1 + Q4*e2)
189 :
190 1038 : lnT = log(T4)*0.25d0
191 :
192 1038 : return
193 :
194 : end subroutine eval_Krishna_Swamy
195 :
196 :
197 1134 : subroutine eval_Trampedach_solar (tau, Teff, lnT)
198 :
199 : real(dp), intent(in) :: tau
200 : real(dp), intent(in) :: Teff
201 : real(dp), intent(out) :: lnT
202 :
203 : real(dp), parameter :: c0 = 0.6887302005929656_dp
204 : real(dp), parameter :: c1 = 0.0668697860833449_dp
205 : real(dp), parameter :: a = 0.9262126497691250_dp
206 : real(dp), parameter :: v = 0.7657856893402466_dp
207 : real(dp), parameter :: b = 0.1148742902769433_dp
208 :
209 : real(dp) :: x
210 1134 : real(dp) :: Teff4
211 1134 : real(dp) :: T4
212 :
213 : ! Evaluate the T-tau relation by Ball (2021, RNAAS 5, 7),
214 : ! which is a fit to the solar simulation by Trampedach
215 : ! et al. (2014, MNRAS 442, 805–820)
216 :
217 1134 : x = log10(tau)
218 :
219 1134 : if (x >= 0.07407427_dp) then
220 0 : write(*,*) 'WARNING: evaluating Trampedach_solar T-tau relation beyond valid region (log10(tau) < 0.0741):', x
221 : end if
222 :
223 1134 : Teff4 = Teff*Teff*Teff*Teff
224 1134 : T4 = 0.75d0*Teff4*(tau + c0 + c1*(x-b) + v*exp((x-a)/v))
225 :
226 1134 : lnT = log(T4)*0.25d0
227 :
228 1134 : return
229 :
230 : end subroutine eval_Trampedach_solar
231 :
232 :
233 0 : subroutine eval_T_tau_dq_dtau (id, tau, dq_dtau, ierr)
234 :
235 : use atm_def, only: &
236 : ATM_T_TAU_EDDINGTON, &
237 : ATM_T_TAU_SOLAR_HOPF, &
238 : ATM_T_TAU_KRISHNA_SWAMY, &
239 : ATM_T_TAU_TRAMPEDACH_SOLAR
240 :
241 : integer, intent(in) :: id
242 : real(dp), intent(in) :: tau
243 : real(dp), intent(out) :: dq_dtau
244 : integer, intent(out) :: ierr
245 :
246 0 : ierr = 0
247 :
248 : ! For a T(τ) relation of the standard form,
249 : !
250 : ! (T/Teff)⁴ = (3/4)[τ + q(τ)]
251 : !
252 : ! evaluate the derivative q'(τ).
253 :
254 0 : select case (id)
255 : case (ATM_T_TAU_EDDINGTON)
256 0 : call eval_Eddington_dq_dtau(tau, dq_dtau)
257 : case (ATM_T_TAU_SOLAR_HOPF)
258 0 : call eval_solar_Hopf_dq_dtau(tau, dq_dtau)
259 : case (ATM_T_TAU_KRISHNA_SWAMY)
260 0 : call eval_Krishna_Swamy_dq_dtau(tau, dq_dtau)
261 : case (ATM_T_TAU_TRAMPEDACH_SOLAR)
262 0 : call eval_Trampedach_solar_dq_dtau(tau, dq_dtau)
263 : case default
264 0 : write(*,*) 'Invalid id in eval_T_tau_dq_dtau: ', id
265 0 : call mesa_error(__FILE__,__LINE__)
266 : end select
267 :
268 0 : return
269 :
270 : end subroutine eval_T_tau_dq_dtau
271 :
272 :
273 0 : subroutine eval_Eddington_dq_dtau (tau, dq_dtau)
274 :
275 : real(dp), intent(in) :: tau
276 : real(dp), intent(out) :: dq_dtau
277 :
278 : ! Evaluate the Eddington q'(τ)
279 :
280 0 : dq_dtau = 0.0_dp
281 :
282 0 : return
283 :
284 : end subroutine eval_Eddington_dq_dtau
285 :
286 :
287 0 : subroutine eval_solar_Hopf_dq_dtau (tau, dq_dtau)
288 :
289 : real(dp), intent(in) :: tau
290 : real(dp), intent(out) :: dq_dtau
291 :
292 : real(dp), parameter :: Q1 = 1.0361_dp
293 : real(dp), parameter :: Q2 = -0.3134_dp
294 : real(dp), parameter :: Q3 = 2.44799995_dp
295 : real(dp), parameter :: Q4 = -0.29589999_dp
296 : real(dp), parameter :: Q5 = 30._dp
297 :
298 0 : real(dp) :: e1
299 0 : real(dp) :: e2
300 :
301 : ! Evaluate q'(τ) for an approximate Hopf function
302 : ! tuned to solar data; see MESA II, Sec. A.5. This is essentially
303 : ! equivalent to the fit given by Sonoi et al. (2019, A&A, 621,
304 : ! 84)
305 :
306 0 : e1 = exp(-Q3*tau)
307 0 : e2 = exp(-Q5*tau)
308 :
309 0 : dq_dtau = - Q2*Q3*e1 - Q4*Q5*e2
310 :
311 0 : return
312 :
313 : end subroutine eval_solar_Hopf_dq_dtau
314 :
315 :
316 0 : subroutine eval_Krishna_Swamy_dq_dtau (tau, dq_dtau)
317 :
318 : real(dp), intent(in) :: tau
319 : real(dp), intent(out) :: dq_dtau
320 :
321 : real(dp), parameter :: Q1 = 1.39_dp
322 : real(dp), parameter :: Q2 = -0.815_dp
323 : real(dp), parameter :: Q3 = 2.54_dp
324 : real(dp), parameter :: Q4 = -0.025_dp
325 : real(dp), parameter :: Q5 = 30._dp
326 :
327 0 : real(dp) :: e1
328 0 : real(dp) :: e2
329 :
330 : ! Evaluate q'(τ) for Krishna-Swamy (1966, ApJ 145, 174–194)
331 :
332 0 : e1 = exp(-Q3*tau)
333 0 : e2 = exp(-Q5*tau)
334 :
335 0 : dq_dtau = - Q2*Q3*e1 - Q4*Q5*e2
336 :
337 0 : return
338 :
339 : end subroutine eval_Krishna_Swamy_dq_dtau
340 :
341 :
342 0 : subroutine eval_Trampedach_solar_dq_dtau (tau, dq_dtau)
343 :
344 : real(dp), intent(in) :: tau
345 : real(dp), intent(out) :: dq_dtau
346 :
347 : real(dp), parameter :: c1 = 0.0668697860833449_dp
348 : real(dp), parameter :: a = 0.9262126497691250_dp
349 : real(dp), parameter :: v = 0.7657856893402466_dp
350 : real(dp), parameter :: b = 0.1148742902769433_dp
351 : real(dp), parameter :: w = 0.0514999047169869_dp
352 :
353 0 : real(dp) :: x
354 :
355 : ! Evaluate q'(τ) for Ball (2021, RNAAS 5, 7)
356 : ! using fit to solar simulation by
357 : ! Trampedach et al. (2014, MNRAS 442, 805–820)
358 :
359 0 : x = log10(tau)
360 :
361 0 : dq_dtau = (c1 + exp((x-a)/v))/(1._dp + exp((x-b)/w))/tau*iln10
362 :
363 0 : return
364 :
365 : end subroutine eval_Trampedach_solar_dq_dtau
366 :
367 : end module atm_T_tau_relations
|
