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 kap_eval
21 :
22 : use utils_lib, only: is_bad, mesa_error
23 : use kap_def
24 : use const_def, only: dp, ln10, sige, me, kev, kerg, amu, clight
25 : use math_lib
26 :
27 : implicit none
28 :
29 : private
30 : public :: get_kap_results
31 : public :: combine_rad_with_conduction
32 : public :: get_kap_results_blend_t
33 : public :: compton_opacity
34 :
35 : contains
36 :
37 86216 : subroutine get_kap_Results( &
38 : rq, zbar, X, Z, XC, XN, XO, XNe, logRho, logT, &
39 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
40 : eta, d_eta_dlnRho, d_eta_dlnT, &
41 : kap_fracs, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
42 : use utils_lib, only: is_bad
43 : use auto_diff
44 :
45 : type (Kap_General_Info), pointer :: rq
46 : real(dp), intent(in) :: X, XC, XN, XO, XNe, Z, zbar
47 : real(dp), intent(in) :: logRho, logT
48 : real(dp), intent(in) :: lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
49 : real(dp), intent(in) :: eta, d_eta_dlnRho, d_eta_dlnT
50 : real(dp), intent(out) :: kap_fracs(num_kap_fracs)
51 : real(dp), intent(out) :: kap ! opacity
52 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
53 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
54 : integer, intent(out) :: ierr ! 0 means AOK.
55 :
56 86216 : real(dp) :: kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT
57 86216 : real(dp) :: kap_compton, dlnkap_compton_dlnRho, dlnkap_compton_dlnT
58 :
59 86216 : real(dp) :: logT_Compton_blend_lo, logT_Compton_blend_hi
60 86216 : real(dp) :: logR_Compton_blend_lo, logR_Compton_blend_hi
61 :
62 : real(dp) :: Rho, T
63 : type(auto_diff_real_2var_order1) :: logT_auto, logRho_auto, logR_auto
64 : type(auto_diff_real_2var_order1) :: logkap_rad, logkap_compton
65 : type(auto_diff_real_2var_order1) :: blend_logT, blend_logR, blend
66 :
67 86216 : real(dp) :: frac_Type2, frac_highT, frac_lowT
68 :
69 : logical :: dbg
70 :
71 : include 'formats'
72 :
73 86216 : dbg = rq% dbg
74 86216 : if (dbg) dbg = & ! check limits
75 : logT >= rq% logT_lo .and. logT <= rq% logT_hi .and. &
76 : logRho >= rq% logRho_lo .and. logRho <= rq% logRho_hi .and. &
77 : X >= rq% X_lo .and. X <= rq% X_hi .and. &
78 0 : Z >= rq% Z_lo .and. Z <= rq% Z_hi
79 :
80 : ! auto diff
81 : ! var1: logRho
82 : ! var2: logT
83 :
84 86216 : Rho = exp10(logRho)
85 86216 : logRho_auto% val = logRho
86 86216 : logRho_auto% d1val1 = 1d0
87 86216 : logRho_auto% d1val2 = 0d0
88 :
89 86216 : T = exp10(logT)
90 86216 : logT_auto% val = logT
91 86216 : logT_auto% d1val1 = 0d0
92 86216 : logT_auto% d1val2 = 1d0
93 :
94 86216 : logR_auto = logRho_auto - 3d0*logT_auto + 18d0
95 :
96 86216 : if (dbg) write(*,1) 'logRho', logRho
97 0 : if (dbg) write(*,1) 'logT', logT
98 86216 : if (dbg) write(*,1) 'logR', logR_auto% val
99 :
100 : ! initialize fracs
101 86216 : kap_fracs = 0
102 86216 : frac_Type2 = 0d0
103 86216 : frac_lowT = 0d0
104 86216 : frac_highT = 0d0
105 :
106 : ! blend to Compton at high T
107 86216 : logT_Compton_blend_hi = rq% logT_Compton_blend_hi
108 86216 : logT_Compton_blend_lo = logT_Compton_blend_hi - 0.50d0
109 :
110 : ! blend in logT
111 86216 : if (logT_auto >= logT_Compton_blend_hi) then
112 0 : blend_logT = 1d0
113 86216 : else if (logT_auto > logT_Compton_blend_lo .and. logT_auto <= logT_Compton_blend_hi) then
114 0 : blend_logT = (logT_auto - logT_Compton_blend_lo) / (logT_Compton_blend_hi - logT_Compton_blend_lo)
115 : else
116 86216 : blend_logT = 0d0
117 : end if
118 : ! quintic smoothing
119 86216 : blend_logT = -blend_logT*blend_logT*blend_logT*(-10d0 + blend_logT*(15d0 - 6d0*blend_logT))
120 :
121 :
122 : ! blend to Compton at low R
123 86216 : logR_Compton_blend_lo = rq% logR_Compton_blend_lo
124 86216 : logR_Compton_blend_hi = logR_Compton_blend_lo + 0.50d0
125 :
126 : ! blend in logR
127 86216 : if (logR_auto <= logR_Compton_blend_lo) then
128 0 : blend_logR = 1d0
129 86216 : else if (logR_auto > logR_Compton_blend_lo .and. logR_auto <= logR_Compton_blend_hi) then
130 0 : blend_logR = (logR_Compton_blend_hi - logR_auto) / (logR_Compton_blend_hi - logR_Compton_blend_lo)
131 : else
132 86216 : blend_logR = 0d0
133 : end if
134 : ! quintic smoothing
135 86216 : blend_logR = -blend_logR*blend_logR*blend_logR*(-10d0 + blend_logR*(15d0 - 6d0*blend_logR))
136 :
137 :
138 : ! smoothly combine blends
139 86216 : blend = 1d0 - (1d0-blend_logT)*(1d0-blend_logR)
140 86216 : kap_fracs(i_frac_Compton) = blend% val
141 :
142 :
143 86216 : if (blend > 0) then ! at least some compton
144 :
145 0 : if (rq % use_other_compton_opacity) then
146 : call rq% other_compton_opacity(rq% handle, Rho, T, &
147 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
148 : eta, d_eta_dlnRho, d_eta_dlnT, &
149 0 : kap_compton, dlnkap_compton_dlnRho, dlnkap_compton_dlnT, ierr)
150 : else
151 : call Compton_Opacity(rq, Rho, T, &
152 : lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
153 : eta, d_eta_dlnRho, d_eta_dlnT, &
154 0 : kap_compton, dlnkap_compton_dlnRho, dlnkap_compton_dlnT, ierr)
155 : end if
156 :
157 0 : if (ierr /= 0) return
158 :
159 : else ! no Compton
160 :
161 86216 : kap_compton = 1d-30
162 86216 : dlnkap_compton_dlnRho = 0d0
163 86216 : dlnkap_compton_dlnT = 0d0
164 :
165 : end if
166 :
167 : ! pack into auto_diff type
168 86216 : logkap_compton = log10(kap_compton)
169 86216 : logkap_compton% d1val1 = dlnkap_compton_dlnRho
170 86216 : logkap_compton% d1val2 = dlnkap_compton_dlnT
171 :
172 :
173 86216 : if (blend < 1) then ! at least some tables
174 :
175 86216 : if (rq% use_other_radiative_opacity) then
176 : call rq% other_radiative_opacity( &
177 : rq% handle, X, Z, XC, XN, XO, XNe, logRho, logT, &
178 : frac_lowT, frac_highT, frac_Type2, &
179 0 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, ierr)
180 : else
181 : call Get_kap_Results_blend_T( &
182 : rq, X, Z, XC, XN, XO, XNe, logRho, logT, &
183 : frac_lowT, frac_highT, frac_Type2, &
184 86216 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, ierr)
185 : end if
186 86216 : if (ierr /= 0) return
187 :
188 : ! revise reported fractions based on Compton blend
189 86216 : frac_lowT = (1d0 - blend% val) * frac_lowT
190 86216 : frac_highT = (1d0 - blend% val) * frac_highT
191 : ! the value of frac_Type2 doesn't need revised if it represents
192 : ! the fraction of highT opacities provided by the Type2 tables
193 : ! frac_Type2 = (1d0 - blend% val) * frac_Type2
194 :
195 : else ! no tables
196 :
197 0 : kap_rad = 1d-30
198 0 : dlnkap_rad_dlnRho = 0d0
199 0 : dlnkap_rad_dlnT = 0d0
200 :
201 : end if
202 :
203 : ! pack into auto_diff type
204 86216 : logkap_rad = log10(kap_rad)
205 86216 : logkap_rad% d1val1 = dlnkap_rad_dlnRho
206 86216 : logkap_rad% d1val2 = dlnkap_rad_dlnT
207 :
208 : ! pack kap_fracs from tables
209 86216 : kap_fracs(i_frac_lowT) = frac_lowT
210 86216 : kap_fracs(i_frac_highT) = frac_highT
211 86216 : kap_fracs(i_frac_Type2) = frac_Type2
212 :
213 : ! do blend
214 86216 : logkap_rad = blend * logkap_compton + (1d0 - blend) * logkap_rad
215 :
216 : ! unpack auto_diff
217 86216 : kap_rad = exp10(logkap_rad% val)
218 86216 : dlnkap_rad_dlnRho = logkap_rad% d1val1
219 86216 : dlnkap_rad_dlnT = logkap_rad% d1val2
220 :
221 : call combine_rad_with_conduction( &
222 : rq, logRho, logT, zbar, &
223 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, &
224 86216 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
225 :
226 86216 : end subroutine get_kap_Results
227 :
228 :
229 86416 : subroutine get_kap_Results_blend_T( &
230 : rq, X, Z, XC, XN, XO, XNe, logRho_in, logT_in, &
231 : frac_lowT, frac_highT, frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
232 :
233 86216 : use kap_def
234 : use utils_lib, only: is_bad
235 :
236 : ! INPUT
237 : type (Kap_General_Info), pointer :: rq
238 : real(dp), intent(in) :: X, Z, XC, XN, XO, XNe ! composition
239 : real(dp), intent(in) :: logRho_in ! density
240 : real(dp), intent(in) :: logT_in ! temperature
241 :
242 : ! OUTPUT
243 : real(dp), intent(out) :: frac_lowT, frac_highT, frac_Type2
244 : real(dp), intent(out) :: kap ! opacity
245 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
246 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
247 : integer, intent(out) :: ierr ! 0 means AOK.
248 :
249 86216 : real(dp) :: logRho, Rho, logT, T
250 : real(dp) :: frac_Type1, dZ, frac_X, frac_dZ, &
251 : dXC, dXO, fC, fN, fO, fNe, fHeavy, ZHeavy, dXsum
252 :
253 86216 : real(dp) :: lowT_logT_max, logT_min
254 :
255 : logical :: clipped_Rho
256 :
257 : logical :: dbg
258 :
259 86216 : real(dp) :: alfa, beta, &
260 86216 : logKap, logKap1, logKap2, &
261 86216 : kap1, dlnkap1_dlnRho, dlnkap1_dlnT, &
262 86216 : kap2, dlnkap2_dlnRho, dlnkap2_dlnT, &
263 86216 : lower_bdy, upper_bdy, alfa0, d_alfa0_dlnT, &
264 86216 : d_alfa_dlnT, d_beta_dlnT, dkap_dlnT
265 :
266 :
267 : include 'formats'
268 :
269 86216 : dbg = .false.
270 :
271 86216 : logRho = logRho_in; Rho = exp10(logRho)
272 86216 : logT = logT_in; T = exp10(logT)
273 :
274 : if (dbg) write(*,1) 'Get_kap_blend_logT logT', logT
275 :
276 86216 : clipped_Rho = .false.
277 86216 : if (logRho < kap_min_logRho) then
278 0 : logRho = kap_min_logRho
279 0 : rho = exp10(kap_min_logRho)
280 0 : clipped_Rho = .true.
281 : end if
282 :
283 86216 : frac_lowT = 0d0
284 86216 : frac_highT = 0d0
285 : frac_Type2 = 0d0
286 :
287 86216 : if (rq% use_Type2_opacities .and. &
288 : associated(kap_co_z_tables(rq% kap_CO_option)% ar)) then
289 86214 : logT_min = kap_co_z_tables(rq% kap_CO_option)% ar(1)% x_tables(1)% logT_min
290 : else
291 2 : logT_min = kap_z_tables(rq% kap_option)% ar(1)% x_tables(1)% logT_min
292 : end if
293 :
294 86216 : lower_bdy = max(rq% kap_blend_logT_lower_bdy, logT_min)
295 : if (dbg) write(*,1) 'Get_kap_blend_logT lower_bdy', lower_bdy, &
296 : rq% kap_blend_logT_lower_bdy, logT_min
297 86216 : if (logT <= lower_bdy) then ! all lowT
298 : if (dbg) write(*,1) 'all lowT'
299 : call Get_kap_lowT_Results(rq, &
300 : X, Z, XC, XN, XO, XNe, logRho, logT, &
301 5852 : frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
302 5852 : if (clipped_Rho) dlnkap_dlnRho = 0d0
303 5852 : frac_lowT = 1d0
304 86116 : return
305 : end if
306 :
307 :
308 80364 : select case (rq% kap_lowT_option)
309 : case (kap_lowT_kapCN)
310 0 : lowT_logT_max = kapCN_max_logT
311 : case (kap_lowT_AESOPUS)
312 0 : lowT_logT_max = kA % max_logT
313 : case default
314 80364 : lowT_logT_max = kap_lowT_z_tables(rq% kap_lowT_option)% ar(1)% x_tables(1)% logT_max
315 : end select
316 :
317 :
318 80364 : upper_bdy = min(rq% kap_blend_logT_upper_bdy, lowT_logT_max)
319 : if (dbg) write(*,1) 'Get_kap_blend_logT upper_bdy', upper_bdy, &
320 : rq% kap_blend_logT_upper_bdy, lowT_logT_max
321 80364 : if (logT >= upper_bdy) then ! no lowT
322 : if (dbg) write(*,1) 'no lowT'
323 : call Get_kap_highT_Results(rq, &
324 : X, Z, XC, XN, XO, XNe, logRho, logT, &
325 80264 : frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
326 80264 : if (clipped_Rho) dlnkap_dlnRho = 0d0
327 80264 : frac_highT = 1d0
328 80264 : return
329 : end if
330 :
331 : ! in blend region
332 :
333 : call Get_kap_lowT_Results(rq, &
334 : X, Z, XC, XN, XO, XNe, logRho, logT, &
335 100 : frac_Type2, kap2, dlnkap2_dlnRho, dlnkap2_dlnT, ierr)
336 100 : if (clipped_Rho) dlnkap2_dlnRho = 0d0
337 100 : if (ierr /= 0) return
338 :
339 : call Get_kap_highT_Results(rq, &
340 : X, Z, XC, XN, XO, XNe, logRho, logT, &
341 100 : frac_Type2, kap1, dlnkap1_dlnRho, dlnkap1_dlnT, ierr)
342 100 : if (clipped_Rho) dlnkap1_dlnRho = 0d0
343 100 : if (ierr /= 0) return
344 :
345 100 : alfa0 = (logT - lower_bdy) / (upper_bdy - lower_bdy)
346 100 : d_alfa0_dlnT = 1d0/(upper_bdy - lower_bdy)/ln10
347 :
348 : ! must smooth the transitions near alfa = 0.0 and 1.0
349 : ! Rich Townsend's smoothing function for this
350 100 : alfa = -alfa0*alfa0*alfa0*(-10d0 + alfa0*(15d0 - 6d0*alfa0))
351 100 : d_alfa_dlnT = 30d0*(alfa0 - 1d0)*(alfa0 - 1d0)*alfa0*alfa0*d_alfa0_dlnT
352 :
353 100 : beta = 1d0 - alfa
354 100 : d_beta_dlnT = -d_alfa_dlnT
355 :
356 :
357 100 : logKap1 = log10(kap1); logKap2 = log10(kap2)
358 :
359 100 : logKap = alfa*logKap1 + beta*logKap2
360 100 : dlnkap_dlnRho = alfa*dlnkap1_dlnRho + beta*dlnkap2_dlnRho
361 : dlnkap_dlnT = &
362 : alfa*dlnkap1_dlnT + beta*dlnkap2_dlnT + &
363 100 : d_alfa_dlnT*logKap1*ln10 + d_beta_dlnT*logKap2*ln10
364 100 : kap = exp10(logKap)
365 :
366 100 : frac_lowT = beta
367 100 : frac_highT = alfa
368 :
369 : if (dbg) then
370 : write(*,1) 'alfa', alfa
371 : write(*,1) 'kap1 (lowT)', kap1
372 : write(*,1) 'beta', beta
373 : write(*,1) 'kap2 (highT)', kap2
374 : write(*,1) 'kap', kap
375 : write(*,'(A)')
376 : end if
377 :
378 : if (dbg) write(*,1) 'Get_kap_blend_logT kap', kap
379 :
380 : end subroutine get_kap_Results_blend_T
381 :
382 :
383 5952 : subroutine get_kap_lowT_Results( &
384 : rq, &
385 : X, Z, XC, XN, XO, XNe, logRho_in, logT_in, &
386 : frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
387 :
388 : use kap_def
389 :
390 : use kap_eval_fixed
391 : use kap_eval_co
392 : use kapcn
393 : use kap_aesopus
394 :
395 : use utils_lib, only: is_bad
396 :
397 : ! INPUT
398 : type (Kap_General_Info), pointer :: rq
399 : real(dp), intent(in) :: X, Z, XC, XN, XO, XNe ! composition
400 : real(dp), intent(in) :: logRho_in ! density
401 : real(dp), intent(in) :: logT_in ! temperature
402 : ! free_e := total combined number per nucleon of free electrons and positrons
403 :
404 : ! OUTPUT
405 : real(dp), intent(out) :: frac_Type2
406 : real(dp), intent(out) :: kap ! opacity
407 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
408 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
409 : integer, intent(out) :: ierr ! 0 means AOK.
410 :
411 : real(dp) :: logRho, Rho, logT, T
412 : real(dp) :: frac_Type1, Zbase, dZ, frac_X, frac_dZ, &
413 5952 : dXC, dXO, fCO, fC, fN, fO, fNe, fHeavy, ZHeavy, dXsum
414 :
415 5952 : real(dp) :: logKap
416 :
417 : logical, parameter :: dbg = .false.
418 :
419 5952 : ierr = -1 ! should be set by each case, otherwise something is wrong
420 5952 : frac_Type2 = 0d0
421 :
422 5952 : logRho = logRho_in; Rho = exp10(logRho)
423 5952 : logT = logT_in; T = exp10(logT)
424 :
425 5952 : Zbase = rq% Zbase
426 :
427 5952 : select case (rq% kap_lowT_option)
428 :
429 : case (kap_lowT_kapCN)
430 : if (dbg) write(*,*) 'Calling kapCN for lowT'
431 :
432 0 : if (Zbase < 0d0) then
433 0 : write(*,*) 'must supply Zbase for kapCN opacities', Zbase
434 0 : call mesa_error(__FILE__,__LINE__)
435 0 : ierr = -1
436 0 : return
437 : end if
438 :
439 0 : fC = XC / (kapCN_ZC * Zbase)
440 0 : fN = XN / (kapCN_ZN * Zbase)
441 : call kapCN_get(Zbase, X, fC, fN, logRho, logT, &
442 2 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
443 :
444 : case (kap_lowT_AESOPUS)
445 : if (dbg) write(*,*) 'Calling AESOPUS for lowT'
446 :
447 2 : if (Zbase < 0d0) then
448 0 : write(*,*) 'must supply Zbase for AESOPUS opacities', Zbase
449 0 : ierr = -1
450 0 : return
451 : end if
452 :
453 : call AESOPUS_get(Zbase, X, XC, XN, XO, logRho, logT, &
454 2 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
455 :
456 :
457 : case default
458 : if (rq% kap_lowT_option == kap_lowT_test) then
459 : if (dbg) write(*,*) 'Calling test for lowT'
460 : else if (rq% kap_lowT_option == kap_lowT_Freedman11) then
461 : if (dbg) write(*,*) 'Calling Freedman for lowT'
462 : else
463 : if (dbg) write(*,*) 'Calling Ferg for lowT'
464 : end if
465 :
466 : call Get1_kap_fixed_metal_Results( &
467 : kap_lowT_z_tables(rq% kap_lowT_option)% ar, num_kap_lowT_Zs(rq% kap_lowT_option), rq, Z, X, Rho, logRho, T, logT, &
468 5950 : logKap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
469 5954 : kap = exp10(logKap)
470 : end select
471 :
472 : return
473 :
474 5952 : end subroutine get_kap_lowT_Results
475 :
476 :
477 80364 : subroutine get_kap_highT_Results(rq, &
478 : X, Z, XC_in, XN_in, XO_in, XNe_in, logRho_in, logT_in, &
479 : frac_Type2, kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
480 :
481 5952 : use kap_def
482 :
483 : use kap_eval_fixed
484 : use kap_eval_co
485 :
486 : use utils_lib, only: is_bad
487 :
488 : ! INPUT
489 : type (Kap_General_Info), pointer :: rq
490 : real(dp), intent(in) :: X, Z, XC_in, XN_in, XO_in, XNe_in ! composition
491 : real(dp), intent(in) :: logRho_in ! density
492 : real(dp), intent(in) :: logT_in ! temperature
493 :
494 : ! OUTPUT
495 : real(dp), intent(out) :: frac_Type2
496 : real(dp), intent(out) :: kap ! opacity
497 : real(dp), intent(out) :: dlnkap_dlnRho ! partial derivative at constant T
498 : real(dp), intent(out) :: dlnkap_dlnT ! partial derivative at constant Rho
499 : integer, intent(out) :: ierr ! 0 means AOK.
500 :
501 : real(dp) :: logRho, Rho, logT, T
502 80364 : real(dp) :: XC, XN, XO, XNe
503 80364 : real(dp) :: frac_Type1, Zbase, dZ, frac_X, frac_dZ, &
504 80364 : dXC, dXO, fC, fN, fO, fNe, fHeavy, ZHeavy, dXsum
505 :
506 80364 : real(dp) :: max_frac_Type2
507 :
508 80364 : real(dp) :: beta, kap_beta, logkap_beta, dlnkap_beta_dlnRho, dlnkap_beta_dlnT, &
509 80364 : alfa, kap_alfa, logkap_alfa, dlnkap_alfa_dlnRho, dlnkap_alfa_dlnT
510 :
511 : logical, parameter :: dbg = .false.
512 :
513 : include 'formats'
514 :
515 80364 : logRho = logRho_in; Rho = exp10(logRho)
516 80364 : logT = logT_in; T = exp10(logT)
517 :
518 80364 : Zbase = rq% Zbase
519 :
520 80364 : if (rq% use_Type2_opacities) then
521 80362 : if (Zbase < 0d0) then
522 0 : write(*,1) 'must supply Zbase for Type2 opacities', Zbase
523 0 : call mesa_error(__FILE__,__LINE__)
524 0 : ierr = -1
525 0 : return
526 : end if
527 : end if
528 :
529 80364 : if (rq% use_Type2_opacities) then
530 80362 : max_frac_Type2 = 1d0
531 80362 : XC = XC_in
532 80362 : XN = XN_in
533 80362 : XO = XO_in
534 80362 : XNe = XNe_in
535 : else
536 2 : max_frac_Type2 = 0d0
537 2 : XC = 0d0
538 2 : XN = 0d0
539 2 : XO = 0d0
540 2 : XNe = 0d0
541 : end if
542 :
543 :
544 80364 : dXC = 0d0
545 80364 : dxO = 0d0
546 80364 : frac_Type2 = max_frac_Type2
547 :
548 :
549 : ! Type2 opacities use a base metallicity and enhancements to C and O
550 : ! from the 3 definitions
551 : ! Z = Zbase + dXC + dXO ! total metals
552 : ! XC = dXC + base_fC*Zbase ! total mass fraction of carbon
553 : ! XO = dXO + base_fO*Zbase ! total mass fraction of oxygen
554 : ! we get expressions for the 3 parameters, Zbase, dXC, and dXO
555 : ! using the base values for fractions of C and O
556 :
557 80364 : if (frac_Type2 > 0d0 .and. &
558 : associated(kap_co_z_tables(rq% kap_CO_option)% ar)) then
559 :
560 80362 : fC = kap_co_z_tables(rq% kap_CO_option)% ar(1)% Zfrac_C
561 80362 : fN = kap_co_z_tables(rq% kap_CO_option)% ar(1)% Zfrac_N
562 80362 : fO = kap_co_z_tables(rq% kap_CO_option)% ar(1)% Zfrac_O
563 80362 : fNe = kap_co_z_tables(rq% kap_CO_option)% ar(1)% Zfrac_Ne
564 :
565 80362 : dXC = max(0d0, xC - fC*Zbase)
566 80362 : dXO = max(0d0, xO - fO*Zbase)
567 :
568 80362 : if (X >= rq% kap_Type2_full_off_X) then
569 80360 : frac_X = 0d0
570 2 : else if (X <= rq% kap_Type2_full_on_X) then
571 2 : frac_X = 1d0
572 : else ! blend
573 : frac_X = (rq% kap_Type2_full_off_X - X) / &
574 0 : (rq% kap_Type2_full_off_X - rq% kap_Type2_full_on_X)
575 : end if
576 :
577 80362 : dZ = Z - Zbase
578 80362 : if (dZ <= rq% kap_Type2_full_off_dZ) then
579 80360 : frac_dZ = 0d0
580 2 : else if (dZ >= rq% kap_Type2_full_on_dZ) then
581 2 : frac_dZ = 1d0
582 : else ! blend
583 : frac_dZ = (rq% kap_Type2_full_off_dZ - dZ) / &
584 0 : (rq% kap_Type2_full_off_dZ - rq% kap_Type2_full_on_dZ)
585 : end if
586 :
587 80362 : frac_Type2 = frac_Type2*frac_X*frac_dZ
588 80362 : if (frac_Type2 == 0d0) then
589 80360 : dXC = 0d0
590 80360 : dXO = 0d0
591 : end if
592 :
593 80362 : if (is_bad(dXC) .or. is_bad(dXO)) then
594 0 : write(*,1) 'X', X
595 0 : write(*,1) 'Z', Z
596 0 : write(*,1) 'Zbase', Zbase
597 0 : write(*,1) 'XC', XC
598 0 : write(*,1) 'XN', XN
599 0 : write(*,1) 'XO', XO
600 0 : write(*,1) 'XNe', XNe
601 0 : write(*,1) 'logRho', logRho
602 0 : write(*,1) 'logT', logT
603 0 : write(*,1) 'dXC', dXC
604 0 : write(*,1) 'dXO', dXO
605 0 : write(*,1) 'frac_X', frac_X
606 0 : write(*,1) 'frac_dZ', frac_dZ
607 0 : write(*,1) 'frac_Type2', frac_Type2
608 0 : ierr = -1
609 0 : return
610 : end if
611 :
612 : end if
613 :
614 80364 : frac_Type1 = 1d0 - frac_Type2
615 : if (dbg) then
616 : write(*,1) 'max_frac_Type2', max_frac_Type2
617 : write(*,1) 'frac_Type2', frac_Type2
618 : write(*,1) 'frac_Type1', frac_Type1
619 : write(*,'(A)')
620 : write(*,1) 'X', X
621 : write(*,1) 'dXC', dXC
622 : write(*,1) 'dXO', dXO
623 : write(*,1) 'Zbase', Zbase
624 : end if
625 :
626 : ! do blend
627 80364 : beta = frac_Type1
628 80364 : alfa = frac_Type2
629 :
630 80364 : if (beta > 0) then ! get value from fixed metal tables
631 80362 : if (rq% use_Type2_opacities .and. rq% use_Zbase_for_Type1) then
632 : ! when Z > Zbase, use Zbase instead of Z.
633 : call Get1_kap_fixed_metal_Results( &
634 : kap_z_tables(rq% kap_option)% ar, num_kap_Zs(rq% kap_option), rq, min(Z,Zbase), X, Rho, logRho, T, logT, &
635 80360 : logkap_beta, dlnkap_beta_dlnRho, dlnkap_beta_dlnT, ierr)
636 : else
637 : call Get1_kap_fixed_metal_Results( &
638 : kap_z_tables(rq% kap_option)% ar, num_kap_Zs(rq% kap_option), rq, Z, X, Rho, logRho, T, logT, &
639 2 : logkap_beta, dlnkap_beta_dlnRho, dlnkap_beta_dlnT, ierr)
640 : end if
641 80362 : kap_beta = exp10(logkap_beta)
642 : if (dbg) then
643 : write(*,1) 'Get_kap_fixed_metal_Results'
644 : write(*,1) 'Z', Z
645 : write(*,1) 'X', X
646 : write(*,1) 'Rho', Rho
647 : write(*,1) 'logRho', logRho
648 : write(*,1) 'T', T
649 : write(*,1) 'logT', logT
650 : write(*,1) 'kap_beta', kap_beta
651 : write(*,1) 'beta', beta
652 : write(*,'(A)')
653 : end if
654 : else
655 2 : kap_beta = 0d0; dlnkap_beta_dlnRho = 0d0; dlnkap_beta_dlnT = 0d0
656 : end if
657 :
658 80364 : if (alfa > 0d0) then ! get value from C/O enhanced tables
659 : call Get1_kap_CO_Results( &
660 : rq, Zbase, X, dXC, dXO, Rho, logRho, T, logT, &
661 2 : logkap_alfa, dlnkap_alfa_dlnRho, dlnkap_alfa_dlnT, ierr)
662 2 : kap_alfa = exp10(logkap_alfa)
663 : if (dbg) then
664 : write(*,1) 'Get_kap_CO_Results'
665 : write(*,1) 'Z', Z
666 : write(*,1) 'X', X
667 : write(*,1) 'Rho', Rho
668 : write(*,1) 'logRho', logRho
669 : write(*,1) 'T', T
670 : write(*,1) 'logT', logT
671 : write(*,1) 'kap_alfa', kap_alfa
672 : write(*,1) 'alfa', alfa
673 : write(*,'(A)')
674 : end if
675 : else
676 80362 : kap_alfa = 0d0; dlnkap_alfa_dlnRho = 0d0; dlnkap_alfa_dlnT = 0d0
677 : end if
678 :
679 80364 : kap = alfa*kap_alfa + beta*kap_beta
680 : if (dbg) then
681 : write(*,1) 'alfa', alfa
682 : write(*,1) 'kap_alfa (Type2)', kap_alfa
683 : write(*,1) 'beta', beta
684 : write(*,1) 'kap_beta (Type1)', kap_beta
685 : write(*,1) 'kap', kap
686 : end if
687 :
688 80364 : if (kap < 1d-30) then
689 0 : kap = 1d-30
690 0 : dlnkap_dlnRho = 0d0
691 0 : dlnkap_dlnT = 0d0
692 : else
693 80364 : dlnkap_dlnRho = (alfa*kap_alfa*dlnkap_alfa_dlnRho + beta*kap_beta*dlnkap_beta_dlnRho) / kap
694 80364 : dlnkap_dlnT = (alfa*kap_alfa*dlnkap_alfa_dlnT + beta*kap_beta*dlnkap_beta_dlnT) / kap
695 : end if
696 :
697 80364 : end subroutine get_kap_highT_Results
698 :
699 :
700 86216 : subroutine combine_rad_with_conduction( &
701 : rq, logRho, logT, zbar, &
702 : kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT, &
703 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
704 :
705 80364 : use condint, only: do_electron_conduction
706 : type (Kap_General_Info), pointer :: rq
707 : real(dp), intent(in) :: logRho, logT, zbar
708 : real(dp), intent(inout) :: kap_rad, dlnkap_rad_dlnRho, dlnkap_rad_dlnT
709 : real(dp), intent(out) :: kap, dlnkap_dlnRho, dlnkap_dlnT
710 : integer, intent(out) :: ierr ! 0 means AOK.
711 :
712 86216 : real(dp) :: kap_ec, dlnkap_ec_dlnRho, dlnkap_ec_dlnT
713 : logical, parameter :: dbg = .false.
714 :
715 : include 'formats'
716 :
717 86216 : ierr = 0
718 :
719 86216 : if (.not. rq% include_electron_conduction) then
720 0 : kap = kap_rad
721 0 : dlnkap_dlnRho = dlnkap_rad_dlnRho
722 0 : dlnkap_dlnT = dlnkap_rad_dlnT
723 0 : return
724 : end if
725 :
726 86216 : if (rq% use_other_elect_cond_opacity) then
727 : call rq% other_elect_cond_opacity( &
728 : rq% handle, zbar, logRho, logT, &
729 0 : kap_ec, dlnkap_ec_dlnRho, dlnkap_ec_dlnT, ierr)
730 : else
731 : call do_electron_conduction( &
732 : rq, zbar, logRho, logT, &
733 86216 : kap_ec, dlnkap_ec_dlnRho, dlnkap_ec_dlnT, ierr)
734 : end if
735 86216 : if (ierr /= 0) return
736 :
737 86216 : if (is_bad(kap_ec)) then
738 0 : write(*,*) 'kap_ec', kap_ec
739 0 : call mesa_error(__FILE__,__LINE__,'combine_rad_with_conduction')
740 : end if
741 : if (dbg) write(*,1) 'kap_ec', kap_ec
742 : if (dbg) write(*,1) 'dlnkap_ec_dlnRho', dlnkap_ec_dlnRho
743 : if (dbg) write(*,1) 'dlnkap_ec_dlnT', dlnkap_ec_dlnT
744 : if (dbg) write(*,*)
745 :
746 86216 : kap = 1d0 / (1d0/kap_rad + 1d0/kap_ec)
747 : if (dbg) write(*,1) 'kap_rad', kap_rad
748 : if (dbg) write(*,1) 'kap', kap
749 : if (dbg) write(*,1) 'log10(kap)', log10(kap)
750 :
751 86216 : if (is_bad(kap)) then
752 0 : ierr = -1; return
753 : write(*,1) 'kap', kap
754 : call mesa_error(__FILE__,__LINE__,'Get_kap_Results')
755 : end if
756 :
757 86216 : dlnkap_dlnRho = (kap/kap_rad) * dlnkap_rad_dlnRho + (kap/kap_ec) * dlnkap_ec_dlnRho
758 :
759 86216 : if (is_bad(dlnkap_dlnRho)) then
760 0 : ierr = -1; return
761 : write(*,1) 'dlnkap_dlnRho', dlnkap_dlnRho
762 : write(*,1) 'kap', kap
763 : write(*,1) 'dkap_dlnRho', kap * dlnkap_dlnRho
764 : write(*,1) 'dkap_ec_dlnRho', kap_ec * dlnkap_ec_dlnRho
765 : write(*,1) 'dkap_rad_dlnRho', kap_rad * dlnkap_rad_dlnRho
766 : write(*,1) 'kap_rad', kap_rad
767 : write(*,1) 'kap_ec', kap_ec
768 : call mesa_error(__FILE__,__LINE__,'combine_rad_with_conduction')
769 : end if
770 :
771 86216 : dlnkap_dlnT = (kap/kap_rad) * dlnkap_rad_dlnT + (kap/kap_ec) * dlnkap_ec_dlnT
772 :
773 86216 : if (is_bad(dlnkap_dlnT)) then
774 0 : ierr = -1; return
775 : write(*,1) 'dlnkap_dlnT', dlnkap_dlnT
776 : write(*,1) 'kap', kap
777 : write(*,1) 'dkap_dlnT', kap * dlnkap_dlnT
778 : write(*,1) 'dkap_ec_dlnT', kap_ec * dlnkap_ec_dlnT
779 : write(*,1) 'dkap_rad_dlnT', kap_rad * dlnkap_rad_dlnT
780 : write(*,1) 'kap_rad', kap_rad
781 : write(*,1) 'kap_ec', kap_ec
782 : call mesa_error(__FILE__,__LINE__,'combine_rad_with_conduction')
783 : end if
784 :
785 : if (dbg) write(*,1) 'dlnkap_dlnRho', dlnkap_dlnRho
786 : if (dbg) write(*,1) 'dlnkap_dlnT', dlnkap_dlnT
787 : if (dbg) call mesa_error(__FILE__,__LINE__,'combine_rad_with_conduction')
788 :
789 :
790 86216 : end subroutine combine_rad_with_conduction
791 :
792 :
793 0 : subroutine Compton_Opacity(rq, &
794 : Rho_in, T_in, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT, &
795 : eta_in, d_eta_dlnRho, d_eta_dlnT, &
796 : kap, dlnkap_dlnRho, dlnkap_dlnT, ierr)
797 86216 : use eos_lib
798 : use eos_def
799 : use auto_diff
800 :
801 : type (Kap_General_Info), pointer :: rq
802 :
803 : ! evaluates the poutanen 2017, apj 835, 119 fitting formula for the compton opacity.
804 : ! coefficients from table 1 between 2 and 300 kev
805 :
806 : real(dp), intent(in) :: Rho_in, T_in, lnfree_e, d_lnfree_e_dlnRho, d_lnfree_e_dlnT
807 : real(dp), intent(in) :: eta_in, d_eta_dlnRho, d_eta_dlnT
808 : real(dp), intent(out) :: kap, dlnkap_dlnRho, dlnkap_dlnT
809 : integer, intent(out) :: ierr
810 :
811 : type(auto_diff_real_2var_order1) :: T, rho, free_e, eta, kap_auto
812 : type(auto_diff_real_2var_order1) :: zeta, f1, f2, f3, alpha, tbr, theta, tkev, mfp
813 :
814 : real(dp) :: free_e_val
815 :
816 : ! poutanen table 1
817 : real(dp), parameter :: &
818 : t0 = 43.3d0, &
819 : alpha0 = 0.885d0, &
820 : c01 = 0.682d0, &
821 : c02 = -0.0454d0, &
822 : c11 = 0.240d0, &
823 : c12 = 0.0043d0, &
824 : c21 = 0.050d0, &
825 : c22 = -0.0067d0, &
826 : c31 = -0.037d0, &
827 : c32 = 0.0031d0
828 :
829 : include 'formats'
830 :
831 0 : ierr = 0
832 :
833 : ! set up auto diff
834 : ! var1: Rho
835 : ! var2: T
836 :
837 0 : Rho = Rho_in
838 : Rho% d1val1 = 1d0
839 : Rho% d1val2 = 0d0
840 :
841 0 : T = T_in
842 0 : T% d1val1 = 0d0
843 0 : T% d1val2 = 1d0
844 :
845 0 : free_e_val = exp(lnfree_e)
846 0 : free_e = free_e_val
847 0 : free_e % d1val1 = free_e_val*d_lnfree_e_dlnRho/Rho_in
848 0 : free_e % d1val2 = free_e_val*d_lnfree_e_dlnT/T_in
849 :
850 0 : eta = eta_in
851 0 : eta % d1val1 = d_eta_dlnRho/Rho_in
852 0 : eta % d1val2 = d_eta_dlnT/T_in
853 :
854 0 : theta = T * kerg / (me * clight * clight)
855 0 : tkev = T * kev * 1.0d-3
856 0 : zeta = exp(c01 * eta + c02 * eta*eta)
857 0 : f1 = 1.0d0 + c11 * zeta + c12 * zeta * zeta
858 0 : f2 = 1.0d0 + c21 * zeta + c22 * zeta * zeta
859 0 : f3 = 1.0d0 + c31 * zeta + c32 * zeta * zeta
860 0 : alpha = alpha0 * f3
861 0 : tbr = t0 * f2
862 0 : mfp = f1 * (1.0d0 + pow(tkev/tbr,alpha))
863 :
864 : ! equation 31
865 0 : kap_auto = sige*(free_e/amu)/mfp
866 :
867 : ! unpack auto_diff
868 0 : kap = kap_auto% val
869 0 : dlnkap_dlnRho = Rho% val * kap_auto% d1val1 / kap
870 0 : dlnkap_dlnT = T% val * kap_auto% d1val2 / kap
871 :
872 0 : end subroutine Compton_Opacity
873 :
874 : end module kap_eval
|
