Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2011 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 ion_table_plot
21 :
22 : use const_def, only: dp
23 : use ion_tables_eval
24 : use math_lib
25 : use utils_lib, only: mesa_error, mkdir
26 :
27 : implicit none
28 :
29 : contains
30 :
31 0 : subroutine do_create_table_plot_files
32 :
33 : character (len=256) :: dir
34 :
35 0 : real(dp) :: lgT_min, lgT_max, lgRho_min, lgRho_max, dlgT, &
36 0 : dlgRho, lgRho, lgT, Rho, T, Z, X, lgQ_min, lgQ_max
37 :
38 : integer :: lgT_points, lgRho_points
39 : integer :: i, j, k, ierr, io, io_first, io_last, io_params, io_rho, io_tmp, num_vals
40 :
41 : integer, parameter :: io_unit0 = 40
42 :
43 0 : real(dp), allocatable :: output_values(:,:,:)
44 :
45 0 : Z = 0.018d0
46 0 : X = 0.72d0
47 :
48 : !..set the sample size
49 0 : lgT_points = 300
50 0 : lgRho_points = 300
51 :
52 : !lgT_points = 100
53 : !lgRho_points = 100
54 :
55 : !lgT_points = 2
56 : !lgRho_points = 2
57 :
58 : !..set the ranges
59 :
60 : ! check opal/scvh
61 : lgT_max = 7.7d0
62 : lgT_min = 2.0d0
63 : lgRho_min = -5d0
64 : lgRho_max = 5.5d0
65 :
66 : ! table full range
67 : lgT_max = 8.2d0
68 : lgT_min = 2.1d0
69 0 : lgQ_min = -10d0
70 0 : lgQ_max = 5.69d0
71 : lgRho_min = -9d0 ! lgQ_min + 2*lgT_min - 12
72 : lgRho_max = 8d0 ! lgQ_max + 2*lgT_max - 12
73 :
74 : ! test
75 0 : lgT_max = 7.5d0
76 0 : lgT_min = 3d0
77 0 : lgRho_max = 3d0
78 0 : lgRho_min = -7d0
79 :
80 0 : io_params = io_unit0
81 0 : io_rho = io_unit0+1
82 0 : io_tmp = io_unit0+2
83 0 : io_first = io_unit0+3
84 :
85 0 : dir = 'plot_data'
86 0 : call mkdir(dir)
87 0 : call open_plot_files(io_first, io_last, io_params, io_rho, io_tmp, dir)
88 0 : write(io_params, '(2(f10.6),2(i7))') Z, X, lgRho_points, lgT_points
89 0 : close(io_params)
90 0 : num_vals = io_last - io_first + 1
91 0 : allocate(output_values(lgRho_points,lgT_points,num_vals))
92 :
93 0 : dlgT = (lgT_max - lgT_min)/(lgT_points-1)
94 0 : dlgRho = (lgRho_max - lgRho_min)/(lgRho_points-1)
95 :
96 0 : do j=1, lgT_points
97 0 : lgT = lgT_min + dlgT*(j-1)
98 0 : T = exp10(lgT)
99 0 : do i=1,lgRho_points
100 0 : lgRho = lgRho_min + dlgRho*(i-1)
101 0 : Rho = exp10(lgRho)
102 : call plot_one( &
103 : i, j, Z, X, lgT, T, lgRho, Rho, output_values, &
104 0 : num_vals, lgRho_points, lgT_points, ierr)
105 0 : if (ierr /= 0) call mesa_error(__FILE__,__LINE__)
106 : end do
107 : end do
108 :
109 0 : write(*,*) 'write plot files'
110 0 : do k = 1, num_vals
111 0 : write(*,*) k
112 0 : write(io_first+k-1,'(e24.16)') output_values(1:lgRho_points,1:lgT_points,k)
113 : end do
114 :
115 0 : do i = 1, lgT_points
116 0 : lgT = lgT_min + dlgT*(i-1)
117 0 : write(io_tmp,*) lgT
118 : end do
119 0 : close(io_tmp)
120 :
121 0 : do j=1,lgRho_points
122 0 : lgRho = lgRho_min + dlgRho*(j-1)
123 0 : write(io_rho,*) lgRho
124 : end do
125 0 : close(io_rho)
126 :
127 0 : do io=io_first,io_last
128 0 : close(io)
129 : end do
130 :
131 0 : deallocate(output_values)
132 :
133 0 : end subroutine do_create_table_plot_files
134 :
135 :
136 0 : subroutine plot_one( &
137 0 : i, j, Z, X, lgT, T, lgRho, Rho, output_values, num_vals, lgRho_points, lgT_points, ierr)
138 : integer, intent(in) :: i, j, num_vals, lgRho_points, lgT_points
139 : real(dp), intent(in) :: Z, X, lgT, T, lgRho, Rho
140 : real(dp), intent(inout) :: output_values(lgRho_points,lgT_points,num_vals)
141 : integer, intent(out) :: ierr
142 :
143 0 : real(dp), dimension(num_ion_vals) :: res
144 : integer :: k
145 :
146 : include 'formats'
147 :
148 : ierr = 0
149 0 : call Get_ion_Results(Z, X, Rho, lgRho, T, lgT, res, ierr)
150 0 : if (ierr /= 0) then
151 0 : write(*,1) 'Get_ion_Results failed for lgRho, lgT', lgRho, lgT
152 0 : stop
153 : end if
154 :
155 0 : if (ierr == 0) then
156 :
157 0 : k = 0
158 0 : call save1('fneut_H', res(ion_ifneut_H))
159 0 : call save1('fneut_He', res(ion_ifneut_He))
160 0 : call save1('fneut_C', res(ion_ifneut_C))
161 0 : call save1('fneut_N', res(ion_ifneut_N))
162 0 : call save1('fneut_O', res(ion_ifneut_O))
163 0 : call save1('fneut_Ne', res(ion_ifneut_Ne))
164 0 : call save1('fneut_Mg', res(ion_ifneut_Mg))
165 0 : call save1('fneut_Si', res(ion_ifneut_Si))
166 0 : call save1('fneut_Fe', res(ion_ifneut_Fe))
167 0 : call save1('z_H', res(ion_iZ_H))
168 0 : call save1('z_He', res(ion_iZ_He))
169 0 : call save1('z_C', res(ion_iZ_C))
170 0 : call save1('z_N', res(ion_iZ_N))
171 0 : call save1('z_O', res(ion_iZ_O))
172 0 : call save1('z_Ne', res(ion_iZ_Ne))
173 0 : call save1('z_Mg', res(ion_iZ_Mg))
174 0 : call save1('z_Si', res(ion_iZ_Si))
175 0 : call save1('z_Fe', res(ion_iZ_Fe))
176 0 : call save1('logpp_H', res(ion_ilogpp_H))
177 0 : call save1('logpp_He', res(ion_ilogpp_He))
178 0 : call save1('logpp_C', res(ion_ilogpp_C))
179 0 : call save1('logpp_N', res(ion_ilogpp_N))
180 0 : call save1('logpp_O', res(ion_ilogpp_O))
181 0 : call save1('logpp_Ne', res(ion_ilogpp_Ne))
182 0 : call save1('logpp_Mg', res(ion_ilogpp_Mg))
183 0 : call save1('logpp_Si', res(ion_ilogpp_Si))
184 0 : call save1('logpp_Fe', res(ion_ilogpp_Fe))
185 :
186 : end if
187 :
188 : contains
189 :
190 0 : subroutine save1(str,v)
191 : character (len=*), intent(in) :: str
192 : real(dp), intent(in) :: v
193 0 : k = k+1; output_values(i,j,k) = v
194 0 : end subroutine save1
195 :
196 : end subroutine plot_one
197 :
198 :
199 0 : subroutine open_plot_files(io_first, io_last, io_params, io_rho, io_tmp, dir)
200 : integer, intent(IN) :: io_first, io_params, io_rho, io_tmp
201 : integer, intent(OUT) :: io_last
202 : character (len=256), intent(IN) :: dir
203 : character (len=256) :: fname
204 : integer :: io
205 :
206 0 : fname = trim(dir) // '/params.data'
207 0 : open(unit=io_params,file=trim(fname))
208 :
209 0 : fname = trim(dir) // '/logRho.data'
210 0 : open(unit=io_rho,file=trim(fname))
211 :
212 0 : fname = trim(dir) // '/logT.data'
213 0 : open(unit=io_tmp,file=trim(fname))
214 :
215 0 : io = io_first-1
216 : !call open1('logP')
217 : !call open1('logPgas')
218 : !call open1('logS')
219 : !call open1('chiRho')
220 : !call open1('chiT')
221 : !call open1('logCp')
222 : !call open1('logCv')
223 : !call open1('dlnE_dlnRho')
224 : !call open1('dlnS_dlnT')
225 : !call open1('dlnS_dlnRho')
226 : !call open1('mu')
227 : !call open1('log_free_e')
228 : !call open1('gamma1')
229 : !call open1('gamma3')
230 : !call open1('grad_ad')
231 : !call open1('eta')
232 0 : call open1('fneut_H')
233 0 : call open1('fneut_He')
234 0 : call open1('fneut_C')
235 0 : call open1('fneut_N')
236 0 : call open1('fneut_O')
237 0 : call open1('fneut_Ne')
238 0 : call open1('fneut_Mg')
239 0 : call open1('fneut_Si')
240 0 : call open1('fneut_Fe')
241 0 : call open1('z_H')
242 0 : call open1('z_He')
243 0 : call open1('z_C')
244 0 : call open1('z_N')
245 0 : call open1('z_O')
246 0 : call open1('z_Ne')
247 0 : call open1('z_Mg')
248 0 : call open1('z_Si')
249 0 : call open1('z_Fe')
250 0 : call open1('logpp_H')
251 0 : call open1('logpp_He')
252 0 : call open1('logpp_C')
253 0 : call open1('logpp_N')
254 0 : call open1('logpp_O')
255 0 : call open1('logpp_Ne')
256 0 : call open1('logpp_Mg')
257 0 : call open1('logpp_Si')
258 0 : call open1('logpp_Fe')
259 : !call open1('logE')
260 : !call open1('logW')
261 0 : io_last = io
262 :
263 : contains
264 :
265 0 : subroutine open1(name)
266 : character (len=*), intent(in) :: name
267 0 : fname = trim(dir) // '/' // trim(name) // '.data'
268 0 : io = io+1; open(unit=io,file=trim(fname))
269 0 : end subroutine open1
270 :
271 : end subroutine open_plot_files
272 :
273 : end module ion_table_plot
|
