Line data Source code
1 : module test_ionization_support
2 :
3 : use ionization_lib
4 : use chem_def
5 : use chem_lib, only: chem_init
6 : use const_lib, only: const_init
7 : use math_lib
8 : use utils_lib, only: mesa_error
9 :
10 : implicit none
11 :
12 : integer :: ierr
13 : logical, parameter :: use_cache = .true.
14 : character(len=256) :: my_mesa_dir, file_prefix, Z1_suffix, &
15 : in_dir, out_dir_ion, out_dir_eosDT, out_dir_eosPT
16 :
17 : contains
18 :
19 4 : subroutine do_test(quietly)
20 : logical, intent(in) :: quietly
21 :
22 2 : file_prefix = 'ion'
23 2 : Z1_suffix = '_CO_1'
24 2 : ierr = 0
25 :
26 2 : my_mesa_dir = '../..'
27 2 : call const_init(my_mesa_dir, ierr)
28 2 : if (ierr /= 0) then
29 0 : write (*, *) 'const_init failed'
30 0 : call mesa_error(__FILE__, __LINE__)
31 : end if
32 :
33 2 : call math_init()
34 :
35 2 : call chem_init('isotopes.data', ierr)
36 2 : if (ierr /= 0) then
37 0 : write (*, *) 'FATAL ERROR: failed in chem_init'
38 0 : call mesa_error(__FILE__, __LINE__)
39 : end if
40 :
41 2 : call ionization_init(file_prefix, Z1_suffix, '', use_cache, ierr)
42 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
43 :
44 : if (.false.) then
45 : in_dir = 'eos_macdonald'
46 : out_dir_ion = '' ! 'ionization_data'
47 : out_dir_eosDT = 'eosDT_data'
48 : out_dir_eosPT = 'eosPT_data'
49 : call create_ion_table_files( &
50 : in_dir, out_dir_ion, out_dir_eosDT, out_dir_eosPT)
51 : stop
52 : end if
53 :
54 : !call create_table_plot_files; stop
55 : !call Build_Plots; stop
56 :
57 2 : call test_fe56_in_he4(quietly)
58 2 : call do_test_eval_ionization(quietly)
59 :
60 2 : end subroutine do_test
61 :
62 2 : subroutine test_fe56_in_he4(quietly)
63 : logical, intent(in) :: quietly
64 : integer :: ierr
65 : 1 format(a40, f16.7)
66 2 : ierr = 0
67 2 : if (ierr /= 0) then
68 0 : write (*, *) 'eval_charge_of_Fe56_in_He4 failed'
69 0 : call mesa_error(__FILE__, __LINE__)
70 : end if
71 2 : if (.not. quietly) then
72 1 : write (*, *)
73 1 : write (*, *) 'test_fe56_in_he4'
74 : end if
75 2 : call do1_fe56_in_he4(30.5d0, 7.9d0, quietly)
76 2 : call do1_fe56_in_he4(26d0, 7d0, quietly)
77 2 : call do1_fe56_in_he4(24.5d0, 6.1d0, quietly)
78 2 : end subroutine test_fe56_in_he4
79 :
80 6 : subroutine do1_fe56_in_he4(log10_ne, log10_T, quietly)
81 : real(dp), intent(in) :: log10_ne, log10_T
82 : logical, intent(in) :: quietly
83 6 : real(dp) :: z
84 : integer :: ierr
85 : 1 format(a40, 3f16.7)
86 6 : ierr = 0
87 6 : z = eval_charge_of_Fe56_in_He4(log10_ne, log10_T, ierr)
88 6 : if (ierr /= 0) then
89 0 : write (*, 1) 'eval_charge_of_Fe56_in_He4 failed', log10_ne, log10_T
90 0 : call mesa_error(__FILE__, __LINE__)
91 : end if
92 6 : if (quietly) return
93 3 : write (*, 1) 'z for Fe56 in He4', z, log10_ne, log10_T
94 : end subroutine do1_fe56_in_he4
95 :
96 4 : subroutine do_test_eval_ionization(quietly)
97 : use ionization_def, only: num_ion_vals, ion_result_names
98 : logical, intent(in) :: quietly
99 4 : real(dp) :: Z, X, Rho, log10Rho, T, log10T
100 58 : real(dp) :: res(num_ion_vals)
101 : integer :: ierr, i
102 :
103 : include 'formats'
104 :
105 2 : ierr = 0
106 2 : Z = 0.0188d0
107 2 : X = 0.725d0
108 2 : T = 3.2345529591587989D+04
109 2 : log10T = log10(T)
110 2 : rho = 9.0768775206067858D-06
111 2 : log10Rho = log10(rho)
112 :
113 2 : if (.not. quietly) then
114 1 : write (*, *)
115 1 : write (*, *) 'do_test_eval_ionization'
116 1 : write (*, 1) 'log10T', log10T
117 1 : write (*, 1) 'log10rho', log10rho
118 1 : write (*, 1) 'Z', Z
119 1 : write (*, 1) 'X', X
120 : end if
121 :
122 2 : call eval_ionization(Z, X, Rho, log10Rho, T, log10T, res, ierr)
123 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
124 :
125 2 : if (quietly) return
126 :
127 1 : write (*, *)
128 29 : do i = 1, num_ion_vals
129 29 : write (*, 1) trim(ion_result_names(i)), res(i)
130 : end do
131 1 : write (*, *)
132 1 : write (*, *) 'done'
133 1 : write (*, *)
134 :
135 : end subroutine do_test_eval_ionization
136 :
137 0 : subroutine Build_Plots
138 : use utils_lib, only: mkdir
139 : character(len=256) :: data_dir, dir
140 :
141 0 : real(dp) :: log_ne_min, log_ne_max, logT_min, logT_max, dlog_ne, dlogT, log_ne, logT
142 :
143 : integer :: log_ne_points, logT_points
144 : integer :: i, j, k, info, io, io_first, io_last, io_log_ne, io_logT, num_vals
145 :
146 : integer, parameter :: io_unit0 = 40
147 :
148 : !real(dp) :: M, X, Z, kap
149 :
150 0 : real(dp), allocatable :: output_values(:, :, :)
151 :
152 0 : data_dir = '../../data'
153 :
154 : !..set the sample size
155 0 : log_ne_points = 300
156 0 : logT_points = 300
157 :
158 : !..set the ranges
159 :
160 0 : log_ne_max = 31
161 0 : log_ne_min = 24
162 0 : logT_min = 6
163 0 : logT_max = 8
164 :
165 : !..open the output files
166 0 : io_log_ne = io_unit0
167 0 : io_logT = io_unit0 + 2
168 0 : io_first = io_unit0 + 3
169 :
170 0 : dir = 'plot_data'
171 0 : call mkdir(dir)
172 0 : call Open_Plot_Outfiles(io_first, io_last, io_log_ne, io_logT, dir)
173 0 : num_vals = io_last - io_first + 1
174 0 : allocate (output_values(logT_points, log_ne_points, num_vals))
175 :
176 : !..get the results
177 :
178 0 : dlog_ne = (log_ne_max - log_ne_min)/(log_ne_points - 1)
179 0 : dlogT = (logT_max - logT_min)/(logT_points - 1)
180 :
181 0 : do j = 1, logT_points
182 0 : logT = logT_min + dlogT*(j - 1)
183 0 : do i = 1, log_ne_points
184 0 : log_ne = log_ne_min + dlog_ne*(i - 1)
185 : call Plot_one( &
186 : i, j, log_ne, logT, &
187 0 : output_values, num_vals, logT_points, log_ne_points, info)
188 : end do
189 : end do
190 :
191 0 : write (*, *) 'write files'
192 0 : do k = 1, num_vals
193 0 : write (*, *) k
194 0 : write (io_first + k - 1, '(e24.16)') output_values(1:logT_points, 1:log_ne_points, k)
195 : end do
196 :
197 0 : do i = 1, log_ne_points
198 0 : log_ne = log_ne_min + dlog_ne*(i - 1)
199 0 : write (io_log_ne, *) log_ne
200 : end do
201 0 : close (io_log_ne)
202 :
203 0 : do j = 1, logT_points
204 0 : logT = logT_min + dlogT*(j - 1)
205 0 : write (io_logT, *) logT
206 : end do
207 0 : close (io_logT)
208 :
209 0 : do io = io_first, io_last
210 0 : close (io)
211 : end do
212 :
213 0 : deallocate (output_values)
214 :
215 0 : end subroutine Build_Plots
216 :
217 0 : subroutine Plot_one( &
218 : i, j, log_ne, logT, &
219 0 : output_values, num_vals, logT_points, log_ne_points, ierr)
220 : integer, intent(in) :: i, j, num_vals, logT_points, log_ne_points
221 : real(dp), intent(in) :: log_ne, logT
222 : real(dp), intent(out) :: output_values(logT_points, log_ne_points, num_vals)
223 : integer, intent(out) :: ierr
224 :
225 0 : real(dp) :: z
226 : integer :: k
227 :
228 : include 'formats'
229 :
230 0 : ierr = 0
231 :
232 0 : z = eval_charge_of_Fe56_in_He4(log_ne, logT, ierr)
233 0 : if (ierr /= 0) then
234 0 : write (*, 1) 'eval_charge_of_Fe56_in_He4 failed log_ne, logT', log_ne, logT
235 0 : return
236 : end if
237 :
238 0 : k = 0
239 0 : k = k + 1; output_values(i, j, k) = z
240 :
241 : end subroutine Plot_one
242 :
243 0 : subroutine Open_Plot_Outfiles(io_first, io_last, io_log_ne, io_logT, dir)
244 : integer, intent(in) :: io_first, io_log_ne, io_logT
245 : integer, intent(out) :: io_last
246 : character(len=256), intent(in) :: dir
247 : character(len=256) :: fname
248 : integer :: io
249 :
250 0 : fname = trim(dir)//'/log_ne.data'
251 0 : open (unit=io_log_ne, file=trim(fname))
252 :
253 0 : fname = trim(dir)//'/logT.data'
254 0 : open (unit=io_logT, file=trim(fname))
255 :
256 0 : io = io_first - 1
257 :
258 0 : fname = trim(dir)//'/z_fe56_he4.data'
259 0 : io = io + 1; open (unit=io, file=trim(fname))
260 :
261 0 : io_last = io
262 :
263 0 : end subroutine Open_Plot_Outfiles
264 :
265 : end module test_ionization_support
|
