Line data Source code
1 : module test_eos_support
2 :
3 : use eos_def
4 : use eos_lib
5 : use chem_lib
6 : use chem_def
7 : use const_def
8 : use eos_support
9 : use math_lib
10 :
11 : implicit none
12 :
13 : contains
14 :
15 2 : subroutine test1_eosPT_for_ck(quietly)
16 : logical, intent(in) :: quietly
17 : real(dp) :: Z, X, logPgas, logT, logRho, logP
18 : ! pick args for test so that use both HELM and tables to get results.
19 2 : Z = 0.02d0
20 2 : X = 0.7d0
21 2 : logT = 3.0d0
22 2 : logPgas = 4.2d0
23 2 : call test1_eosPT(Z, X, logPgas, logT, .false., quietly, logRho, logP)
24 2 : end subroutine test1_eosPT_for_ck
25 :
26 0 : subroutine test_eosPT(which)
27 : integer, intent(in) :: which
28 : real(dp) :: Z, X, logPgas, logT, logRho, logP
29 0 : Z = 0.02d0
30 0 : X = 0.6d0
31 0 : logT = 4d0
32 0 : logPgas = 5d0
33 0 : call test1_eosPT(Z, X, logPgas, logT, .false., .false., logRho, logP)
34 0 : end subroutine test_eosPT
35 :
36 2 : subroutine test1_eosPT(Z, X, logPgas, logT, do_compare, quietly, logRho, logP)
37 : logical, intent(in) :: quietly
38 : real(dp), intent(in) :: Z, X, logPgas, logT
39 : real(dp), intent(out) :: logRho, logP
40 : logical, intent(in) :: do_compare
41 : real(dp) :: &
42 4 : P, Pgas, Prad, T, Rho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
43 108 : res(num_eos_basic_results), d_dlnd(num_eos_basic_results), &
44 54 : d_dlnT(num_eos_basic_results), &
45 108 : res2(num_eos_basic_results), d_dlnd2(num_eos_basic_results), &
46 44 : d_dxa2(num_eos_d_dxa_results, species), &
47 54 : d_dlnT2(num_eos_basic_results)
48 : integer:: ierr, i
49 : character(len=eos_name_length) :: names(num_eos_basic_results)
50 :
51 : include 'formats'
52 :
53 2 : ierr = 0
54 :
55 2 : call Init_Composition(X, Z, 0d0, 0d0) ! sets abar and zbar
56 :
57 : !if (.false.) then ! TESTING
58 : ! Z = 0d0
59 : ! X = 0.72d0
60 : ! call Init_Composition(X, Z, 0d0, 0d0) ! sets abar and zbar
61 : ! abar = 1.2966413082679851D+00
62 : ! zbar = 1.1021433867453336D+00
63 : ! logPgas = 4.8066181993619859D+00
64 : ! logT = 3.7569035961895620D+00
65 : !end if
66 :
67 2 : T = exp10(logT)
68 2 : Pgas = exp10(logPgas)
69 :
70 2 : if (.not. quietly) then
71 1 : write (*, *) 'test1_eosPT'
72 1 : write (*, 1) 'logT', logT
73 1 : write (*, 1) 'logPgas', logPgas
74 1 : write (*, 1) 'T', T
75 1 : write (*, 1) 'Pgas', Pgas
76 1 : write (*, 1) 'Z', Z
77 1 : write (*, 1) 'X', X
78 1 : write (*, 1) 'abar', abar
79 1 : write (*, 1) 'zbar', zbar
80 1 : write (*, '(A)')
81 : end if
82 :
83 : call eosPT_get( &
84 : handle, &
85 : species, chem_id, net_iso, xa, &
86 : Pgas, logPgas, T, logT, &
87 : Rho, logRho, dlnRho_dlnPgas_c_T, dlnRho_dlnT_c_Pgas, &
88 2 : res, d_dlnd, d_dlnT, d_dxa, ierr)
89 2 : if (ierr /= 0) then
90 0 : write (*, *) 'ierr in eosPT_get for test1_eosPT'
91 0 : call mesa_error(__FILE__, __LINE__)
92 : end if
93 :
94 2 : Prad = crad*T*T*T*T/3
95 2 : P = Pgas + Prad
96 2 : logP = log10(P)
97 :
98 3 : if (quietly) return
99 :
100 1 : write (*, '(A)')
101 1 : write (*, 1) 'rho', rho
102 1 : write (*, 1) 'logRho', logRho
103 1 : write (*, 1) 'logT', logT
104 1 : write (*, 1) 'logPgas', logPgas
105 1 : write (*, '(A)')
106 :
107 27 : names = eosDT_result_names
108 :
109 1 : if (.not. do_compare) then ! simple form of output
110 1 : write (*, 1) 'dlnRho_dlnPgas_c_T', dlnRho_dlnPgas_c_T
111 1 : write (*, 1) 'dlnRho_dlnT_c_Pgas', dlnRho_dlnT_c_Pgas
112 1 : write (*, '(A)')
113 1 : write (*, 1) 'P', P
114 1 : write (*, 1) 'E', exp(res(i_lnE))
115 1 : write (*, 1) 'S', exp(res(i_lnS))
116 1 : write (*, '(A)')
117 24 : do i = 4, num_eos_basic_results
118 24 : write (*, 1) trim(names(i)), res(i)
119 : end do
120 1 : write (*, '(A)')
121 1 : return
122 : end if
123 :
124 : call eosDT_get( &
125 : handle, &
126 : species, chem_id, net_iso, xa, &
127 : rho, logRho, T, logT, &
128 0 : res2, d_dlnd2, d_dlnT2, d_dxa2, ierr)
129 0 : if (ierr /= 0) then
130 0 : write (*, *) 'ierr in eosDT_get for test1_eosPT'
131 0 : call mesa_error(__FILE__, __LINE__)
132 : end if
133 :
134 0 : write (*, '(A)')
135 :
136 0 : write (*, 1) 'dlnRho_dlnPgas_c_T', dlnRho_dlnPgas_c_T
137 0 : write (*, 1) 'dlnRho_dlnT_c_Pgas', dlnRho_dlnT_c_Pgas
138 0 : do i = 1, num_eos_basic_results
139 0 : write (*, 1) trim(names(i)), res(i), res2(i), &
140 0 : (res(i) - res2(i))/max(1d0, abs(res(i)), abs(res2(i)))
141 : end do
142 0 : write (*, '(A)')
143 :
144 0 : do i = 1, num_eos_basic_results
145 0 : write (*, 1) 'd_dlnd '//trim(names(i)), d_dlnd(i), d_dlnd2(i), &
146 0 : (d_dlnd(i) - d_dlnd2(i))/max(1d0, abs(d_dlnd(i)), abs(d_dlnd2(i)))
147 : end do
148 0 : write (*, '(A)')
149 :
150 0 : do i = 1, num_eos_basic_results
151 0 : write (*, 1) 'd_dlnT '//trim(names(i)), d_dlnT(i), &
152 0 : d_dlnT2(i), (d_dlnT(i) - d_dlnT2(i))/max(1d0, abs(d_dlnT(i)), abs(d_dlnT2(i)))
153 : end do
154 0 : write (*, '(A)')
155 :
156 : end subroutine test1_eosPT
157 :
158 2 : subroutine Do_One(quietly)
159 : logical, intent(in) :: quietly
160 2 : real(dp) :: T, rho
161 54 : real(dp), dimension(num_eos_basic_results) :: res
162 : real(dp) :: dXC
163 :
164 : if (.true.) then
165 : ! pure Helium
166 2 : X = 0.00d0
167 2 : Zinit = 0.00d0
168 2 : dXO = 0.00d0
169 2 : dXC = 0.00d0
170 2 : call doit('pure Helium')
171 : ! pure Hydrogen
172 2 : X = 1.00d0
173 2 : Zinit = 0.00d0
174 2 : dXO = 0.00d0
175 2 : dXC = 0.00d0
176 2 : call doit('pure Hydrogen')
177 : ! mixed Z with H&He 3:1 ratio
178 2 : Zinit = 0.03d0
179 2 : X = 0.75d0*(1 - Zinit)
180 2 : dXO = 0.00d0
181 2 : dXC = 0.00d0
182 2 : call doit('mixed Z with H&He 3:1 ratio')
183 : ! solar
184 2 : X = 0.70d0
185 2 : Zinit = 0.02d0
186 2 : dXO = 0.00d0
187 2 : dXC = 0.00d0
188 2 : call doit('solar')
189 : end if
190 :
191 2 : if (.true.) then ! do get_Rho and get_T
192 2 : X = 0.70d+00
193 2 : Zinit = 0.02d0
194 2 : dXO = 0.00d0
195 2 : dXC = 0.00d0
196 2 : T = exp10(4.8d0)
197 2 : rho = 1d-7
198 2 : Z = Zinit + dXC + dXO
199 2 : Y = 1 - (X + Z)
200 2 : call Init_Composition(X, Zinit, dXC, dXO)
201 2 : res(i_lnS) = log(2.9680645120000000d+09)
202 2 : call test_get_Rho_T
203 2 : if (.not. quietly) write (*, *)
204 : end if
205 :
206 : contains
207 :
208 40 : subroutine doit(str)
209 : character(len=*), intent(in) :: str
210 :
211 : if (.false.) then
212 : T = 2d8; rho = 100d0
213 : call Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res) ! scvh
214 : stop
215 : end if
216 :
217 8 : if (.not. quietly) write (*, *) trim(str)
218 :
219 8 : T = 1d6; rho = 1d-2
220 8 : call Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res) ! opal
221 8 : T = 1d4; rho = 1d-1
222 8 : call Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res) ! scvh
223 8 : T = 1d5; rho = 1d-1
224 8 : call Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res) ! opal-scvh overlap
225 8 : T = 2d8; rho = 1d2
226 8 : call Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res) ! helm
227 :
228 8 : if (.not. quietly) write (*, *)
229 :
230 8 : end subroutine doit
231 :
232 2 : subroutine test_get_Rho_T ! using most recent values from subroutine Do_One_TRho
233 2 : real(dp) :: tol, othertol, &
234 2 : result, result_log10, log10_T, log10_rho, lnS, Prad, Pgas, logP, &
235 2 : logRho_guess, logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
236 2 : logT_guess, logT_bnd1, logT_bnd2
237 : integer :: max_iter, eos_calls, ierr
238 : real(dp), dimension(num_eos_basic_results) :: &
239 106 : d_dlnd, d_dlnT
240 : real(dp), dimension(num_eos_d_dxa_results, species) :: &
241 44 : d_dxa
242 :
243 2 : if (.not. quietly) write (*, *)
244 :
245 2 : log10_rho = log10(rho)
246 2 : log10_T = log10(T)
247 2 : lnS = res(i_lnS)
248 :
249 2 : ierr = 0
250 2 : max_iter = 100
251 2 : tol = 1d-5
252 2 : othertol = 1d-12
253 :
254 : 1 format(a30, 1pe24.16)
255 :
256 2 : if (.not. quietly) then
257 1 : write (*, '(A)')
258 1 : write (*, 1) ' tolerance', tol
259 : end if
260 2 : if (.not. quietly) write (*, *)
261 2 : result = rho*0.5d0 ! initial guess
262 2 : result_log10 = log10(result)
263 2 : res = 0
264 2 : logRho_guess = result_log10
265 2 : logRho_bnd1 = arg_not_provided
266 2 : logRho_bnd2 = arg_not_provided
267 2 : other_at_bnd1 = arg_not_provided
268 2 : other_at_bnd2 = arg_not_provided
269 : call eosDT_get_Rho( &
270 : handle, &
271 : species, chem_id, net_iso, xa, &
272 : log10_T, i_lnS, lnS, &
273 : tol, othertol, max_iter, logRho_guess, &
274 : logRho_bnd1, logRho_bnd2, other_at_bnd1, other_at_bnd2, &
275 : result_log10, res, d_dlnd, d_dlnT, &
276 2 : d_dxa, eos_calls, ierr)
277 2 : result = exp10(result_log10)
278 2 : if (ierr /= 0) then
279 0 : write (*, *) 'ierr in test_get_Rho_T'
280 0 : call mesa_error(__FILE__, __LINE__)
281 : end if
282 2 : if (.not. quietly) then
283 1 : write (*, 1) 'actual logRho', log10_rho
284 1 : write (*, 1) ' guess logRho', logRho_guess
285 1 : write (*, 1) ' found logRho', result_log10
286 1 : write (*, 1) ' wanted logS', lnS/ln10
287 1 : write (*, 1) ' got logS', res(i_lnS)/ln10
288 1 : write (*, '(A)')
289 : end if
290 :
291 2 : result = T*2d0 ! initial guess
292 2 : result_log10 = log10(result)
293 2 : res = 0
294 2 : logT_guess = result_log10
295 2 : logT_bnd1 = 3
296 2 : logT_bnd2 = 9
297 : call eosDT_get_T( &
298 : handle, &
299 : species, chem_id, net_iso, xa, &
300 : log10_rho, i_lnS, lnS, &
301 : tol, othertol, max_iter, logT_guess, &
302 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
303 : result_log10, res, d_dlnd, d_dlnT, &
304 2 : d_dxa, eos_calls, ierr)
305 2 : result = exp10(result_log10)
306 2 : if (ierr /= 0) then
307 0 : write (*, *) 'ierr in test_get_Rho_T'
308 0 : call mesa_error(__FILE__, __LINE__)
309 : end if
310 2 : if (.not. quietly) then
311 1 : write (*, '(A)')
312 1 : write (*, 1) 'actual logT', log10_T
313 1 : write (*, 1) ' guess logT', logT_guess
314 1 : write (*, 1) ' found logT', result_log10
315 1 : write (*, 1) ' wanted logS', lnS/ln10
316 1 : write (*, 1) ' got logS', res(i_lnS)/ln10
317 1 : write (*, '(A)')
318 : end if
319 :
320 2 : T = exp10(log10_T)
321 2 : Prad = crad*T*T*T*T/3
322 2 : Pgas = exp(res(i_lnPgas))
323 2 : logP = log10(Prad + Pgas)
324 2 : result = T*2d0 ! initial guess
325 2 : result_log10 = log10(result)
326 2 : res = 0
327 2 : logT_guess = result_log10
328 : logT_bnd1 = 3
329 : logT_bnd2 = 9
330 : call eosDT_get_T( &
331 : handle, &
332 : species, chem_id, net_iso, xa, &
333 : log10_rho, i_logPtot, logP, &
334 : tol, othertol, max_iter, logT_guess, &
335 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
336 : result_log10, res, d_dlnd, d_dlnT, &
337 2 : d_dxa, eos_calls, ierr)
338 2 : result = exp10(result_log10)
339 2 : if (ierr /= 0) then
340 0 : write (*, *) 'ierr in test_get_Rho_T (eosDT_get_T_given_Ptotal)'
341 0 : call mesa_error(__FILE__, __LINE__)
342 : end if
343 2 : if (.not. quietly) then
344 1 : write (*, '(A)')
345 1 : write (*, 1) 'actual logT', log10_T
346 1 : write (*, 1) ' guess logT', logT_guess
347 1 : write (*, 1) ' found logT', result_log10
348 1 : write (*, 1) ' wanted logP', logP
349 1 : T = exp10(result_log10)
350 1 : Prad = crad*T*T*T*T/3
351 1 : Pgas = exp(res(i_lnPgas))
352 1 : logP = log10(Prad + Pgas)
353 1 : write (*, 1) ' got logP', logP
354 1 : write (*, '(A)')
355 : end if
356 :
357 2 : end subroutine test_get_Rho_T
358 :
359 : end subroutine Do_One
360 :
361 0 : subroutine test1_eosPT_get_T
362 :
363 : real(dp) :: &
364 0 : energy, abar, zbar, X, Z, logPgas, logT_tol, other_tol, other, &
365 0 : logT_guess, logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
366 0 : logT_result, &
367 0 : res(num_eos_basic_results), d_dlnd(num_eos_basic_results), &
368 0 : d_dxa(num_eos_d_dxa_results, species), &
369 0 : d_dlnT(num_eos_basic_results), &
370 0 : Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas
371 : integer:: ierr, which_other, eos_calls, max_iter
372 :
373 : 1 format(a40, 1pe26.16)
374 :
375 0 : call Setup_eos
376 :
377 : ierr = 0
378 :
379 0 : write (*, *) 'test1_eosPT_get_T'
380 :
381 0 : Z = 0.02d0
382 0 : X = 0.70d0
383 0 : abar = 1.2966353559153956d0
384 0 : zbar = 1.1021400447995373d0
385 0 : logPgas = 15d0
386 0 : energy = exp(3.5034294596213336d+01)
387 0 : logT_tol = 1d-6
388 0 : other_tol = ln10*1d-6
389 0 : logT_guess = 7d0
390 0 : logT_bnd1 = arg_not_provided
391 0 : logT_bnd2 = arg_not_provided
392 0 : other_at_bnd1 = arg_not_provided
393 0 : other_at_bnd2 = arg_not_provided
394 :
395 0 : which_other = i_lnE
396 0 : other = log(energy)
397 0 : max_iter = 100
398 :
399 0 : write (*, 1) 'logPgas', logPgas
400 0 : write (*, 1) 'logT_guess', logT_guess
401 0 : write (*, 1) 'logT_bnd1', logT_bnd1
402 0 : write (*, 1) 'logT_bnd2', logT_bnd2
403 0 : write (*, 1) 'energy', energy
404 0 : write (*, 1) 'other', other
405 0 : write (*, 1) 'Z', Z
406 0 : write (*, 1) 'X', X
407 0 : write (*, 1) 'abar', abar
408 0 : write (*, 1) 'zbar', zbar
409 0 : write (*, 1) 'logT_tol', logT_tol
410 0 : write (*, 1) 'other_tol', other_tol
411 0 : write (*, '(A)')
412 :
413 : call eosPT_get_T( &
414 : handle, &
415 : species, chem_id, net_iso, xa, &
416 : logPgas, which_other, other, &
417 : logT_tol, other_tol, max_iter, logT_guess, &
418 : logT_bnd1, logT_bnd2, other_at_bnd1, other_at_bnd2, &
419 : logT_result, Rho, log10Rho, dlnRho_dlnPgas_const_T, dlnRho_dlnT_const_Pgas, &
420 : res, d_dlnd, d_dlnT, d_dxa, &
421 0 : eos_calls, ierr)
422 0 : if (ierr /= 0) then
423 0 : write (*, *) 'ierr in test1_eosPT_get_T'
424 0 : call mesa_error(__FILE__, __LINE__)
425 : end if
426 0 : write (*, '(A)')
427 0 : write (*, 1) 'guess logT', logT_guess
428 0 : write (*, 1) 'found logT', logT_result
429 0 : write (*, 1) 'wanted logE', other/ln10
430 0 : write (*, 1) 'got logE', res(i_lnE)/ln10
431 0 : write (*, '(A)')
432 0 : write (*, *) 'eos_calls', eos_calls
433 0 : write (*, '(A)')
434 :
435 0 : end subroutine test1_eosPT_get_T
436 :
437 0 : subroutine test_components
438 : real(dp) :: &
439 : X_test, Z, XC_test, XO_test, &
440 0 : logT, logRho, Pgas, energy, entropy
441 : real(dp), dimension(num_eos_basic_results) :: &
442 0 : res, d_dlnd, d_dlnT
443 : integer:: ierr
444 : include 'formats'
445 :
446 0 : write (*, *) 'test_components'
447 :
448 0 : call Setup_eos
449 :
450 0 : ierr = 0
451 :
452 0 : X_test = 0.12d0
453 0 : Z = 0.03d0
454 0 : XC_test = 0d0
455 0 : XO_test = 0d0
456 :
457 0 : call Init_Composition(X_test, Z, XC_test, XO_test)
458 :
459 0 : logT = 6.0d0
460 0 : logRho = 3.5d0
461 :
462 0 : write (*, 1) 'logT', logT
463 0 : write (*, 1) 'logRho', logRho
464 0 : write (*, '(A)')
465 0 : write (*, 1) 'xa h1', xa(h1)
466 0 : write (*, 1) 'xa he4', xa(he4)
467 0 : write (*, 1) 'xa c12', xa(c12)
468 0 : write (*, 1) 'xa n14', xa(n14)
469 0 : write (*, 1) 'xa o16', xa(o16)
470 0 : write (*, 1) 'xa ne20', xa(ne20)
471 0 : write (*, 1) 'xa mg24', xa(mg24)
472 0 : write (*, '(A)')
473 0 : write (*, 1) 'X', X
474 0 : write (*, 1) 'Y', 1d0 - (X + Z)
475 0 : write (*, 1) 'Z', Z
476 0 : write (*, '(A)')
477 0 : write (*, 1) 'abar', abar
478 0 : write (*, 1) 'zbar', zbar
479 0 : write (*, '(A)')
480 0 : write (*, '(a55,9a26)') ' ', 'Pgas', 'energy', 'entropy'
481 0 : call test1(1, 'opal_scvh')
482 0 : call test1(2, 'helm')
483 0 : call test1(4, 'pc')
484 0 : write (*, '(A)')
485 :
486 : contains
487 :
488 0 : subroutine test1(which_eos, str)
489 : integer, intent(in) :: which_eos
490 : character(len=*), intent(in) :: str
491 : include 'formats'
492 : call eosDT_get_component( &
493 : handle, which_eos, &
494 : species, chem_id, net_iso, xa, &
495 : exp10(logRho), logRho, exp10(logT), logT, &
496 0 : res, d_dlnd, d_dlnT, d_dxa, ierr)
497 0 : if (ierr /= 0) then
498 0 : write (*, 1) trim(str)//' no results'
499 : else
500 0 : write (*, 1) trim(str), Pgas, energy, entropy
501 : end if
502 0 : end subroutine test1
503 :
504 : end subroutine test_components
505 :
506 0 : subroutine test1_eosDT_get_T_given_egas
507 :
508 : real(dp) :: &
509 0 : X, Z, abar, zbar, logRho, egas_want, egas_tol, logT_tol, logT_guess, &
510 0 : logT_bnd1, logT_bnd2, egas_at_bnd1, egas_at_bnd2, logT_result, erad, egas, energy, &
511 0 : res(num_eos_basic_results), d_dlnd(num_eos_basic_results), &
512 0 : d_dxa(num_eos_d_dxa_results, species), &
513 0 : d_dlnT(num_eos_basic_results)
514 : integer:: ierr, eos_calls, max_iter
515 :
516 : 1 format(a40, 1pe26.16)
517 :
518 0 : call Setup_eos
519 :
520 : ierr = 0
521 :
522 0 : Z = 0.7D-02
523 0 : X = 7.3D-01
524 :
525 0 : call Init_Composition(X, Z, 0d0, 0d0) ! sets abar and zbar
526 :
527 0 : write (*, *) 'test1_eosDT_get_T_given_egas'
528 :
529 0 : abar = 1.2559567378472252D+00
530 0 : zbar = 1.0864043570945732D+00
531 0 : logRho = -9.4201625429594529D+00
532 :
533 0 : egas_want = 2.0596457989663662D+12
534 0 : egas_tol = egas_want*1d-11
535 0 : logT_tol = 1d-11
536 0 : logT_guess = 3.6962155439999007D+00
537 :
538 0 : logT_bnd1 = arg_not_provided
539 0 : logT_bnd2 = arg_not_provided
540 0 : egas_at_bnd1 = arg_not_provided
541 0 : egas_at_bnd2 = arg_not_provided
542 :
543 0 : max_iter = 100
544 :
545 0 : write (*, 1) 'logRho', logRho
546 0 : write (*, 1) 'logT_guess', logT_guess
547 0 : write (*, 1) 'egas_want', egas_want
548 0 : write (*, 1) 'Z', Z
549 0 : write (*, 1) 'X', X
550 0 : write (*, 1) 'abar', abar
551 0 : write (*, 1) 'zbar', zbar
552 0 : write (*, 1) 'logT_tol', logT_tol
553 0 : write (*, 1) 'egas_tol', egas_tol
554 0 : write (*, '(A)')
555 :
556 : call eosDT_get_T( &
557 : handle, &
558 : species, chem_id, net_iso, xa, &
559 : logRho, i_egas, egas_want, &
560 : logT_tol, egas_tol, max_iter, logT_guess, &
561 : logT_bnd1, logT_bnd2, egas_at_bnd1, egas_at_bnd2, &
562 : logT_result, res, d_dlnd, d_dlnT, d_dxa, &
563 0 : eos_calls, ierr)
564 0 : if (ierr /= 0) then
565 0 : write (*, *) 'ierr in eosDT_get_T_given_egas'
566 0 : call mesa_error(__FILE__, __LINE__)
567 : end if
568 0 : energy = exp(res(i_lnE))
569 0 : erad = crad*exp10(logT_result)**4/exp10(logRho)
570 0 : egas = energy - erad
571 0 : write (*, '(A)')
572 0 : write (*, 1) 'guess logT', logT_guess
573 0 : write (*, 1) 'found logT', logT_result
574 0 : write (*, 1) 'wanted egas', egas_want
575 0 : write (*, 1) 'got egas', egas
576 0 : write (*, 1) '(want - got)/got', (egas_want - egas)/egas
577 0 : write (*, '(A)')
578 0 : write (*, *) 'eos_calls', eos_calls
579 0 : write (*, '(A)')
580 :
581 0 : end subroutine test1_eosDT_get_T_given_egas
582 :
583 32 : subroutine Do_One_TRho(quietly, T, Rho, X, Zinit, dXC, dXO, Y, Z, res)
584 : logical, intent(in) :: quietly
585 : real(dp), intent(in) :: T, Rho, X, Zinit, dXC, dXO
586 : real(dp), intent(out) :: Y, Z
587 : real(dp), intent(out), dimension(num_eos_basic_results) :: res
588 :
589 : real(dp), dimension(num_eos_basic_results) :: &
590 1696 : d_dlnd, d_dlnT
591 : real(dp), dimension(num_eos_d_dxa_results, species) :: &
592 704 : d_dxa
593 : integer :: info, i
594 32 : real(dp) :: Prad, Pgas, P
595 :
596 : 101 format(a30, 4x, 1pe24.16)
597 : 102 format(a30, 3x, 1pe24.16)
598 :
599 32 : Z = Zinit + dXC + dXO
600 32 : Y = 1 - (X + Z)
601 :
602 32 : call Init_Composition(X, Zinit, dXC, dXO)
603 :
604 32 : if (.not. quietly) then
605 16 : write (*, '(A)')
606 16 : write (*, '(A)')
607 16 : write (*, 102) 'X', X
608 16 : write (*, 102) 'Y', Y
609 16 : write (*, 102) 'Z', Z
610 16 : write (*, 102) 'abar', abar
611 16 : write (*, 102) 'zbar', zbar
612 16 : write (*, 102) 'logRho', log10(Rho)
613 16 : write (*, 102) 'logT', log10(T)
614 16 : write (*, 102) 'T6', T*1d-6
615 16 : write (*, '(A)')
616 : end if
617 :
618 : call eosDT_get( &
619 : handle, &
620 : species, chem_id, net_iso, xa, &
621 : Rho, arg_not_provided, T, arg_not_provided, &
622 32 : res, d_dlnd, d_dlnT, d_dxa, info)
623 32 : if (info /= 0) then
624 0 : write (*, *) 'info', info, 'Rho', Rho, 'T', T
625 0 : write (*, *) 'failed in Do_One_TRho'
626 0 : call mesa_error(__FILE__, __LINE__)
627 : end if
628 :
629 32 : if (.not. quietly) then
630 :
631 16 : write (*, *) 'eosDT_get'
632 16 : Prad = crad*T*T*T*T/3
633 16 : Pgas = exp(res(i_lnPgas))
634 16 : P = Pgas + Prad
635 16 : write (*, 101) 'P', P
636 16 : write (*, 101) 'E', exp(res(i_lnE))
637 16 : write (*, 101) 'S', exp(res(i_lnS))
638 112 : do i = 4, 9
639 112 : write (*, 101) trim(eos_names(i)), res(i)
640 : end do
641 16 : write (*, 101) trim(eos_names(i_gamma1)), res(i_gamma1)
642 16 : write (*, 101) trim(eos_names(i_gamma3)), res(i_gamma3)
643 16 : write (*, 101) trim(eos_names(i_eta)), res(i_eta)
644 :
645 : if (.false.) then ! debugging
646 : do i = 1, num_eos_basic_results
647 : write (*, 101) 'd_dlnd '//trim(eos_names(i)), d_dlnd(i)
648 : end do
649 : write (*, '(A)')
650 : do i = 1, num_eos_basic_results
651 : write (*, 101) 'd_dlnT '//trim(eos_names(i)), d_dlnT(i)
652 : end do
653 : write (*, '(A)')
654 : end if
655 :
656 : end if
657 :
658 32 : end subroutine Do_One_TRho
659 :
660 0 : subroutine test_dirac_integrals
661 0 : real(dp) :: T, eta, theta, fdph, fdmh, fdeta, fdtheta, theta_e
662 : 1 format(a40, 1pe26.16)
663 0 : eta = 1.46722890948893d0
664 0 : T = 11327678.5183021d0
665 0 : theta = (kerg*T)/(me*clight*clight)
666 : !zt 1.55623520289424d0
667 0 : theta_e = 0.929542529701454d0
668 0 : call eos_fermi_dirac_integral(-0.5d0, eta, theta, fdmh, fdeta, fdtheta)
669 0 : call eos_fermi_dirac_integral(0.5d0, eta, theta, fdph, fdeta, fdtheta)
670 0 : write (*, '(A)')
671 0 : write (*, *) 'test_dirac_integrals'
672 0 : write (*, 1) 'calculated theta_e', fdmh/fdph
673 0 : write (*, '(A)')
674 0 : stop
675 : end subroutine test_dirac_integrals
676 :
677 : end module test_eos_support
|
