Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015-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_rti
21 :
22 : use star_private_def
23 : use const_def, only: dp, Msun
24 : use pgstar_support
25 : use star_pgstar
26 :
27 : implicit none
28 :
29 : contains
30 :
31 0 : subroutine rti_Plot(id, device_id, ierr)
32 : integer, intent(in) :: id, device_id
33 : integer, intent(out) :: ierr
34 : type (star_info), pointer :: s
35 :
36 : ierr = 0
37 0 : call get_star_ptr(id, s, ierr)
38 0 : if (ierr /= 0) return
39 :
40 0 : call pgslct(device_id)
41 0 : call pgbbuf()
42 0 : call pgeras()
43 :
44 : call do_rti_Plot(s, id, device_id, &
45 : s% pg% rti_xleft, s% pg% rti_xright, &
46 : s% pg% rti_ybot, s% pg% rti_ytop, .false., &
47 0 : s% pg% rti_title, s% pg% rti_txt_scale, ierr)
48 0 : call pgebuf()
49 :
50 : end subroutine rti_Plot
51 :
52 :
53 0 : subroutine do_rti_Plot(s, id, device_id, &
54 : vp_xleft, vp_xright, vp_ybot, vp_ytop, subplot, title, txt_scale, ierr)
55 : use chem_def
56 : use net_def
57 : use utils_lib
58 :
59 : type (star_info), pointer :: s
60 : integer, intent(in) :: id, device_id
61 : real, intent(in) :: vp_xleft, vp_xright, vp_ybot, vp_ytop, txt_scale
62 : logical, intent(in) :: subplot
63 : character (len=*), intent(in) :: title
64 : integer, intent(out) :: ierr
65 :
66 : integer :: i, n, step_min, step_max
67 : real :: xmin, xmax, ymin_mass_axis, ymax_mass_axis
68 0 : real, allocatable, dimension(:) :: xvec, star_mass, star_M_center, log_xmstar, &
69 0 : he_core_mass, &
70 0 : c_core_mass, &
71 0 : o_core_mass, &
72 0 : si_core_mass, &
73 0 : fe_core_mass, &
74 0 : shock_mass
75 : logical :: have_star_mass, have_log_xmstar, &
76 : have_he_core_mass, &
77 : have_c_core_mass, &
78 : have_o_core_mass, &
79 : have_si_core_mass, &
80 : have_fe_core_mass, &
81 : have_shock_mass
82 :
83 : integer :: ix,k
84 : real :: xleft,xright,now
85 : real, parameter :: dxmin = -1
86 :
87 : include 'formats'
88 :
89 0 : ierr = 0
90 :
91 0 : step_min = s% pg% rti_step_xmin
92 0 : if (step_min <= 0) step_min = 1
93 0 : step_max = s% pg% rti_step_xmax
94 0 : if (step_max <= 0 .or. step_max > s% model_number) step_max = s% model_number
95 :
96 0 : if (step_min >= s% model_number) step_min = 1
97 :
98 0 : if (s% pg% rti_max_width > 0) &
99 0 : step_min = max(step_min, step_max - s% pg% rti_max_width)
100 :
101 0 : n = count_hist_points(s, step_min, step_max)
102 0 : if (n <= 1) return
103 :
104 0 : call integer_dict_lookup(s% history_names_dict, s% pg% rti_xaxis_name, ix, ierr)
105 0 : if (ierr /= 0) ix = -1
106 0 : if (ix <= 0) then
107 0 : write(*,'(A)')
108 : write(*,*) 'ERROR: failed to find ' // &
109 0 : trim(s% pg% rti_xaxis_name) // ' in rti data'
110 0 : write(*,'(A)')
111 0 : ierr = -1
112 0 : return
113 : end if
114 :
115 : allocate(xvec(n), &
116 : star_mass(n), &
117 : log_xmstar(n), &
118 : star_M_center(n), &
119 : he_core_mass(n), &
120 : c_core_mass(n), &
121 : o_core_mass(n), &
122 : si_core_mass(n), &
123 : fe_core_mass(n), &
124 : shock_mass(n), &
125 0 : stat=ierr)
126 0 : if (ierr /= 0) then
127 0 : write(*,*) 'allocate failed for PGSTAR rti'
128 0 : return
129 : end if
130 :
131 0 : call get_hist_points(s, step_min, step_max, n, ix, xvec, ierr)
132 0 : if (ierr /= 0) then
133 0 : write(*,*) 'pgstar get_hist_points failed ' // trim(s% pg% rti_xaxis_name)
134 0 : call dealloc
135 0 : ierr = 0
136 0 : return
137 : end if
138 :
139 0 : if (s% pg% rti_xaxis_in_seconds .and. s% pg% rti_xaxis_name=='star_age') then
140 0 : do k=1,n
141 0 : xvec(k) = xvec(k)*secyer
142 : end do
143 0 : else if (s% pg% rti_xaxis_in_Myr .and. s% pg% rti_xaxis_name=='star_age') then
144 0 : do k=1,n
145 0 : xvec(k) = xvec(k)*1d-6
146 : end do
147 : end if
148 :
149 0 : now=xvec(n)
150 0 : if (s% pg% rti_xaxis_time_from_present .and. s% pg% rti_xaxis_name=='star_age') then
151 0 : do k=1,n
152 0 : xvec(k) = xvec(k)-now
153 : end do
154 : end if
155 :
156 0 : if (s% pg% rti_xaxis_log) then
157 0 : do k=1,n
158 0 : xvec(k) = real(log10(max(1d-50,abs(dble(xvec(k))))))
159 : end do
160 : end if
161 :
162 0 : xmin = s% pg% rti_xmin
163 0 : xmax = s% pg% rti_xmax
164 :
165 : call set_xleft_xright( &
166 : n, xvec, xmin, xmax, s% pg% rti_xmargin, &
167 0 : s% pg% rti_xaxis_reversed, dxmin, xleft, xright)
168 :
169 0 : have_star_mass = get1_yvec('star_mass', star_mass)
170 0 : if (.not. have_star_mass) then
171 0 : write(*,*) 'failed to find star_mass in history data'
172 0 : ierr = -1
173 : end if
174 0 : if (ierr /= 0) then
175 0 : call dealloc
176 0 : return
177 : end if
178 :
179 0 : have_log_xmstar = get1_yvec('log_xmstar', log_xmstar)
180 0 : if (have_log_xmstar) then
181 0 : do i = 1, n
182 : star_M_center(i) = &
183 0 : star_mass(i) - real(exp10(dble(log_xmstar(i)))/Msun)
184 : end do
185 : else
186 0 : star_M_center(:) = 0
187 : end if
188 :
189 0 : have_he_core_mass = get1_yvec('he_core_mass', he_core_mass)
190 0 : have_c_core_mass = get1_yvec('c_core_mass', c_core_mass)
191 0 : have_o_core_mass = get1_yvec('o_core_mass', o_core_mass)
192 0 : have_si_core_mass = get1_yvec('si_core_mass', si_core_mass)
193 0 : have_fe_core_mass = get1_yvec('fe_core_mass', fe_core_mass)
194 0 : have_shock_mass = get1_yvec('shock_mass', shock_mass)
195 :
196 0 : call pgsave
197 0 : call pgsch(txt_scale)
198 :
199 0 : call init_rti_plot
200 0 : call setup_mass_yaxis
201 0 : call plot_total_mass_line
202 0 : call plot_rti_data
203 0 : if (ierr /= 0) then
204 0 : call pgunsa
205 0 : call dealloc
206 0 : return
207 : end if
208 0 : call plot_mass_lines
209 :
210 : call show_annotations(s, &
211 : s% pg% show_rti_annotation1, &
212 : s% pg% show_rti_annotation2, &
213 0 : s% pg% show_rti_annotation3)
214 :
215 0 : call pgsch(txt_scale)
216 0 : call finish_rti_plot
217 :
218 0 : call pgunsa
219 :
220 0 : call dealloc
221 :
222 :
223 : contains
224 :
225 0 : subroutine dealloc
226 0 : deallocate(xvec, star_mass, star_M_center, log_xmstar, &
227 0 : he_core_mass, &
228 0 : c_core_mass, &
229 0 : o_core_mass, &
230 0 : si_core_mass, &
231 0 : fe_core_mass, &
232 0 : shock_mass)
233 0 : end subroutine dealloc
234 :
235 0 : subroutine plot_total_mass_line
236 : use pgstar_colors
237 0 : call pgsave
238 0 : call pgsch(txt_scale*0.8)
239 0 : call pgsci(clr_Gray)
240 0 : call pgslw(2)
241 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,1.0,1.0,'M\dtotal\u',-0.8)
242 0 : call pgslw(s% pg% pgstar_lw)
243 0 : call plot_mass_line(star_mass)
244 0 : call pgunsa
245 0 : end subroutine plot_total_mass_line
246 :
247 :
248 0 : subroutine plot_mass_line(yvec)
249 : real, intent(in) :: yvec(:)
250 0 : if (any(yvec > 1e-2*ymax_mass_axis)) then
251 0 : call pgsave
252 0 : call pgslw(s% pg% rti_line_weight)
253 0 : call pgline(n, xvec, yvec)
254 0 : call pgunsa
255 : end if
256 0 : end subroutine plot_mass_line
257 :
258 :
259 0 : subroutine setup_mass_yaxis
260 : use pgstar_colors, only: clr_Foreground
261 : real :: dy, ymin, ymax
262 : include 'formats'
263 0 : ymax = s% pg% rti_mass_max
264 0 : if (ymax <= 0) ymax = maxval(star_mass)
265 0 : ymin = s% pg% rti_mass_min
266 0 : if (ymin < 0) ymin = minval(star_M_center)
267 0 : if (ymax <= ymin) ymax = ymin+1
268 0 : dy = ymax - ymin
269 0 : if (s% pg% rti_mass_min /= -101d0) ymin = ymin - s% pg% rti_mass_margin*dy
270 0 : if (s% pg% rti_mass_max /= -101d0) ymax = ymax + s% pg% rti_mass_margin*dy
271 0 : call pgswin(xleft, xright, ymin, ymax)
272 0 : call pgscf(1)
273 0 : call pgsci(clr_Foreground)
274 0 : ymin_mass_axis = ymin
275 0 : ymax_mass_axis = ymax
276 0 : end subroutine setup_mass_yaxis
277 :
278 :
279 0 : subroutine init_rti_plot
280 0 : call pgsvp(vp_xleft, vp_xright, vp_ybot, vp_ytop)
281 0 : end subroutine init_rti_plot
282 :
283 :
284 0 : subroutine finish_rti_plot
285 : character(len=256) :: xlabel
286 0 : call show_box_pgstar(s,'BCNST1','BCNMSTV1')
287 :
288 0 : xlabel=''
289 0 : if (s% pg% rti_xaxis_log) then
290 0 : xlabel='log '// s% pg% rti_xaxis_name
291 : else
292 0 : xlabel=s% pg% rti_xaxis_name
293 : end if
294 :
295 0 : if (s% pg% rti_xaxis_name =='star_age') then
296 0 : if (s% pg% rti_xaxis_in_seconds) then
297 0 : xlabel=trim(xlabel)//' (s)'
298 0 : else if (s% pg% rti_xaxis_in_Myr) then
299 0 : xlabel=trim(xlabel)//' (Myr)'
300 : else
301 0 : xlabel=trim(xlabel)//' (yr)'
302 : end if
303 : end if
304 :
305 0 : call show_xaxis_label_pgstar(s,trim(xlabel),1.0)
306 :
307 0 : call show_left_yaxis_label_pgstar(s,'M/M\d\(2281)')
308 0 : if (.not. subplot) then
309 0 : call show_model_number_pgstar(s)
310 0 : call show_age_pgstar(s)
311 : end if
312 0 : call show_title_pgstar(s, title)
313 :
314 : call show_pgstar_decorator(id, s% pg% rti_use_decorator, &
315 0 : s% pg% rti_pgstar_decorator, 0, ierr)
316 :
317 0 : end subroutine finish_rti_plot
318 :
319 :
320 0 : logical function get1_yvec(name, vec)
321 : character (len=*) :: name
322 : real, dimension(:), allocatable :: vec
323 0 : get1_yvec = get1_hist_yvec(s, step_min, step_max, n, name, vec)
324 0 : end function get1_yvec
325 :
326 :
327 0 : subroutine plot_rti_data
328 : use history_specs, only: rti_offset
329 : type (pgstar_hist_node), pointer :: pg
330 : integer :: i_rti_type_first, i_rti_type_last
331 : integer :: k, cnt, num_specs, step, j
332 :
333 : include 'formats'
334 :
335 0 : i_rti_type_first = 0
336 0 : num_specs = size(s% history_column_spec, dim=1)
337 0 : do k = 1, num_specs
338 0 : if (s% history_column_spec(k) == rti_offset + 1) then
339 0 : i_rti_type_first = k
340 0 : exit
341 : end if
342 : end do
343 0 : if (i_rti_type_first == 0) then
344 0 : write(*,*) "please add 'rti_regions <num>' to your history_columns.list for the rti plot"
345 0 : ierr = -1
346 0 : return
347 : end if
348 :
349 0 : i_rti_type_last = i_rti_type_first
350 0 : cnt = 1
351 0 : do k=i_rti_type_first+1, num_specs
352 0 : cnt = cnt+1
353 0 : if (s% history_column_spec(k) /= rti_offset + cnt) exit
354 0 : i_rti_type_last = k
355 : end do
356 :
357 0 : call pgsave
358 0 : call pgslw(s% pg% rti_line_weight)
359 0 : if (.not. associated(s% pg% pgstar_hist)) then
360 0 : write(*,*) '.not. associated(s% pg% pgstar_hist)'
361 : end if
362 0 : pg => s% pg% pgstar_hist
363 0 : j = n
364 :
365 0 : do
366 0 : if (.not. associated(pg)) exit
367 0 : step = pg% step
368 0 : if (step < step_min) exit
369 0 : if (step <= step_max) then
370 0 : if (j < 1) exit
371 0 : if (mod(step, max(1, s% pg% rti_interval)) == 0) then
372 : call draw_rti_for_step( &
373 0 : pg, i_rti_type_first, i_rti_type_last, xvec(j), star_mass(j), star_M_center(j))
374 : end if
375 0 : j = j-1
376 : end if
377 0 : pg => pg% next
378 : end do
379 0 : call pgunsa
380 :
381 : end subroutine plot_rti_data
382 :
383 :
384 0 : subroutine draw_rti_for_step( &
385 : pg, i_rti_type_first, i_rti_type_last, xval, mass, mass_center)
386 : use pgstar_colors
387 : type (pgstar_hist_node), pointer :: pg
388 : integer, intent(in) :: i_rti_type_first, i_rti_type_last
389 : real, intent(in) :: xval, mass, mass_center
390 : real :: qbot, qtop, mbot, mtop, xmass
391 : integer :: k, rti_type
392 : include 'formats'
393 0 : qbot = 0
394 0 : xmass = mass - mass_center
395 0 : do k = i_rti_type_first, i_rti_type_last-1, 2
396 0 : rti_type = int(pg% vals(k))
397 0 : if (rti_type < 0) exit
398 0 : qtop = pg% vals(k+1)
399 0 : if (rti_type > 0) then
400 0 : mbot = mass_center + xmass*qbot
401 0 : mtop = mass_center + xmass*qtop
402 0 : call draw1(xval, mbot, mtop, clr_Blue)
403 : end if
404 0 : qbot = qtop
405 0 : if (qbot >= 1) exit
406 : end do
407 0 : end subroutine draw_rti_for_step
408 :
409 :
410 0 : subroutine draw1(xval,y1,y2,clr)
411 : real, intent(in) :: xval, y1, y2
412 : integer, intent(in) :: clr
413 : real :: top, bot
414 0 : if (y1 < y2) then
415 0 : bot = y1; top = y2
416 : else
417 0 : bot = y2; top = y1
418 : end if
419 0 : if (top < 1d-50) return
420 0 : call pgsci(clr)
421 0 : call pgmove(xval, bot)
422 0 : call pgdraw(xval, top)
423 : end subroutine draw1
424 :
425 :
426 0 : subroutine plot_mass_lines
427 : use pgstar_colors
428 :
429 : include 'formats'
430 :
431 0 : call pgsave
432 0 : call pgsch(txt_scale*0.8)
433 0 : call pgslw(2)
434 :
435 0 : if (have_he_core_mass) then
436 0 : call pgsci(clr_Teal)
437 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.77,0.5,'He',-0.8)
438 0 : call plot_mass_line(he_core_mass)
439 : else
440 0 : write(*,*) "please add 'he_core_mass' to your history_columns.list"
441 : end if
442 :
443 0 : if (have_c_core_mass) then
444 0 : call pgsci(clr_LightOliveGreen)
445 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.59,0.5,'C',-0.8)
446 0 : call plot_mass_line(c_core_mass)
447 : else
448 0 : write(*,*) "please add 'c_core_mass' to your history_columns.list"
449 : end if
450 :
451 0 : if (have_o_core_mass) then
452 0 : call pgsci(clr_SeaGreen)
453 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.41,0.5,'O',-0.8)
454 0 : call plot_mass_line(o_core_mass)
455 : else
456 0 : write(*,*) "please add 'o_core_mass' to your history_columns.list"
457 : end if
458 :
459 0 : if (have_si_core_mass) then
460 0 : call pgsci(clr_Lilac)
461 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.23,0.5,'Si',-0.8)
462 0 : call plot_mass_line(si_core_mass)
463 : else
464 0 : write(*,*) "please add 'si_core_mass' to your history_columns.list"
465 : end if
466 :
467 0 : if (have_fe_core_mass) then
468 0 : call pgsci(clr_Crimson)
469 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.00,0.0,'Iron',-0.8)
470 0 : call plot_mass_line(fe_core_mass)
471 : else
472 0 : write(*,*) "please add 'fe_core_mass' to your history_columns.list"
473 : end if
474 :
475 0 : if (have_shock_mass) then
476 0 : call pgsci(clr_IndianRed)
477 0 : call show_right_yaxis_label_pgmtxt_pgstar(s,0.5,0.5,'Outward Shock',-0.8)
478 0 : call plot_mass_line(shock_mass)
479 : else
480 0 : write(*,*) "please add 'shock_mass' to your history_columns.list"
481 : end if
482 :
483 0 : call pgunsa
484 :
485 0 : end subroutine plot_mass_lines
486 :
487 : end subroutine do_rti_Plot
488 :
489 : end module pgstar_rti
|
