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 pgstar_kipp_residuals
21 :
22 : use const_def, only: dp, Msun, Rsun
23 : use star_private_def
24 : use star_pgstar
25 : use pgstar_colors
26 :
27 : implicit none
28 :
29 : private
30 : public :: Kipp_residuals_plot, do_Kipp_residuals_plot
31 :
32 : real(dp), parameter :: resid_floor = 1d-20
33 :
34 : ! Grid points with no data (e.g, because of mass loss)
35 : ! are not drawn and show the background color.
36 : ! The upper half of pgstar's own colormap (pgstar_colors), which
37 : ! runs from near-white to pure red. These indices are set once per device
38 : ! by set_device_colors, so no color representation is changed here (safe
39 : ! for PNG and other file devices). Small residuals -> white, large -> red.
40 : ! Grid points with no data (e.g, because of mass loss)
41 : ! are not drawn and show the background color.
42 : integer, parameter :: c1 = colormap_offset + (colormap_length + 1)/2
43 : integer, parameter :: c2 = colormap_offset + colormap_length
44 :
45 : contains
46 :
47 : ! Kippenhahn-style log10(max(|res|)) across all equations for a
48 : ! given space-time location
49 0 : subroutine Kipp_residuals_plot(id, device_id, ierr)
50 : integer, intent(in) :: id, device_id
51 : integer, intent(out) :: ierr
52 :
53 : type (star_info), pointer :: s
54 :
55 : ierr = 0
56 0 : call get_star_ptr(id, s, ierr)
57 0 : if (ierr /= 0) return
58 :
59 0 : call pgslct(device_id)
60 0 : call pgbbuf()
61 0 : call pgeras()
62 :
63 : call do_Kipp_residuals_plot(s, id, &
64 : s% pg% Kipp_residuals_xleft, s% pg% Kipp_residuals_xright, &
65 : s% pg% Kipp_residuals_ybot, s% pg% Kipp_residuals_ytop, .false., &
66 : s% pg% Kipp_residuals_title, s% pg% Kipp_residuals_txt_scale, &
67 0 : s% pg% Kipp_residuals_max_width, s% pg% Kipp_residuals_yaxis_name, ierr)
68 0 : if (ierr /= 0) return
69 :
70 0 : call pgebuf()
71 :
72 : end subroutine Kipp_residuals_plot
73 :
74 :
75 0 : subroutine do_Kipp_residuals_plot(s, id, &
76 : win_xleft, win_xright, win_ybot, win_ytop, subplot, title, txt_scale, &
77 : max_width, yaxis_name, ierr)
78 : type (star_info), pointer :: s
79 : integer, intent(in) :: id, max_width
80 : real, intent(in) :: win_xleft, win_xright, win_ybot, win_ytop, txt_scale
81 : logical, intent(in) :: subplot
82 : character(len=*), intent(in) :: title, yaxis_name
83 : integer, intent(out) :: ierr
84 :
85 : integer :: k, neq, initial_width, n, x_size, y_size, x_new
86 : integer, save :: init_model
87 :
88 : real, allocatable, dimension(:, :), save :: nr_resid_buf, nr_ycoord_buf
89 : integer, allocatable, dimension(:), save :: nr_model_buf, nr_zone_buf
90 :
91 0 : real, allocatable :: tmp_resid_buf(:, :), tmp_ycoord_buf(:, :)
92 0 : integer, allocatable :: tmp_model_buf(:), tmp_zone_buf(:)
93 :
94 0 : ierr = 0
95 0 : neq = size(s% equ, 1) ! number of equations per cell
96 :
97 0 : if (.not. allocated(nr_resid_buf)) then
98 0 : if (max_width < 0) then
99 : initial_width = 10
100 : else
101 0 : initial_width = max_width
102 : end if
103 : ! 1.5 times the current nz; grown below if remeshing increases nz
104 0 : allocate(nr_resid_buf(initial_width, int(1.5 * s% nz)))
105 0 : allocate(nr_ycoord_buf(initial_width, int(1.5 * s% nz)))
106 0 : allocate(nr_model_buf(initial_width))
107 0 : allocate(nr_zone_buf(initial_width))
108 0 : nr_resid_buf = 0.0
109 0 : nr_ycoord_buf = 0.0
110 0 : nr_model_buf = 0
111 0 : nr_zone_buf = 0 ! 0 zones = empty column, skipped when plotting
112 :
113 0 : init_model = s% model_number - 1
114 0 : else if (s% nz > 0.95 * size(nr_resid_buf, 2)) then
115 : ! remeshing increased nz beyond 95% of the current array size
116 0 : x_size = size(nr_resid_buf, 1)
117 0 : y_size = size(nr_resid_buf, 2)
118 0 : call move_alloc(nr_resid_buf, tmp_resid_buf)
119 0 : call move_alloc(nr_ycoord_buf, tmp_ycoord_buf)
120 :
121 0 : allocate(nr_resid_buf(x_size, 5 * s% nz))
122 0 : allocate(nr_ycoord_buf(x_size, 5 * s% nz))
123 0 : nr_resid_buf = 0.0
124 0 : nr_ycoord_buf = 0.0
125 :
126 0 : nr_resid_buf(:, 1:y_size) = tmp_resid_buf(:, 1:y_size)
127 0 : nr_ycoord_buf(:, 1:y_size) = tmp_ycoord_buf(:, 1:y_size)
128 :
129 0 : deallocate(tmp_resid_buf, tmp_ycoord_buf)
130 : end if
131 :
132 0 : if (max_width < 0) then
133 0 : n = s% model_number - init_model
134 : else
135 0 : n = min(s% model_number - init_model, max_width)
136 : end if
137 0 : if (n < 1) return
138 :
139 : ! grow dimension
140 0 : if (n > size(nr_resid_buf, 1)) then
141 0 : x_size = size(nr_resid_buf, 1)
142 0 : y_size = size(nr_resid_buf, 2)
143 0 : x_new = max(2 * x_size, n)
144 :
145 0 : call move_alloc(nr_resid_buf, tmp_resid_buf)
146 0 : call move_alloc(nr_ycoord_buf, tmp_ycoord_buf)
147 0 : call move_alloc(nr_zone_buf, tmp_zone_buf)
148 0 : call move_alloc(nr_model_buf, tmp_model_buf)
149 :
150 0 : allocate(nr_resid_buf(x_new, y_size))
151 0 : allocate(nr_ycoord_buf(x_new, y_size))
152 0 : allocate(nr_zone_buf(x_new))
153 0 : allocate(nr_model_buf(x_new))
154 :
155 0 : nr_resid_buf = 0.0
156 0 : nr_ycoord_buf = 0.0
157 0 : nr_zone_buf = 0
158 0 : nr_model_buf = 0
159 :
160 0 : nr_resid_buf(1:x_size, :) = tmp_resid_buf(1:x_size, :)
161 0 : nr_ycoord_buf(1:x_size, :) = tmp_ycoord_buf(1:x_size, :)
162 0 : nr_zone_buf(1:x_size) = tmp_zone_buf(1:x_size)
163 0 : nr_model_buf(1:x_size) = tmp_model_buf(1:x_size)
164 :
165 0 : deallocate(tmp_resid_buf, tmp_ycoord_buf, tmp_zone_buf, tmp_model_buf)
166 : end if
167 :
168 0 : if (max_width > 1 .and. s% model_number - init_model > max_width) then
169 0 : nr_model_buf(1:max_width-1) = nr_model_buf(2:max_width)
170 0 : nr_resid_buf(1:max_width-1, :) = nr_resid_buf(2:max_width, :)
171 0 : nr_ycoord_buf(1:max_width-1, :) = nr_ycoord_buf(2:max_width, :)
172 0 : nr_zone_buf(1:max_width-1) = nr_zone_buf(2:max_width)
173 : end if
174 :
175 0 : nr_model_buf(n) = s% model_number
176 0 : nr_zone_buf(n) = s% nz
177 :
178 : ! store the current model; all y coordinates are cell centered
179 0 : do k = 1, s% nz
180 : nr_resid_buf(n, k) = real(log10(max(resid_floor, &
181 0 : maxval(abs(s% equ(1:neq, k))))))
182 0 : select case (trim(yaxis_name))
183 : case ('logR')
184 0 : nr_ycoord_buf(n, k) = real(log10(max(1d-99, s% rmid(k)) / Rsun))
185 : case ('tau')
186 0 : nr_ycoord_buf(n, k) = real(log10(max(1d-99, s% tau(k))))
187 : case default ! 'mass'
188 0 : nr_ycoord_buf(n, k) = real((s% m(k) - 0.5d0 * s% dm(k)) / Msun)
189 : end select
190 : end do
191 :
192 0 : call pgsave
193 0 : call pgsci(0) ! color index 0 = background
194 0 : call pgsvp(win_xleft, win_xright, win_ybot, win_ytop)
195 0 : call pgrect(0.0, 1.0, 0.0, 1.0) ! pseudo erase
196 : call render(n, nr_model_buf, nr_zone_buf, nr_resid_buf, nr_ycoord_buf, &
197 : max(2, s% nz), yaxis_name, title, txt_scale, &
198 0 : s% pg% Kipp_residuals_min_residual, s% pg% Kipp_residuals_max_residual)
199 0 : call pgunsa
200 :
201 : end subroutine do_Kipp_residuals_plot
202 :
203 :
204 0 : subroutine render(nx, model_buf, zone_buf, resid_buf, ycoord_buf, &
205 : ny, yaxis_name, title, txt_scale, cbar_min, cbar_max)
206 : integer, intent(in) :: nx, ny
207 : integer, intent(in) :: model_buf(:), zone_buf(:)
208 : real, intent(in) :: resid_buf(:, :), ycoord_buf(:, :)
209 : character(len=*), intent(in) :: yaxis_name, title
210 : real, intent(in) :: txt_scale
211 : real, intent(in) :: cbar_min, cbar_max
212 0 : real, allocatable :: col(:, :), coord_grid(:)
213 : real :: tr(6), fg, bg, missing, xlo, xhi, xl, xr
214 : real :: clo, chi, dcoord, resid_hi, resid_lo
215 : integer :: i, j, j1, j2, nz_i, ip, cir_lo_save, cir_hi_save
216 : logical :: y_reverse, have_data
217 :
218 : ! --- data ranges over all stored columns ---
219 0 : have_data = .false.
220 0 : chi = -huge(1.0)
221 0 : clo = huge(1.0)
222 0 : resid_hi = -huge(1.0); resid_lo = huge(1.0)
223 0 : do i = 1, nx
224 0 : nz_i = zone_buf(i)
225 0 : if (nz_i < 1) cycle
226 0 : have_data = .true.
227 0 : chi = max(chi, maxval(ycoord_buf(i, 1:nz_i)))
228 0 : clo = min(clo, minval(ycoord_buf(i, 1:nz_i)))
229 0 : resid_hi = max(resid_hi, maxval(resid_buf(i, 1:nz_i)))
230 0 : resid_lo = min(resid_lo, minval(resid_buf(i, 1:nz_i)))
231 : end do
232 0 : if (.not. have_data) return
233 0 : if (chi <= clo) chi = clo + 1.0
234 :
235 0 : y_reverse = (trim(yaxis_name) == 'tau')
236 :
237 : ! uniform y grid, ascending with j; axis reversal is done in pgswin
238 0 : allocate(coord_grid(ny), col(1, ny))
239 0 : dcoord = (chi - clo) / real(ny - 1)
240 0 : do j = 1, ny
241 0 : coord_grid(j) = clo + real(j - 1) * dcoord
242 : end do
243 :
244 : ! --- color scale ---
245 0 : if (cbar_max /= -101.0) then
246 0 : fg = cbar_max
247 : else
248 0 : fg = resid_hi
249 : end if
250 0 : if (cbar_min /= -101.0) then
251 0 : bg = cbar_min
252 : else
253 0 : bg = resid_lo
254 : end if
255 0 : if (fg - bg < 1.0) then
256 0 : if (cbar_min /= -101.0 .and. cbar_max == -101.0) then
257 0 : fg = bg + 1.0
258 : else
259 0 : bg = fg - 1.0
260 : end if
261 : end if
262 0 : missing = bg - 1.0 ! marker only, never drawn
263 :
264 0 : call pgqcir(cir_lo_save, cir_hi_save)
265 0 : call pgscir(c1, c2)
266 :
267 : ! --- x range from stored model numbers ---
268 0 : xlo = huge(1.0); xhi = -huge(1.0)
269 0 : ip = 0
270 0 : do i = 1, nx
271 0 : if (zone_buf(i) < 1) cycle
272 0 : call column_edges(i, ip, nx, model_buf, zone_buf, xl, xr)
273 0 : xlo = min(xlo, xl)
274 0 : xhi = max(xhi, xr)
275 0 : ip = i
276 : end do
277 :
278 0 : if (.not. y_reverse) then
279 0 : call pgswin(xlo, xhi, clo - 0.5*dcoord, chi + 0.5*dcoord)
280 : else
281 0 : call pgswin(xlo, xhi, chi + 0.5*dcoord, clo - 0.5*dcoord)
282 : end if
283 :
284 : ! fill the plot area with grey first: points without data (e.g. above
285 : ! the surface after mass loss) stay grey, distinct from white (small
286 : ! residuals) and from the background, whether black or white
287 0 : call pgsci(clr_Gray)
288 0 : call pgsfs(1)
289 0 : call pgrect(xlo, xhi, clo - 0.5*dcoord, chi + 0.5*dcoord)
290 :
291 : ! --- one pgimag call per stored model, placed at its model number ---
292 0 : tr(3) = 0.0
293 0 : tr(4) = clo - dcoord ! j = 1 -> clo
294 0 : tr(5) = 0.0
295 0 : tr(6) = dcoord
296 0 : ip = 0
297 0 : do i = 1, nx
298 0 : nz_i = zone_buf(i)
299 0 : if (nz_i < 1) cycle
300 0 : call column_edges(i, ip, nx, model_buf, zone_buf, xl, xr)
301 0 : ip = i
302 : call interp_column(nz_i, ycoord_buf(i, 1:nz_i), resid_buf(i, 1:nz_i), &
303 0 : ny, coord_grid, dcoord, bg, fg, missing, col(1, :))
304 0 : tr(2) = xr - xl
305 0 : tr(1) = 0.5*(xl + xr) - tr(2) ! array index 1 -> column center
306 : ! the column coordinate is monotone, so missing points can only be
307 : ! at the ends; draw just the contiguous range with data
308 0 : j1 = 1
309 0 : do while (j1 <= ny)
310 0 : if (col(1, j1) > missing) exit
311 0 : j1 = j1 + 1
312 : end do
313 0 : j2 = ny
314 0 : do while (j2 >= j1)
315 0 : if (col(1, j2) > missing) exit
316 0 : j2 = j2 - 1
317 : end do
318 0 : if (j1 > j2) cycle
319 0 : call pgimag(col, 1, ny, 1, 1, j1, j2, bg, fg, tr)
320 : end do
321 :
322 : ! --- axes and labels ---
323 0 : call pgsci(clr_Foreground)
324 0 : call pgsch(txt_scale)
325 0 : call pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0)
326 0 : call pgmtxt('B', 2.5, 0.5, 0.5, 'Model Number')
327 0 : select case (trim(yaxis_name))
328 : case ('logR')
329 0 : call pgmtxt('L', 3.0, 0.5, 0.5, 'log(R/R\d\(2281)\u)')
330 : case ('tau')
331 0 : call pgmtxt('L', 3.0, 0.5, 0.5, 'log(\gt)')
332 : case default
333 0 : call pgmtxt('L', 3.0, 0.5, 0.5, 'm/M\d\(2281)\u')
334 : end select
335 0 : call pgmtxt('T', 1.0, 0.5, 0.5, trim(title))
336 :
337 : ! --- color wedge ---
338 0 : call pgwedg('RI', 0.5, 3.0, bg, fg, 'log\d10\u max |res|')
339 :
340 : ! restore the index range; no color representations were modified
341 0 : call pgscir(cir_lo_save, cir_hi_save)
342 :
343 0 : deallocate(col, coord_grid)
344 :
345 : end subroutine render
346 :
347 :
348 : ! x extent of column i: halfway to the neighboring non-empty columns
349 0 : subroutine column_edges(i, iprev, nx, model_buf, zone_buf, xl, xr)
350 : integer, intent(in) :: i, iprev, nx
351 : integer, intent(in) :: model_buf(:), zone_buf(:)
352 : real, intent(out) :: xl, xr
353 : integer :: inext, j
354 : real :: xc
355 :
356 0 : xc = real(model_buf(i))
357 0 : inext = 0
358 0 : do j = i + 1, nx
359 0 : if (zone_buf(j) > 0) then
360 : inext = j
361 : exit
362 : end if
363 : end do
364 :
365 0 : if (iprev > 0) then
366 0 : xl = 0.5 * (real(model_buf(iprev)) + xc)
367 0 : else if (inext > 0) then
368 0 : xl = xc - 0.5 * (real(model_buf(inext)) - xc)
369 : else
370 0 : xl = xc - 0.5
371 : end if
372 0 : if (inext > 0) then
373 0 : xr = 0.5 * (real(model_buf(inext)) + xc)
374 0 : else if (iprev > 0) then
375 0 : xr = xc + 0.5 * (xc - real(model_buf(iprev)))
376 : else
377 0 : xr = xc + 0.5
378 : end if
379 0 : end subroutine column_edges
380 :
381 0 : subroutine interp_column(m, ycol, rcol, ny, coord_grid, dcoord, bg, fg, missing, out)
382 : integer, intent(in) :: m, ny
383 : real, intent(in) :: ycol(m), rcol(m), coord_grid(ny), dcoord, bg, fg, missing
384 : real, intent(out) :: out(ny)
385 0 : real :: yc(m), rv(m), y, v, dy
386 : integer :: j, kk
387 :
388 0 : if (ycol(1) <= ycol(m)) then
389 0 : yc = ycol; rv = rcol
390 : else
391 0 : yc = ycol(m:1:-1); rv = rcol(m:1:-1)
392 : end if
393 :
394 : kk = 1
395 0 : do j = 1, ny
396 0 : y = coord_grid(j)
397 0 : if (y < yc(1)) then
398 : ! allow up to one grid step beyond the outermost cell center
399 0 : if (y < yc(1) - dcoord) then
400 0 : out(j) = missing
401 0 : cycle
402 : end if
403 0 : v = rv(1)
404 0 : else if (y > yc(m)) then
405 0 : if (y > yc(m) + dcoord) then
406 0 : out(j) = missing
407 0 : cycle
408 : end if
409 0 : v = rv(m)
410 0 : else if (m == 1) then
411 0 : v = rv(1)
412 : else
413 0 : do while (kk < m - 1)
414 0 : if (yc(kk + 1) >= y) exit
415 0 : kk = kk + 1
416 : end do
417 0 : dy = yc(kk + 1) - yc(kk)
418 0 : if (dy > 0.0) then
419 0 : v = rv(kk) + (y - yc(kk)) / dy * (rv(kk + 1) - rv(kk))
420 : else
421 0 : v = 0.5 * (rv(kk) + rv(kk + 1))
422 : end if
423 : end if
424 0 : out(j) = min(fg, max(bg, v))
425 : end do
426 0 : end subroutine interp_column
427 :
428 : end module pgstar_kipp_residuals
|
