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
21 :
22 : use star_private_def
23 : use const_def, only: dp, msun, anonymous_mixing, minimum_mixing, rayleigh_taylor_mixing
24 : use pgstar_support
25 : use star_pgstar
26 : use pgstar_colors
27 :
28 : implicit none
29 : private
30 :
31 : public :: Kipp_Plot, do_Kipp_Plot
32 :
33 : contains
34 :
35 0 : subroutine Kipp_Plot(id, device_id, ierr)
36 : integer, intent(in) :: id, device_id
37 : integer, intent(out) :: ierr
38 : type (star_info), pointer :: s
39 :
40 : ierr = 0
41 0 : call get_star_ptr(id, s, ierr)
42 0 : if (ierr /= 0) return
43 :
44 0 : call pgslct(device_id)
45 0 : call pgbbuf()
46 0 : call pgeras()
47 :
48 : call do_Kipp_Plot(s, id, device_id, &
49 : s% pg% Kipp_xleft, s% pg% Kipp_xright, &
50 : s% pg% Kipp_ybot, s% pg% Kipp_ytop, .false., &
51 0 : s% pg% Kipp_title, s% pg% Kipp_txt_scale, ierr)
52 0 : if (ierr /= 0) return
53 :
54 0 : call pgebuf()
55 :
56 : end subroutine Kipp_Plot
57 :
58 :
59 0 : subroutine do_Kipp_Plot(s, id, device_id, &
60 : vp_xleft, vp_xright, vp_ybot, vp_ytop, subplot, title, txt_scale, ierr)
61 : use chem_def
62 : use net_def
63 : use utils_lib
64 :
65 : type (star_info), pointer :: s
66 : integer, intent(in) :: id, device_id
67 : real, intent(in) :: vp_xleft, vp_xright, vp_ybot, vp_ytop, txt_scale
68 : logical, intent(in) :: subplot
69 : character (len=*), intent(in) :: title
70 : integer, intent(out) :: ierr
71 :
72 : integer :: i, n, step_min, step_max
73 0 : real :: xmin, xmax, ymin_L_axis, ymax_L_axis, &
74 0 : ymin_mass_axis, ymax_mass_axis, dx, burn_type_cutoff
75 0 : real, allocatable, dimension(:) :: xvec, &
76 0 : log_L, &
77 0 : log_Lneu, &
78 0 : log_LH, &
79 0 : log_LHe, &
80 0 : star_mass, &
81 0 : log_xmstar, &
82 0 : star_M_center, &
83 0 : he_core_mass, &
84 0 : c_core_mass, &
85 0 : o_core_mass, &
86 0 : si_core_mass, &
87 0 : fe_core_mass
88 : logical :: &
89 : have_log_L, &
90 : have_log_Lneu, &
91 : have_log_LH, &
92 : have_log_LHe, &
93 : have_star_mass, &
94 : have_log_xmstar, &
95 : have_he_core_mass, &
96 : have_c_core_mass, &
97 : have_o_core_mass, &
98 : have_si_core_mass, &
99 : have_fe_core_mass
100 : integer, parameter :: max_mix_type = leftover_convective_mixing
101 : integer :: mix_clr(max_mix_type), max_mix_type_to_show
102 : logical :: showed_this_mix_type(max_mix_type)
103 :
104 : integer :: ix,k
105 0 : real :: xleft,xright,now
106 : real, save :: dxmin = -1.d0
107 :
108 : include 'formats'
109 :
110 0 : ierr = 0
111 :
112 0 : mix_clr(convective_mixing) = clr_convection
113 0 : mix_clr(overshoot_mixing) = clr_overshoot
114 0 : mix_clr(semiconvective_mixing) = clr_semiconvection
115 0 : mix_clr(thermohaline_mixing) = clr_thermohaline
116 0 : mix_clr(rotation_mixing) = clr_rotation
117 0 : mix_clr(rayleigh_taylor_mixing) = clr_rayleigh_taylor
118 0 : mix_clr(minimum_mixing) = clr_minimum
119 0 : mix_clr(anonymous_mixing) = clr_anonymous
120 0 : mix_clr(leftover_convective_mixing) = clr_leftover_convection
121 :
122 0 : max_mix_type_to_show = thermohaline_mixing
123 0 : showed_this_mix_type = .false.
124 :
125 0 : step_min = s% pg% Kipp_step_xmin
126 0 : if (step_min <= 0) step_min = 1
127 0 : step_max = s% pg% Kipp_step_xmax
128 0 : if (step_max <= 0 .or. step_max > s% model_number) step_max = s% model_number
129 :
130 0 : if (step_min >= s% model_number) step_min = 1
131 :
132 0 : if (s% pg% Kipp_max_width > 0) &
133 0 : step_min = max(step_min, step_max - s% pg% Kipp_max_width)
134 :
135 0 : n = count_hist_points(s, step_min, step_max)
136 0 : if (n <= 1) return
137 0 : step_min = max(step_min, step_max-n+1)
138 :
139 0 : call integer_dict_lookup(s% history_names_dict, s% pg% kipp_xaxis_name, ix, ierr)
140 0 : if (ierr /= 0) ix = -1
141 0 : if (ix <= 0) then
142 0 : write(*,'(A)')
143 : write(*,*) 'ERROR: failed to find ' // &
144 0 : trim(s% pg% kipp_xaxis_name) // ' in kipp data'
145 0 : write(*,'(A)')
146 0 : ierr = -1
147 : end if
148 :
149 : allocate(xvec(n), &
150 : log_L(n), &
151 : log_Lneu(n), &
152 : log_LH(n), &
153 : log_LHe(n), &
154 : star_mass(n), &
155 : log_xmstar(n), &
156 : star_M_center(n), &
157 : he_core_mass(n), &
158 : c_core_mass(n), &
159 : o_core_mass(n), &
160 : si_core_mass(n), &
161 : fe_core_mass(n), &
162 0 : stat=ierr)
163 0 : if (ierr /= 0) then
164 0 : write(*,*) 'allocate failed for PGSTAR Kipp'
165 0 : return
166 : end if
167 :
168 0 : call get_hist_points(s, step_min, step_max, n, ix, xvec, ierr)
169 0 : if (ierr /= 0) then
170 0 : write(*,*) 'pgstar get_hist_points failed ' // trim(s% pg% kipp_xaxis_name)
171 0 : call dealloc
172 0 : ierr = 0
173 0 : return
174 : end if
175 :
176 0 : if (s% pg% kipp_xaxis_in_seconds .and. s% pg% kipp_xaxis_name=='star_age')THEN
177 0 : do k=1,n
178 0 : xvec(k) = xvec(k)*secyer
179 : end do
180 0 : else if (s% pg% kipp_xaxis_in_Myr .and. s% pg% kipp_xaxis_name=='star_age')THEN
181 0 : do k=1,n
182 0 : xvec(k) = xvec(k)*1d-6
183 : end do
184 : end if
185 :
186 0 : now=xvec(n)
187 0 : if (s% pg% kipp_xaxis_time_from_present .and. s% pg% kipp_xaxis_name=='star_age') then
188 0 : do k=1,n
189 0 : xvec(k) = xvec(k)-now
190 : end do
191 : end if
192 :
193 0 : if (s% pg% kipp_xaxis_log) then
194 0 : do k=1,n
195 0 : xvec(k) = log10(max(tiny(xvec(k)),abs(xvec(k))))
196 : end do
197 : end if
198 :
199 0 : if(s% pg% kipp_xmin<-100d0) s% pg% kipp_xmin=xvec(1)
200 0 : if(s% pg% kipp_xmax<-100d0) s% pg% kipp_xmax=xvec(n)
201 :
202 0 : xmin=max(s% pg% kipp_xmin,xvec(1))
203 0 : xmax=min(s% pg% kipp_xmax,xvec(n))
204 :
205 0 : burn_type_cutoff = s% pg% Kipp_burn_type_cutoff
206 :
207 : call set_xleft_xright( &
208 : n, xvec, xmin, xmax, s% pg% kipp_xmargin, &
209 : s% pg% kipp_xaxis_reversed, dxmin, xleft, xright)
210 :
211 0 : have_star_mass = get1_yvec('star_mass', star_mass)
212 0 : if (.not. have_star_mass) then
213 0 : write(*,*) 'PGSTAR Kipp failed to find star_mass in history data'
214 0 : ierr = -1
215 : end if
216 0 : if (ierr /= 0) return
217 0 : have_log_xmstar = get1_yvec('log_xmstar', log_xmstar)
218 0 : if (have_log_xmstar) then
219 0 : do i = 1, n
220 : star_M_center(i) = &
221 0 : star_mass(i) - real(exp10(dble(log_xmstar(i)))/Msun)
222 : end do
223 : else
224 0 : star_M_center(:) = 0
225 : end if
226 :
227 0 : if (s% pg% Kipp_show_luminosities) then
228 0 : have_log_L = get1_yvec('log_L', log_L)
229 0 : have_log_Lneu = get1_yvec('log_Lneu', log_Lneu)
230 0 : have_log_LH = get1_yvec('log_LH', log_LH)
231 0 : have_log_LHe = get1_yvec('log_LHe', log_LHe)
232 : else
233 0 : have_log_L = .false.
234 0 : have_log_Lneu = .false.
235 0 : have_log_LH = .false.
236 0 : have_log_LHe = .false.
237 : end if
238 :
239 0 : if (s% pg% Kipp_show_mass_boundaries) then
240 0 : have_he_core_mass = get1_yvec('he_core_mass', he_core_mass)
241 0 : have_c_core_mass = get1_yvec('c_core_mass', c_core_mass)
242 0 : have_o_core_mass = get1_yvec('o_core_mass', o_core_mass)
243 0 : have_si_core_mass = get1_yvec('si_core_mass', si_core_mass)
244 0 : have_fe_core_mass = get1_yvec('fe_core_mass', fe_core_mass)
245 : else
246 0 : have_he_core_mass = .false.
247 0 : have_c_core_mass = .false.
248 0 : have_o_core_mass = .false.
249 0 : have_si_core_mass = .false.
250 0 : have_fe_core_mass = .false.
251 : end if
252 :
253 0 : call pgsave
254 0 : call pgsch(txt_scale)
255 :
256 0 : dx = (xmax - xmin)/250.0
257 0 : call init_Kipp_plot
258 0 : call setup_mass_yaxis
259 0 : call plot_total_mass_line
260 0 : if (s% pg% Kipp_show_burn) then
261 0 : call plot_burn_data(dx)
262 0 : if (s% pg% Kipp_show_mixing) call plot_mix_data
263 0 : else if (s% pg% Kipp_show_mixing) then
264 0 : call plot_mix_data
265 : end if
266 0 : if (s% pg% Kipp_show_mass_boundaries) call plot_mass_lines
267 0 : if (s% pg% Kipp_show_luminosities) call plot_L_lines
268 :
269 : call show_annotations(s, &
270 : s% pg% show_Kipp_annotation1, &
271 : s% pg% show_Kipp_annotation2, &
272 0 : s% pg% show_Kipp_annotation3)
273 :
274 0 : call pgsch(txt_scale)
275 0 : call finish_Kipp_plot
276 :
277 0 : call pgunsa
278 :
279 0 : call dealloc
280 :
281 :
282 : contains
283 :
284 0 : subroutine dealloc
285 0 : deallocate(xvec, &
286 0 : log_L, &
287 0 : log_Lneu, &
288 0 : log_LH, &
289 0 : log_LHe, &
290 0 : star_mass, &
291 0 : log_xmstar, &
292 0 : star_M_center, &
293 0 : he_core_mass, &
294 0 : c_core_mass, &
295 0 : o_core_mass, &
296 0 : si_core_mass, &
297 0 : fe_core_mass)
298 0 : end subroutine dealloc
299 :
300 0 : subroutine init_Kipp_plot
301 0 : call pgsvp(vp_xleft, vp_xright, vp_ybot, vp_ytop)
302 0 : end subroutine init_Kipp_plot
303 :
304 :
305 0 : subroutine finish_Kipp_plot
306 : character(len=256) :: xlabel
307 0 : if (s% pg% Kipp_show_luminosities) then
308 0 : call setup_L_yaxis
309 0 : call show_box_pgstar(s,'','CMSTV')
310 0 : call show_right_yaxis_label_pgstar(s,'log (L\d\(2281)\u)')
311 : end if
312 0 : call setup_mass_yaxis
313 0 : if (s% pg% Kipp_show_luminosities) then
314 0 : call show_box_pgstar(s,'BCNST1','BNSTV1')
315 : else
316 0 : call show_box_pgstar(s,'BCNST1','BCNMSTV1')
317 : end if
318 :
319 0 : xlabel=''
320 0 : if (s% pg% kipp_xaxis_log) then
321 0 : xlabel='log '// s% pg% kipp_xaxis_name
322 : else
323 0 : xlabel=s% pg% kipp_xaxis_name
324 : end if
325 :
326 0 : if (s% pg% kipp_xaxis_name =='star_age') then
327 0 : if (s% pg% kipp_xaxis_in_seconds) then
328 0 : xlabel=trim(xlabel)//' (s)'
329 0 : else if (s% pg% Kipp_xaxis_in_Myr) then
330 0 : xlabel=trim(xlabel)//' (Myr)'
331 : else
332 0 : xlabel=trim(xlabel)//' (yr)'
333 : end if
334 : end if
335 :
336 0 : call show_xaxis_label_pgstar(s,trim(xlabel),1.0)
337 :
338 0 : call show_left_yaxis_label_pgstar(s,'M/M\d\(2281)')
339 0 : if (.not. subplot) then
340 0 : call show_model_number_pgstar(s)
341 0 : call show_age_pgstar(s)
342 : end if
343 0 : call show_title_pgstar(s, title)
344 0 : call show_mix_legend
345 0 : call show_burn_legend
346 :
347 0 : call show_pgstar_decorator(s%id, s% pg% kipp_use_decorator,s% pg% kipp_pgstar_decorator, 0, ierr)
348 :
349 :
350 0 : end subroutine finish_Kipp_plot
351 :
352 :
353 0 : logical function get1_yvec(name, vec)
354 : character (len=*) :: name
355 : real, dimension(:), allocatable :: vec
356 0 : get1_yvec = get1_hist_yvec(s, step_min, step_max, n, name, vec)
357 0 : end function get1_yvec
358 :
359 :
360 0 : subroutine setup_L_yaxis
361 0 : real :: dy, ymin, ymax
362 0 : ymin = -2
363 0 : ymax = ymin
364 0 : if (have_log_L) ymax = maxval(log_L)
365 0 : if (have_log_LH) ymax = max(ymax, maxval(log_LH))
366 0 : if (have_log_LHe) ymax = max(ymax, maxval(log_LHe))
367 0 : if (ymax <= ymin) ymax = ymin+1
368 0 : if (s% pg% Kipp_lgL_min /= -101d0) ymin = s% pg% Kipp_lgL_min
369 0 : if (s% pg% Kipp_lgL_max /= -101d0) ymax = s% pg% Kipp_lgL_max
370 0 : dy = ymax - ymin
371 0 : if (s% pg% Kipp_lgL_min /= -101d0) ymin = ymin - s% pg% Kipp_lgL_margin*dy
372 0 : if (s% pg% Kipp_lgL_max /= -101d0) ymax = ymax + s% pg% Kipp_lgL_margin*dy
373 0 : call pgswin(xleft, xright, ymin, ymax)
374 0 : call pgscf(1)
375 0 : call pgsci(clr_Foreground)
376 0 : ymin_L_axis = ymin
377 0 : ymax_L_axis = ymax
378 0 : end subroutine setup_L_yaxis
379 :
380 :
381 0 : subroutine setup_mass_yaxis
382 0 : real :: dy, ymin, ymax
383 0 : ymax = s% pg% Kipp_mass_max
384 0 : if (ymax <= 0) ymax = maxval(star_mass)
385 0 : ymin = s% pg% Kipp_mass_min
386 0 : if (ymin < 0) ymin = minval(star_M_center)
387 0 : if (ymax <= ymin) ymax = ymin+1
388 0 : dy = ymax - ymin
389 0 : if (s% pg% Kipp_mass_min /= -101d0) ymin = ymin - s% pg% Kipp_mass_margin*dy
390 0 : if (s% pg% Kipp_mass_max /= -101d0) ymax = ymax + s% pg% Kipp_mass_margin*dy
391 0 : call pgswin(xleft, xright, ymin, ymax)
392 0 : call pgscf(1)
393 0 : call pgsci(clr_Foreground)
394 0 : ymin_mass_axis = ymin
395 0 : ymax_mass_axis = ymax
396 0 : end subroutine setup_mass_yaxis
397 :
398 :
399 0 : subroutine plot_burn_data(dx)
400 : use history_specs, only: burning_offset
401 : real, intent(in) :: dx
402 : type (pgstar_hist_node), pointer :: pg
403 0 : real :: burn_max, burn_min, burn_type, x, xnext
404 : integer :: i_burn_type_first, i_burn_type_last
405 : integer :: k, cnt, num_specs, step
406 :
407 : include 'formats'
408 :
409 0 : if (.not. s% pg% Kipp_show_burn) return
410 0 : i_burn_type_first = 0
411 0 : num_specs = size(s% history_column_spec, dim=1)
412 0 : do k = 1, num_specs
413 0 : if (s% history_column_spec(k) == burning_offset + 1) then
414 0 : i_burn_type_first = k
415 0 : exit
416 : end if
417 : end do
418 0 : if (i_burn_type_first == 0) return
419 :
420 0 : i_burn_type_last = 0
421 0 : cnt = 1
422 0 : do k=i_burn_type_first+1, num_specs
423 0 : i_burn_type_last = k-1
424 0 : cnt = cnt+1
425 0 : if (s% history_column_spec(k) /= burning_offset + cnt) exit
426 : end do
427 :
428 0 : burn_max = 0.1
429 0 : burn_min = -0.1
430 0 : pg => s% pg% pgstar_hist
431 0 : do
432 0 : if (.not. associated(pg)) exit
433 0 : step = pg% step
434 0 : if (step < step_min) exit
435 0 : if (step <= step_max) then
436 0 : do k = i_burn_type_first, i_burn_type_last, 2
437 0 : burn_type = pg% vals(k)
438 0 : if (burn_type < -9990) exit
439 0 : if (burn_type /= 0) then
440 0 : if (burn_type > 0) then
441 0 : burn_max = max(burn_max, burn_type)
442 0 : else if (burn_type < -1) then
443 0 : burn_min = min(burn_min, burn_type)
444 : end if
445 : end if
446 0 : if (pg% vals(k+1) >= 1d0) exit
447 : end do
448 :
449 : end if
450 0 : pg => pg% next
451 : end do
452 :
453 0 : call pgsave
454 0 : call pgslw(s% pg% Kipp_burn_line_weight)
455 0 : pg => s% pg% pgstar_hist
456 0 : xnext = xmax
457 0 : do
458 0 : if (.not. associated(pg)) exit
459 0 : step = pg% step
460 0 : if (step < step_min) exit
461 0 : x = real(xvec(step-step_min+1))
462 0 : if (step <= step_max .and. x <= xnext) then
463 : call draw_burn_for_step( &
464 : pg, i_burn_type_first, i_burn_type_last, &
465 : burn_max, burn_min, burn_type_cutoff, x, &
466 0 : star_mass(step-step_min+1), star_M_center(step-step_min+1))
467 0 : xnext = max(x - dx, xmin)
468 : end if
469 0 : pg => pg% next
470 : end do
471 0 : call pgunsa
472 0 : end subroutine plot_burn_data
473 :
474 :
475 0 : subroutine draw_burn_for_step( &
476 : pg, i_burn_type_first, i_burn_type_last, &
477 : bmax, bmin, burn_type_cutoff, xval, mass, mass_center)
478 :
479 : type (pgstar_hist_node), pointer :: pg
480 : integer, intent(in) :: i_burn_type_first, i_burn_type_last
481 : real, intent(in) :: bmax, bmin, burn_type_cutoff, xval, mass, mass_center
482 :
483 0 : real :: burn_qbot, burn_qtop, mbot, mtop, xmass
484 0 : real :: burn_type, color_frac
485 : integer :: k, colormap_index, clr, mid_map
486 : include 'formats'
487 0 : xmass = mass - mass_center
488 0 : burn_qbot = 0
489 0 : mid_map = colormap_length/2
490 0 : do k = i_burn_type_first, i_burn_type_last, 2
491 0 : burn_type = pg% vals(k)
492 0 : burn_qtop = pg% vals(k+1)
493 0 : if (burn_type < -9990) exit
494 0 : if (abs(burn_type) > burn_type_cutoff) then
495 0 : mbot = mass_center + xmass*burn_qbot
496 0 : mtop = mass_center + xmass*burn_qtop
497 0 : if (burn_type > 0.0) then
498 0 : color_frac = 1.0 - max(0.0, min(1.0, burn_type/bmax))
499 : colormap_index = &
500 0 : colormap_length - int(0.6*color_frac*(colormap_length - mid_map))
501 : else ! burn_type < 0.0
502 0 : color_frac = 1.0 - max(0.0, min(1.0, burn_type/bmin))
503 0 : colormap_index = 1 + int(0.6*color_frac*mid_map)
504 : end if
505 0 : clr = colormap_offset + colormap_index
506 0 : call pgsci(clr)
507 0 : call draw1(xval, mbot, mtop, clr)
508 : end if
509 0 : burn_qbot = burn_qtop
510 0 : if (burn_qbot >= 1d0) exit
511 : end do
512 0 : end subroutine draw_burn_for_step
513 :
514 :
515 0 : subroutine plot_L_lines
516 : integer :: cnt, n
517 0 : real :: coords(4), fjusts(4)
518 :
519 : logical, parameter :: dbg = .false.
520 :
521 : include 'formats'
522 :
523 0 : cnt = 0
524 0 : if (have_log_L) cnt = cnt + 1
525 0 : if (have_log_Lneu) cnt = cnt + 1
526 0 : if (have_log_LH) cnt = cnt + 1
527 0 : if (have_log_LHe) cnt = cnt + 1
528 0 : select case(cnt)
529 : case (1)
530 0 : coords(1) = 0.5; fjusts(1) = 0.5
531 : case (2)
532 0 : coords(1) = 0.80; fjusts(1) = 1.0
533 0 : coords(2) = 0.20; fjusts(2) = 0.0
534 : case (3)
535 0 : coords(1) = 0.90; fjusts(1) = 1.0
536 0 : coords(2) = 0.50; fjusts(2) = 0.5
537 0 : coords(3) = 0.10; fjusts(3) = 0.0
538 : case (4)
539 0 : coords(1) = 0.95; fjusts(1) = 1.0
540 0 : coords(2) = 0.75; fjusts(2) = 0.7
541 0 : coords(3) = 0.25; fjusts(3) = 0.5
542 0 : coords(4) = 0.05; fjusts(4) = 0.0
543 : case default
544 0 : return
545 : end select
546 :
547 0 : call pgsave
548 :
549 0 : call setup_L_yaxis
550 :
551 0 : call pgsch(txt_scale*0.8)
552 0 : call pgslw(2)
553 :
554 0 : n = 0
555 :
556 0 : if (have_log_L) then
557 0 : n = n+1
558 0 : call pgsci(clr_Crimson)
559 0 : call show_right_yaxis_label_pgmtxt_pgstar(s,coords(n),fjusts(n),'logL',-1.2)
560 0 : call plot_L_line(log_L)
561 : end if
562 :
563 0 : if (have_log_Lneu) then
564 0 : n = n+1
565 0 : call pgsci(clr_Tan)
566 0 : call show_right_yaxis_label_pgmtxt_pgstar(s,coords(n),fjusts(n),'logL\d\gn',-1.2)
567 0 : call plot_L_line(log_Lneu)
568 : end if
569 :
570 0 : if (have_log_LH) then
571 0 : n = n+1
572 0 : call pgsci(clr_Goldenrod)
573 0 : call show_right_yaxis_label_pgmtxt_pgstar(s,coords(n),fjusts(n),'logL\dHe',-1.2)
574 0 : call plot_L_line(log_LHe)
575 : end if
576 :
577 0 : if (have_log_LHe) then
578 0 : n = n+1
579 0 : call pgsci(clr_Silver)
580 0 : call show_right_yaxis_label_pgmtxt_pgstar(s,coords(n),fjusts(n),'logL\dH',-1.2)
581 0 : call plot_L_line(log_LH)
582 : end if
583 :
584 0 : call pgunsa
585 :
586 : end subroutine plot_L_lines
587 :
588 :
589 0 : subroutine plot_total_mass_line
590 :
591 0 : call pgsave
592 0 : call pgsch(txt_scale*0.8)
593 :
594 0 : call pgsci(clr_Gray)
595 0 : call pgslw(2)
596 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,1.0,1.0,'M\dtotal\u',-0.8)
597 0 : call pgslw(s% pg% pgstar_lw)
598 0 : call plot_mass_line(star_mass)
599 :
600 0 : call pgunsa
601 :
602 0 : end subroutine plot_total_mass_line
603 :
604 :
605 0 : subroutine plot_mass_lines
606 :
607 : include 'formats'
608 :
609 0 : call pgsave
610 0 : call pgsch(txt_scale*0.8)
611 0 : call pgslw(2)
612 :
613 0 : if (have_he_core_mass) then
614 0 : call pgsci(clr_Teal)
615 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.77,0.5,'He',-0.8)
616 0 : call plot_mass_line(he_core_mass)
617 : end if
618 :
619 0 : if (have_c_core_mass) then
620 0 : call pgsci(clr_LightOliveGreen)
621 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.59,0.5,'C',-0.8)
622 0 : call plot_mass_line(c_core_mass)
623 : end if
624 :
625 0 : if (have_o_core_mass) then
626 0 : call pgsci(clr_SeaGreen)
627 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.41,0.5,'O',-0.8)
628 0 : call plot_mass_line(o_core_mass)
629 : end if
630 :
631 0 : if (have_si_core_mass) then
632 0 : call pgsci(clr_Lilac)
633 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.23,0.5,'Si',-0.8)
634 0 : call plot_mass_line(si_core_mass)
635 : end if
636 :
637 0 : if (have_fe_core_mass) then
638 0 : call pgsci(clr_Crimson)
639 0 : call show_left_yaxis_label_pgmtxt_pgstar(s,0.00,0.0,'Iron',-0.8)
640 0 : call plot_mass_line(fe_core_mass)
641 : end if
642 :
643 0 : call pgunsa
644 :
645 0 : end subroutine plot_mass_lines
646 :
647 :
648 0 : subroutine show_mix_legend
649 : integer :: mix_type
650 :
651 0 : call pgsave
652 0 : call pgslw(2)
653 :
654 0 : mix_type = convective_mixing
655 0 : call pgsci(mix_clr(mix_type))
656 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.05, 0.0, 'conv', 0.0)
657 :
658 0 : mix_type = leftover_convective_mixing
659 0 : call pgsci(mix_clr(mix_type))
660 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.275, 0.5, 'left', 0.0)
661 :
662 0 : mix_type = overshoot_mixing
663 0 : call pgsci(mix_clr(mix_type))
664 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.5, 0.5, 'over', 0.0)
665 :
666 0 : mix_type = semiconvective_mixing
667 0 : call pgsci(mix_clr(mix_type))
668 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.725, 0.5, 'semi', 0.0)
669 :
670 0 : mix_type = thermohaline_mixing
671 0 : call pgsci(mix_clr(mix_type))
672 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.95, 1.0, 'thrm', 0.0)
673 :
674 0 : mix_type = rayleigh_taylor_mixing
675 0 : call pgsci(mix_clr(mix_type))
676 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.85, 0.75, 'RTI', 0.0)
677 :
678 0 : call pgunsa
679 :
680 0 : end subroutine show_mix_legend
681 :
682 :
683 0 : subroutine show_burn_legend
684 : integer :: colormap_index
685 :
686 0 : call pgsave
687 0 : call pgslw(2)
688 :
689 0 : colormap_index = int(colormap_length*0.85)
690 0 : call pgsci(colormap_offset + colormap_index)
691 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.17, 0.5, 'burning', 1.0)
692 :
693 0 : colormap_index = int(colormap_length*0.15)
694 0 : call pgsci(colormap_offset + colormap_index)
695 0 : call show_xaxis_label_pgmtxt_pgstar(s, 0.82, 0.5, 'cooling', 1.0)
696 :
697 0 : call pgunsa
698 :
699 0 : end subroutine show_burn_legend
700 :
701 :
702 0 : subroutine plot_mass_line(yvec)
703 : real, intent(in) :: yvec(:)
704 0 : if (any(yvec > 1e-2*ymax_mass_axis)) then
705 0 : call pgsave
706 0 : call pgslw(s% pg% Kipp_masses_line_weight)
707 0 : call pgline(n, xvec, yvec)
708 0 : call pgunsa
709 : end if
710 0 : end subroutine plot_mass_line
711 :
712 :
713 0 : subroutine plot_L_line(yvec)
714 : real, intent(in) :: yvec(:)
715 0 : if (any(yvec > ymin_L_axis)) then
716 0 : call pgsave
717 0 : call pgslw(s% pg% Kipp_luminosities_line_weight)
718 0 : call pgline(n, xvec, yvec)
719 0 : call pgunsa
720 : end if
721 0 : end subroutine plot_L_line
722 :
723 :
724 0 : subroutine plot_mix_data
725 : use history_specs, only: mixing_offset
726 : type (pgstar_hist_node), pointer :: pg
727 : integer :: i_mix_type_first, i_mix_type_last
728 : integer :: k, cnt, num_specs, step
729 :
730 : include 'formats'
731 :
732 0 : i_mix_type_first = 0
733 0 : num_specs = size(s% history_column_spec, dim=1)
734 0 : do k = 1, num_specs
735 0 : if (s% history_column_spec(k) == mixing_offset + 1) then
736 0 : i_mix_type_first = k
737 0 : exit
738 : end if
739 : end do
740 0 : if (i_mix_type_first == 0) then
741 : !write(*,*) 'i_mix_type_first == 0'
742 : return
743 : end if
744 :
745 0 : i_mix_type_last = 0
746 0 : cnt = 1
747 0 : do k=i_mix_type_first+1, num_specs
748 0 : i_mix_type_last = k-1
749 0 : cnt = cnt+1
750 0 : if (s% history_column_spec(k) /= mixing_offset + cnt) exit
751 : end do
752 :
753 0 : call pgsave
754 0 : call pgslw(s% pg% Kipp_mix_line_weight)
755 0 : if (.not. associated(s% pg% pgstar_hist)) then
756 0 : write(*,*) '.not. associated(s% pg% pgstar_hist)'
757 : end if
758 0 : pg => s% pg% pgstar_hist
759 0 : do
760 0 : if (.not. associated(pg)) exit
761 0 : step = pg% step
762 0 : if (step < step_min) exit
763 0 : if (step <= step_max .and. mod(step, s% pg% Kipp_mix_interval) == 0) then
764 : call draw_mix_for_step( &
765 : pg, step, i_mix_type_first, i_mix_type_last, real(xvec(step-step_min+1)), &
766 0 : star_mass(step-step_min+1), star_M_center(step-step_min+1))
767 : end if
768 0 : pg => pg% next
769 : end do
770 0 : call pgunsa
771 0 : end subroutine plot_mix_data
772 :
773 :
774 0 : subroutine draw_mix_for_step( &
775 : pg, step, i_mix_type_first, i_mix_type_last, xval, mass, mass_center)
776 : type (pgstar_hist_node), pointer :: pg
777 : integer, intent(in) :: step, i_mix_type_first, i_mix_type_last
778 : real, intent(in) :: xval, mass, mass_center
779 0 : real :: mix_qbot, mix_qtop, mbot, mtop, xmass
780 : integer :: k, mix_type
781 : include 'formats'
782 0 : mix_qbot = 0
783 0 : xmass = mass - mass_center
784 0 : do k = i_mix_type_first, i_mix_type_last, 2
785 0 : mix_type = int(pg% vals(k))
786 0 : if (mix_type < 0) exit
787 0 : mix_qtop = pg% vals(k+1)
788 0 : if (mix_type > 0 .and. mix_type <= max_mix_type_to_show) then
789 0 : mbot = mass_center + xmass*mix_qbot
790 0 : mtop = mass_center + xmass*mix_qtop
791 0 : call draw1(xval, mbot, mtop, mix_clr(mix_type))
792 0 : showed_this_mix_type(mix_type) = .true.
793 : end if
794 0 : mix_qbot = mix_qtop
795 : end do
796 0 : end subroutine draw_mix_for_step
797 :
798 :
799 0 : subroutine draw1(xval,y1,y2,clr)
800 : real, intent(in) :: xval, y1, y2
801 : integer, intent(in) :: clr
802 0 : real :: top, bot
803 0 : if (y1 < y2) then
804 0 : bot = y1; top = y2
805 : else
806 0 : bot = y2; top = y1
807 : end if
808 0 : if (top < 1d-50) return
809 0 : call pgsci(clr)
810 0 : call pgmove(xval, bot)
811 0 : call pgdraw(xval, top)
812 : end subroutine draw1
813 :
814 : end subroutine do_Kipp_Plot
815 :
816 : end module pgstar_kipp
|
