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_support
21 :
22 : use star_private_def
23 : use const_def, only: dp, secday, secyer, mesa_data_dir, &
24 : overshoot_mixing, rotation_mixing, thermohaline_mixing, semiconvective_mixing, &
25 : leftover_convective_mixing, convective_mixing, no_mixing
26 : use rates_def, only : i_rate
27 : use utils_lib
28 : use star_pgstar
29 :
30 : implicit none
31 :
32 : integer, parameter :: category_offset = 1000
33 : integer, parameter :: abundance_offset = 2000
34 : integer, parameter :: extras_offset = 3000
35 :
36 : logical :: have_initialized_pgstar = .false.
37 :
38 : real :: sum_dHR_since_last_file_write
39 :
40 :
41 : ! lines for TRho profile
42 : real, dimension(:), allocatable :: hydrogen_burn_logT, hydrogen_burn_logRho
43 : real, dimension(:), allocatable :: helium_burn_logT, helium_burn_logRho
44 : real, dimension(:), allocatable :: carbon_burn_logT, carbon_burn_logRho
45 : real, dimension(:), allocatable :: oxygen_burn_logT, oxygen_burn_logRho
46 : real, dimension(:), allocatable :: psi4_logT, psi4_logRho
47 : real, dimension(:), allocatable :: elect_data_logT, elect_data_logRho
48 : real, dimension(:), allocatable :: gamma_4_thirds_logT, gamma_4_thirds_logRho
49 : real, dimension(:), allocatable :: kap_rad_cond_eq_logT, kap_rad_cond_eq_logRho
50 : real, dimension(:), allocatable :: opal_clip_logT, opal_clip_logRho
51 : real, dimension(:), allocatable :: scvh_clip_logT, scvh_clip_logRho
52 :
53 : ! Tioga line types
54 : integer, parameter :: Line_Type_Solid = 1
55 : integer, parameter :: Line_Type_Dash = 2
56 : integer, parameter :: Line_Type_Dot_Dash = 3
57 : integer, parameter :: Line_Type_Dot = 4
58 :
59 : contains
60 :
61 :
62 0 : subroutine add_to_pgstar_hist(s, pg_hist_new)
63 : type (star_info), pointer :: s
64 : type (pgstar_hist_node), pointer :: pg_hist_new
65 : type (pgstar_hist_node), pointer :: next => null()
66 : integer :: step
67 0 : step = pg_hist_new% step
68 0 : do
69 0 : if (.not. associated(s% pg% pgstar_hist)) then
70 0 : s% pg% pgstar_hist => pg_hist_new
71 0 : nullify(pg_hist_new% next)
72 0 : return
73 : end if
74 0 : if (step > s% pg% pgstar_hist% step) then
75 0 : pg_hist_new% next => s% pg% pgstar_hist
76 0 : s% pg% pgstar_hist => pg_hist_new
77 0 : return
78 : end if
79 : ! discard item
80 0 : next => s% pg% pgstar_hist% next
81 0 : deallocate(s% pg% pgstar_hist% vals)
82 0 : deallocate(s% pg% pgstar_hist)
83 0 : s% pg% pgstar_hist => next
84 : end do
85 : end subroutine add_to_pgstar_hist
86 :
87 :
88 2 : subroutine pgstar_clear(s)
89 : type (star_info), pointer :: s
90 : integer :: i
91 : type (pgstar_win_file_data), pointer :: p
92 : type (pgstar_hist_node), pointer :: pg_hist => null(), next => null()
93 2 : pg_hist => s% pg% pgstar_hist
94 2 : do while(associated(pg_hist))
95 0 : if (associated(pg_hist% vals)) deallocate(pg_hist% vals)
96 0 : next => pg_hist% next
97 0 : deallocate(pg_hist)
98 0 : pg_hist => next
99 : end do
100 2 : nullify(s% pg% pgstar_hist)
101 2 : if (have_initialized_pgstar) return
102 196 : do i = 1, num_pgstar_plots
103 194 : p => s% pg% pgstar_win_file_ptr(i)
104 194 : p% id_win = 0
105 194 : p% have_called_mkdir = .false.
106 196 : p% file_dir_for_previous_mkdir = ''
107 : end do
108 : end subroutine pgstar_clear
109 :
110 :
111 0 : subroutine init_pgstar(ierr)
112 : integer, intent(out) :: ierr
113 :
114 0 : if (have_initialized_pgstar) return
115 :
116 0 : call read_support_info(ierr)
117 0 : if (failed('read_support_info')) return
118 :
119 0 : have_initialized_pgstar = .true.
120 0 : sum_dHR_since_last_file_write = 0
121 :
122 : contains
123 :
124 0 : logical function failed(str)
125 : character (len = *), intent(in) :: str
126 0 : failed = (ierr /= 0)
127 0 : if (failed) then
128 0 : write(*, *) trim(str) // ' ierr', ierr
129 : end if
130 0 : end function failed
131 :
132 : end subroutine init_pgstar
133 :
134 0 : subroutine check_window(s, p, ierr)
135 : type (star_info), pointer :: s
136 : type (pgstar_win_file_data), pointer :: p
137 : integer, intent(out) :: ierr
138 0 : ierr = 0
139 0 : if (p% do_win .and. (.not. p% win_flag)) then
140 0 : p% do_win = .false.
141 0 : if (p% id_win > 0) then
142 0 : call pgslct(p% id_win)
143 0 : call pgclos
144 0 : p% id_win = 0
145 : end if
146 0 : else if (p% win_flag .and. (.not. p% do_win)) then
147 0 : if (p% id_win == 0) &
148 0 : call open_device(s, p, .false., '/xwin', p% id_win, ierr)
149 0 : if (ierr == 0 .and. p% id_win > 0) p% do_win = .true.
150 : end if
151 0 : if (p% do_win .and. p% id_win > 0 .and. &
152 : (p% win_width /= p% prev_win_width .or. &
153 : p% win_aspect_ratio /= p% prev_win_aspect_ratio)) then
154 0 : call pgslct(p% id_win)
155 0 : call pgpap(p% win_width, p% win_aspect_ratio)
156 0 : p% prev_win_width = p% win_width
157 0 : p% prev_win_aspect_ratio = p% win_aspect_ratio
158 : end if
159 0 : end subroutine check_window
160 :
161 :
162 0 : subroutine check_file(s, p, ierr)
163 : use utils_lib, only : mkdir
164 : type (star_info), pointer :: s
165 : type (pgstar_win_file_data), pointer :: p
166 : integer, intent(out) :: ierr
167 : character (len = strlen) :: name
168 0 : ierr = 0
169 0 : if (p% do_file .and. (.not. p% file_flag)) then
170 0 : p% do_file = .false.
171 0 : else if (p% file_flag .and. (.not. p% do_file)) then
172 0 : if (p% id_file == 0) then
173 0 : if (.not. p% have_called_mkdir .or. &
174 : p% file_dir /= p% file_dir_for_previous_mkdir) then
175 0 : if(.not. folder_exists(trim(p% file_dir))) call mkdir(trim(p% file_dir))
176 0 : p% have_called_mkdir = .true.
177 0 : p% file_dir_for_previous_mkdir = p% file_dir
178 : end if
179 0 : call create_file_name(s, p% file_dir, p% file_prefix, name)
180 0 : name = trim(name) // '/' // trim(s% pg% file_device)
181 0 : call open_device(s, p, .true., name, p% id_file, ierr)
182 0 : if (ierr /= 0) return
183 0 : p% most_recent_filename = name
184 : end if
185 0 : p% do_file = .true.
186 : end if
187 : end subroutine check_file
188 :
189 :
190 0 : subroutine create_file_name(s, dir, prefix, name)
191 : use star_utils, only : get_string_for_model_number
192 : type (star_info), pointer :: s
193 : character (len = *), intent(in) :: dir, prefix
194 : character (len = *), intent(out) :: name
195 : character (len = strlen) :: num_str, fstring
196 : character (len = 32) :: file_extension
197 0 : write(fstring, '( "(i",i2.2,".",i2.2,")" )') s% pg% file_digits, s% pg% file_digits
198 0 : write(num_str, fstring) s% model_number
199 0 : if (len_trim(dir) > 0) then
200 0 : name = trim(dir) // '/' // trim(prefix)
201 : else
202 0 : name = prefix
203 : end if
204 0 : if (s%pg%file_device=='vcps') then
205 0 : file_extension = 'ps'
206 : else
207 0 : file_extension = s%pg%file_device ! e.g.: png, ps
208 : end if
209 0 : name = trim(name) // trim(num_str) // '.' // trim(file_extension)
210 0 : end subroutine create_file_name
211 :
212 :
213 0 : subroutine write_plot_to_file(s, p, filename, ierr)
214 : type (star_info), pointer :: s
215 : type (pgstar_win_file_data), pointer :: p
216 : character (len = *), intent(in) :: filename
217 : integer, intent(out) :: ierr
218 : character (len = strlen) :: name
219 0 : ierr = 0
220 : !name = trim(filename) // '/' // trim(s% pg% file_device)
221 0 : name = trim(filename) // '/png'
222 0 : write(*, '(a)') 'write_plot_to_file device: ' // trim(name)
223 0 : call open_device(s, p, .true., trim(name), p% id_file, ierr)
224 0 : if (ierr /= 0) then
225 0 : write(*, *) 'failed in open_device'
226 0 : return
227 : end if
228 0 : call p% plot(s% id, p% id_file, ierr)
229 0 : call pgclos
230 0 : p% id_file = 0
231 0 : p% do_file = .false.
232 : end subroutine write_plot_to_file
233 :
234 :
235 0 : subroutine open_device(s, p, is_file, dev, id, ierr)
236 : use pgstar_colors, only: set_device_colors
237 : type (star_info), pointer :: s
238 : type (pgstar_win_file_data), pointer :: p
239 : logical, intent(in) :: is_file
240 : character (len = *), intent(in) :: dev
241 : integer, intent(out) :: id
242 : integer, intent(out) :: ierr
243 :
244 : integer :: pgopen
245 : character (len = strlen) :: dir
246 : logical :: white_on_black_flag
247 : real :: width, ratio
248 :
249 0 : if (is_file) then
250 : dir = p% file_dir
251 0 : white_on_black_flag = s% pg% file_white_on_black_flag
252 : else
253 : dir = ''
254 0 : white_on_black_flag = s% pg% win_white_on_black_flag
255 : end if
256 :
257 0 : ierr = 0
258 0 : id = -1
259 0 : id = pgopen(trim(dev))
260 0 : if (id <= 0) return
261 :
262 : !write(*,*) 'open device <' // trim(dev) // '> ' // trim(p% name), id
263 0 : if (is_file) then
264 0 : width = p% file_width; if (width < 0) width = p% win_width
265 0 : ratio = p% file_aspect_ratio; if (ratio < 0) ratio = p% win_aspect_ratio
266 0 : call pgpap(width, ratio)
267 : else
268 0 : call pgpap(p% win_width, p% win_aspect_ratio)
269 0 : p% prev_win_width = p% win_width
270 0 : p% prev_win_aspect_ratio = p% win_aspect_ratio
271 : end if
272 0 : call set_device_colors(white_on_black_flag)
273 : end subroutine open_device
274 :
275 :
276 0 : subroutine read_support_info(ierr)
277 : integer, intent(out) :: ierr
278 :
279 : ierr = 0
280 :
281 : call read_TRho_data(&
282 0 : 'hydrogen_burn.data', hydrogen_burn_logT, hydrogen_burn_logRho, ierr)
283 0 : if (ierr /= 0) then
284 0 : write(*, *) 'PGSTAR failed in reading hydrogen burn data'
285 0 : return
286 : end if
287 :
288 : call read_TRho_data(&
289 0 : 'helium_burn.data', helium_burn_logT, helium_burn_logRho, ierr)
290 0 : if (ierr /= 0) then
291 0 : write(*, *) 'PGSTAR failed in reading helium burn data'
292 0 : return
293 : end if
294 :
295 : call read_TRho_data(&
296 0 : 'carbon_burn.data', carbon_burn_logT, carbon_burn_logRho, ierr)
297 0 : if (ierr /= 0) then
298 0 : write(*, *) 'PGSTAR failed in reading carbon burn data'
299 0 : return
300 : end if
301 :
302 : call read_TRho_data(&
303 0 : 'oxygen_burn.data', oxygen_burn_logT, oxygen_burn_logRho, ierr)
304 0 : if (ierr /= 0) then
305 0 : write(*, *) 'PGSTAR failed in reading oxygen burn data'
306 0 : return
307 : end if
308 :
309 : call read_TRho_data(&
310 0 : 'psi4.data', psi4_logT, psi4_logRho, ierr)
311 0 : if (ierr /= 0) then
312 0 : write(*, *) 'PGSTAR failed in reading psi4 data'
313 0 : return
314 : end if
315 :
316 : call read_TRho_data(&
317 0 : 'elect.data', elect_data_logT, elect_data_logRho, ierr)
318 0 : if (ierr /= 0) then
319 0 : write(*, *) 'PGSTAR failed in reading elect data'
320 0 : return
321 : end if
322 :
323 : call read_TRho_data(&
324 0 : 'gamma_4_thirds.data', gamma_4_thirds_logT, gamma_4_thirds_logRho, ierr)
325 0 : if (ierr /= 0) then
326 0 : write(*, *) 'PGSTAR failed in reading gamma_4_thirds data'
327 0 : return
328 : end if
329 :
330 : call read_TRho_data(&
331 0 : 'kap_rad_cond_eq.data', kap_rad_cond_eq_logT, kap_rad_cond_eq_logRho, ierr)
332 0 : if (ierr /= 0) then
333 0 : write(*, *) 'PGSTAR failed in reading kap_rad_cond_eq data'
334 0 : return
335 : end if
336 :
337 : call read_TRho_data(&
338 0 : 'opal_clip.data', opal_clip_logT, opal_clip_logRho, ierr)
339 0 : if (ierr /= 0) then
340 0 : write(*, *) 'PGSTAR failed in reading opal_clip data'
341 0 : return
342 : end if
343 :
344 : call read_TRho_data(&
345 0 : 'scvh_clip.data', scvh_clip_logT, scvh_clip_logRho, ierr)
346 0 : if (ierr /= 0) then
347 0 : write(*, *) 'PGSTAR failed in reading scvh_clip data'
348 0 : return
349 : end if
350 :
351 : end subroutine read_support_info
352 :
353 0 : subroutine read_TRho_data(fname, logTs, logRhos, ierr)
354 : use utils_lib
355 : character (len = *), intent(in) :: fname
356 : real, dimension(:), allocatable :: logTs, logRhos ! will allocate
357 : integer, intent(out) :: ierr
358 :
359 : character (len = strlen) :: filename
360 : real :: logT, logRho
361 : integer :: iounit, i, sz, cnt
362 :
363 0 : filename = trim(mesa_data_dir) // '/star_data/plot_info/' // trim(fname)
364 :
365 0 : open(newunit = iounit, file = trim(filename), status = 'old', action = 'read', iostat = ierr)
366 0 : if (ierr/=0) then
367 0 : write(*, *) 'failed to open ' // trim(filename)
368 0 : call done
369 0 : return
370 : end if
371 :
372 : sz = 0
373 0 : do
374 0 : read(iounit, *, iostat = ierr)
375 0 : if(ierr/=0) exit
376 0 : sz = sz + 1
377 : end do
378 :
379 0 : rewind(iounit)
380 :
381 0 : allocate(logTs(sz), logRhos(sz))
382 :
383 0 : cnt = 0
384 0 : do i = 1, sz
385 0 : read(iounit, *, iostat = ierr) logRho, logT
386 0 : if (ierr /= 0) then
387 0 : ierr = 0; exit
388 : end if
389 0 : logRhos(i) = logRho
390 0 : logTs(i) = logT
391 0 : cnt = i
392 : end do
393 :
394 0 : call done
395 :
396 :
397 : contains
398 :
399 :
400 0 : subroutine done
401 0 : close(iounit)
402 0 : end subroutine done
403 :
404 : end subroutine read_TRho_data
405 :
406 :
407 0 : integer function write_info_line_str(cnt, ypos, xpos0, dxpos, str)
408 : integer, intent(in) :: cnt
409 : real, intent(in) :: ypos, xpos0, dxpos
410 : character (len = *), intent(in) :: str
411 : real :: xpos
412 0 : xpos = cnt * dxpos + xpos0
413 0 : call pgptxt(xpos, ypos, 0.0, 0.5, trim(adjustl(str)))
414 0 : write_info_line_str = cnt + 1
415 0 : end function write_info_line_str
416 :
417 :
418 0 : integer function write_info_line_int(cnt, ypos, xpos0, dxpos, dxval, label, val)
419 : integer, intent(in) :: cnt, val
420 : real, intent(in) :: ypos, xpos0, dxpos, dxval
421 : character (len = *), intent(in) :: label
422 :
423 : character (len = 128) :: str
424 : real :: xpos
425 :
426 0 : write(str, '(a)') trim(label)
427 0 : xpos = cnt * dxpos + xpos0
428 0 : call pgptxt(xpos, ypos, 0.0, 1.0, trim(adjustl(str)))
429 0 : write(str, '(i9)') val
430 0 : xpos = xpos + dxval
431 0 : call pgptxt(xpos, ypos, 0.0, 0.0, trim(adjustl(str)))
432 :
433 0 : write_info_line_int = cnt + 1
434 0 : end function write_info_line_int
435 :
436 :
437 0 : integer function write_info_line_flt(&
438 : cnt, ypos, xpos0, dxpos, dxval, label, val)
439 : integer, intent(in) :: cnt
440 : real, intent(in) :: ypos, xpos0, dxpos, dxval
441 : real(dp), intent(in) :: val
442 : character (len = *), intent(in) :: label
443 :
444 : character (len = 128) :: str
445 : real :: xpos
446 : integer :: ierr
447 :
448 0 : write_info_line_flt = cnt + 1
449 0 : write(str, '(a)') trim(label)
450 0 : xpos = cnt * dxpos + xpos0
451 0 : call pgptxt(xpos, ypos, 0.0, 1.0, trim(adjustl(str)))
452 0 : ierr = 0
453 0 : write(str, '(f12.7)', iostat = ierr) val
454 0 : if (ierr /= 0) then
455 0 : ierr = 0
456 0 : write(str, '(e10.3)', iostat = ierr) val
457 0 : if (ierr /= 0) then
458 0 : write(*, *) trim(label), val
459 0 : write(*, *) 'problem in write_info_line_flt'
460 0 : return
461 : end if
462 : end if
463 0 : xpos = xpos + dxval
464 0 : call pgptxt(xpos, ypos, 0.0, 0.0, trim(adjustl(str)))
465 :
466 0 : end function write_info_line_flt
467 :
468 :
469 0 : integer function write_info_line_flt2(cnt, ypos, xpos0, dxpos, dxval, label, val)
470 : integer, intent(in) :: cnt
471 : real, intent(in) :: ypos, xpos0, dxpos, dxval
472 : real(dp), intent(in) :: val
473 : character (len = *), intent(in) :: label
474 :
475 : character (len = 128) :: str
476 : real :: xpos
477 : integer :: ierr
478 :
479 0 : write_info_line_flt2 = cnt + 1
480 :
481 0 : write(str, '(a)') trim(label)
482 0 : xpos = cnt * dxpos + xpos0
483 0 : call pgptxt(xpos, ypos, 0.0, 1.0, trim(adjustl(str)))
484 0 : ierr = 0
485 0 : write(str, '(f12.3)', iostat = ierr) val
486 0 : if (ierr /= 0) then
487 0 : ierr = 0
488 0 : write(str, '(e10.3)', iostat = ierr) val
489 0 : if (ierr /= 0) then
490 0 : write(*, *) trim(label), val
491 0 : write(*, *) 'problem in write_info_line_flt2'
492 0 : return
493 : end if
494 : end if
495 0 : xpos = xpos + dxval
496 0 : call pgptxt(xpos, ypos, 0.0, 0.0, trim(adjustl(str)))
497 :
498 0 : end function write_info_line_flt2
499 :
500 :
501 0 : integer function write_info_line_exp(cnt, ypos, xpos0, dxpos, dxval, label, val)
502 : integer, intent(in) :: cnt
503 : real, intent(in) :: ypos, xpos0, dxpos, dxval
504 : real(dp), intent(in) :: val
505 : character (len = *), intent(in) :: label
506 :
507 : character (len = 128) :: str
508 : real :: xpos
509 :
510 0 : write(str, '(a)') trim(label)
511 0 : xpos = cnt * dxpos + xpos0
512 0 : call pgptxt(xpos, ypos, 0.0, 1.0, trim(adjustl(str)))
513 0 : write(str, '(1pe10.3)') val
514 0 : xpos = xpos + dxval
515 0 : call pgptxt(xpos, ypos, 0.0, 0.0, trim(adjustl(str)))
516 :
517 0 : write_info_line_exp = cnt + 1
518 0 : end function write_info_line_exp
519 :
520 :
521 0 : integer function count_hist_points(&
522 : s, step_min, step_max) result(numpts)
523 : type (star_info), pointer :: s
524 : integer, intent(in) :: step_min, step_max
525 : type (pgstar_hist_node), pointer :: pg
526 : include 'formats'
527 0 : numpts = 0
528 0 : pg => s% pg% pgstar_hist
529 0 : do ! recall that hist list is decreasing by age (and step)
530 0 : if (.not. associated(pg)) return
531 0 : if (pg% step < step_min) return
532 0 : if (pg% step <= step_max .or. step_max <= 0) numpts = numpts + 1
533 0 : pg => pg% next
534 : end do
535 : end function count_hist_points
536 :
537 :
538 0 : logical function get1_hist_yvec(s, step_min, step_max, n, name, vec)
539 : use utils_lib, only : integer_dict_lookup
540 : type (star_info), pointer :: s
541 : integer, intent(in) :: step_min, step_max, n ! n = count_hist_points
542 : character (len = *) :: name
543 : real, dimension(:), allocatable :: vec
544 : integer :: i, cnt, ierr
545 : character (len = 64) :: key_name
546 : include 'formats'
547 0 : cnt = 0
548 : ierr = 0
549 0 : do i = 1, len(key_name)
550 0 : key_name(i:i) = ' '
551 : end do
552 0 : do i = 1, len_trim(name)
553 0 : if (name(i:i) == ' ') then
554 0 : cnt = cnt + 1
555 0 : key_name(i:i) = '_'
556 : else
557 0 : key_name(i:i) = name(i:i)
558 : end if
559 : end do
560 0 : call integer_dict_lookup(s% history_names_dict, key_name, i, ierr)
561 0 : if (ierr /= 0 .or. i <= 0) then ! didn't find it
562 : get1_hist_yvec = .false.
563 : return
564 : end if
565 0 : call get_hist_points(s, step_min, step_max, n, i, vec, ierr)
566 0 : if (ierr /= 0) then ! didn't get them
567 : get1_hist_yvec = .false.
568 : return
569 : end if
570 : get1_hist_yvec = .true.
571 : end function get1_hist_yvec
572 :
573 :
574 0 : subroutine set_hist_points_steps(&
575 0 : s, step_min, step_max, numpts, vec, ierr)
576 : type (star_info), pointer :: s
577 : integer, intent(in) :: step_min, step_max, numpts
578 : real, intent(out) :: vec(:)
579 : integer, intent(out) :: ierr
580 : integer :: i
581 : type (pgstar_hist_node), pointer :: pg
582 0 : ierr = 0
583 0 : if (numpts == 0) return
584 0 : pg => s% pg% pgstar_hist
585 0 : i = numpts
586 0 : do ! recall that hist list is decreasing by age (and step)
587 0 : if (.not. associated(pg)) then
588 0 : ierr = -1
589 0 : return
590 : end if
591 0 : if (pg% step < step_min) then
592 0 : ierr = -1
593 0 : return
594 : end if
595 0 : if (pg% step <= step_max) then
596 0 : vec(i) = real(pg% step)
597 0 : i = i - 1
598 0 : if (i == 0) return
599 : end if
600 0 : pg => pg% next
601 : end do
602 : end subroutine set_hist_points_steps
603 :
604 :
605 0 : integer function get_hist_index(s, spec) result(index)
606 : type (star_info), pointer :: s
607 : integer, intent(in) :: spec
608 : integer :: i, num
609 : ! note: this doesn't include "extra" columns
610 0 : num = size(s% history_column_spec, dim = 1)
611 0 : do i = 1, num
612 0 : if (s% history_column_spec(i) == spec) then
613 : index = i
614 : return
615 : end if
616 : end do
617 : index = -1
618 : end function get_hist_index
619 :
620 :
621 0 : subroutine get_hist_points(&
622 0 : s, step_min, step_max, numpts, index, vec, ierr)
623 : type (star_info), pointer :: s
624 : integer, intent(in) :: step_min, step_max, numpts, index
625 : real, intent(out) :: vec(:)
626 : integer, intent(out) :: ierr
627 : integer :: i
628 : type (pgstar_hist_node), pointer :: pg => null()
629 : include 'formats'
630 0 : if (numpts == 0) return
631 0 : pg => s% pg% pgstar_hist
632 0 : i = numpts
633 0 : vec = 0
634 0 : ierr = 0
635 0 : do ! recall that hist list is decreasing by age (and step)
636 0 : if (.not. associated(pg)) return
637 0 : if (pg% step < step_min) then
638 : ! this will not happen if have correct numpts
639 : return
640 : end if
641 0 : if (pg% step <= step_max .or. step_max <= 0) then
642 0 : if (.not. associated(pg% vals)) then
643 : !ierr = -1
644 : !write(*,6) 'failed in get_hist_points: not associated', &
645 : ! s% model_number, index, numpts, step_min, step_max
646 : !call mesa_error(__FILE__,__LINE__,'get_hist_points')
647 : return
648 : end if
649 0 : if (size(pg% vals, dim = 1) < index) then
650 : !ierr = -1
651 : !write(*,7) 'failed in get_hist_points: size < index', &
652 : ! s% model_number, size(pg% vals,dim=1), index, numpts, step_min, step_max
653 : !call mesa_error(__FILE__,__LINE__,'get_hist_points')
654 : return
655 : end if
656 0 : vec(i) = pg% vals(index)
657 0 : i = i - 1
658 0 : if (i == 0) return
659 : end if
660 0 : pg => pg% next
661 : end do
662 : end subroutine get_hist_points
663 :
664 :
665 0 : logical function find_shock(s, xaxis_id, xshock) result(found_shock)
666 : use profile, only : get_profile_val
667 : use num_lib, only : find0
668 : type (star_info), pointer :: s
669 : integer, intent(in) :: xaxis_id
670 : real(dp), intent(out) :: xshock
671 : integer :: k, nz
672 : real(dp) :: cs, x00, xp1, ms
673 : real(dp), pointer :: v(:) => null()
674 : include 'formats'
675 0 : nz = s% nz
676 0 : if (s% u_flag) then
677 0 : v => s% u
678 : else
679 0 : v => s% v
680 : end if
681 : ! search in from surface
682 0 : do k = 1, nz - 1
683 0 : cs = s% csound(k) ! cell center
684 0 : if (v(k + 1) >= cs .and. v(k) < cs) then
685 : found_shock = .true.
686 : exit
687 : end if
688 : end do
689 : if (.not. found_shock) then
690 : ! search out from center
691 : do k = nz - 1, 1, -1
692 : cs = s% csound(k) ! cell center
693 : if (v(k + 1) >= -cs .and. v(k) < -cs) then
694 : found_shock = .true.
695 : exit
696 : end if
697 : end do
698 : end if
699 : if (found_shock) then
700 0 : x00 = get_profile_val(s, xaxis_id, k)
701 0 : xp1 = get_profile_val(s, xaxis_id, k + 1)
702 0 : cs = s% csound(k)
703 0 : ms = find0(0d0, s% dm(k), v(k + 1) - cs, v(k) - cs)
704 0 : xshock = xp1 + (x00 - xp1) * ms / s% dm(k)
705 0 : if (is_bad(xshock)) then
706 0 : write(*, 2) 'xshock nz', nz, xshock
707 0 : write(*, 2) 's% dm(k)', k, s% dm(k)
708 0 : write(*, 2) 's% csound(k)', k, s% csound(k)
709 0 : write(*, 2) 'v(k)', k, v(k)
710 0 : write(*, 2) 'v(k+1)', k + 1, v(k + 1)
711 0 : write(*, 1) 'ms', ms
712 0 : write(*, 1) 'x00', x00
713 0 : write(*, 1) 'xp1', xp1
714 0 : nullify(v)
715 0 : call mesa_error(__FILE__, __LINE__, 'find_shock')
716 : end if
717 : end if
718 0 : nullify(v)
719 0 : end function find_shock
720 :
721 :
722 0 : subroutine set_grid_minmax(&
723 0 : nz, xvec, xmin, xmax, xleft, xright, xaxis_by, &
724 : given_xmin, given_xmax, margin_in, reversed, &
725 : grid_min, grid_max, dxmin)
726 : integer, intent(in) :: nz
727 : real, intent(in), dimension(:) :: xvec
728 : real, intent(out) :: xmin, xmax, xleft, xright
729 : character (len = *), intent(in) :: xaxis_by
730 : real, intent(in) :: given_xmin, given_xmax, margin_in
731 : real, intent(in) :: dxmin
732 : logical, intent(in) :: reversed
733 : integer, intent(out) :: grid_min, grid_max
734 : integer :: k
735 : real :: dx, margin
736 : logical :: use_given_xmin, use_given_xmax
737 :
738 : include 'formats'
739 :
740 0 : margin = max(0.0, margin_in)
741 :
742 : ! use given if it isn't = 101
743 0 : use_given_xmin = abs(given_xmin + 101.0) > 1e-6
744 0 : if (xaxis_by == 'mass' .and. given_xmin < 0 .and. use_given_xmin) then
745 0 : xmin = maxval(xvec(1:nz), mask=.not. is_bad(xvec(1:nz))) + given_xmin
746 0 : else if (use_given_xmin) then
747 0 : xmin = given_xmin
748 0 : else if (xaxis_by == 'logxm' .or. xaxis_by == 'logxq') then
749 0 : xmin = minval(xvec(2:nz), mask=.not. is_bad(xvec(2:nz)))
750 : else
751 0 : xmin = minval(xvec(1:nz), mask=.not. is_bad(xvec(1:nz)))
752 : end if
753 :
754 0 : use_given_xmax = abs(given_xmax + 101.0) > 1e-6
755 0 : if (xaxis_by == 'mass' .and. given_xmax < 0 .and. use_given_xmax) then
756 0 : xmax = maxval(xvec(1:nz), mask=.not. is_bad(xvec(1:nz))) + given_xmax
757 0 : else if (use_given_xmax) then
758 0 : xmax = given_xmax
759 : else
760 0 : xmax = maxval(xvec(1:nz), mask=.not. is_bad(xvec(1:nz)))
761 : end if
762 0 : dx = xmax - xmin
763 :
764 0 : if (is_bad(dx)) then
765 0 : xmax = given_xmax
766 0 : xmin = given_xmin
767 0 : dx = xmax - xmin
768 : end if
769 :
770 0 : if (.not. use_given_xmin) xmin = xmin - margin * dx
771 0 : if (.not. use_given_xmax) xmax = xmax + margin * dx
772 :
773 0 : dx = xmax - xmin
774 0 : if (dx < dxmin) then
775 : dx = dxmin
776 0 : xmax = (xmax + xmin) / 2 + dx / 2
777 0 : xmin = xmax - dx
778 : end if
779 :
780 0 : if (xmin == xmax) then
781 0 : xmin = xmin - margin / 2
782 0 : xmax = xmax + margin / 2
783 : end if
784 :
785 0 : if (reversed) then
786 0 : xright = xmin; xleft = xmax
787 : else
788 0 : xright = xmax; xleft = xmin
789 : end if
790 :
791 0 : if (xvec(1) < xvec(nz)) then ! increasing xs
792 0 : grid_max = nz
793 0 : do k = nz - 1, 1, -1 ! in decreasing order
794 0 : if (xvec(k) < xmax) then ! this is the first one < xmax
795 0 : grid_max = k + 1
796 0 : exit
797 : end if
798 : end do
799 0 : grid_min = 1
800 0 : do k = grid_max, 1, -1
801 0 : if (xvec(k) <= xmin) then ! this is the first one <= xmin
802 0 : grid_min = k
803 0 : exit
804 : end if
805 : end do
806 : else ! decreasing
807 0 : grid_min = 1
808 0 : do k = 2, nz ! in decreasing order
809 0 : if (xvec(k) < xmax) then ! this is the first one < xmax
810 0 : grid_min = k - 1
811 0 : exit
812 : end if
813 : end do
814 0 : grid_max = nz
815 0 : do k = grid_min, nz
816 0 : if (xvec(k) <= xmin) then ! this is the first one <= xmin
817 0 : grid_max = k
818 0 : exit
819 : end if
820 : end do
821 : end if
822 :
823 0 : end subroutine set_grid_minmax
824 :
825 :
826 0 : subroutine set_xleft_xright(&
827 0 : npts, xvec, given_xmin, given_xmax, xmargin_in, &
828 : reversed, dxmin, xleft, xright)
829 : integer, intent(in) :: npts
830 : real, intent(in), dimension(:) :: xvec
831 : real, intent(in) :: given_xmin, given_xmax, xmargin_in
832 : logical, intent(in) :: reversed
833 : real, intent(in) :: dxmin
834 : real, intent(out) :: xleft, xright
835 : call set_ytop_ybot(&
836 : npts, xvec, given_xmin, given_xmax, -101.0, xmargin_in, &
837 0 : reversed, dxmin, xleft, xright)
838 0 : end subroutine set_xleft_xright
839 :
840 :
841 0 : subroutine set_ytop_ybot(&
842 0 : npts, yvec, given_ymin, given_ymax, given_ycenter, &
843 : ymargin_in, reversed, dymin, ybot, ytop)
844 : integer, intent(in) :: npts
845 : real, intent(in), dimension(:) :: yvec
846 : real, intent(in) :: given_ymin, given_ymax, given_ycenter, ymargin_in
847 : logical, intent(in) :: reversed
848 : real, intent(in) :: dymin
849 : real, intent(out) :: ybot, ytop
850 :
851 : real :: dy, ymax, ymin, ymargin
852 : logical :: use_given_ymin, use_given_ymax
853 : real, parameter :: dymin_min = 1d-34
854 :
855 : include 'formats'
856 :
857 0 : ymargin = max(0.0, ymargin_in)
858 :
859 0 : use_given_ymin = abs(given_ymin + 101.0) > 1e-6
860 0 : if (use_given_ymin) then
861 : ymin = given_ymin
862 : else
863 0 : ymin = minval(yvec(1:npts), mask=.not. is_bad(yvec(1:npts)))
864 : end if
865 :
866 0 : use_given_ymax = abs(given_ymax + 101.0) > 1e-6
867 0 : if (use_given_ymax) then
868 : ymax = given_ymax
869 : else
870 0 : ymax = maxval(yvec(1:npts), mask=.not. is_bad(yvec(1:npts)))
871 : end if
872 0 : dy = ymax - ymin
873 :
874 0 : if (is_bad(dy)) then
875 0 : ymax = given_ymax
876 0 : ymin = given_ymin
877 0 : dy = ymax - ymin
878 : end if
879 :
880 0 : if (.not. use_given_ymin) ymin = ymin - ymargin * dy
881 0 : if (.not. use_given_ymax) ymax = ymax + ymargin * dy
882 :
883 : if (abs(given_ycenter + 101.0) > 1e-6 .and. &
884 0 : ymax > given_ycenter .and. given_ycenter > ymin) then
885 0 : dy = 2d0 * max(given_ycenter - ymin, ymax - given_ycenter)
886 0 : ymax = given_ycenter + 0.5d0 * dy
887 0 : ymin = given_ycenter - 0.5d0 * dy
888 : end if
889 :
890 0 : dy = ymax - ymin
891 0 : if (dy == 0.0) then
892 0 : ymax = ymax + dy
893 0 : ymin = ymin - dy
894 0 : else if (dy < max(dymin, dymin_min)) then
895 : !dymin_min prevents graphical glitches and segmentation faults when dy is too small
896 : dy = max(dymin, dymin_min)
897 0 : ymax = (ymax + ymin) / 2 + dy / 2
898 0 : ymin = ymax - dy
899 : end if
900 :
901 0 : if (ymin == ymax) then
902 0 : ymin = ymin - max(1.0, ymargin) / 2
903 0 : ymax = ymax + max(1.0, ymargin) / 2
904 : end if
905 :
906 0 : if (ymin == ymax) then ! round off problems
907 0 : dy = 1e-6 * abs(ymax)
908 0 : ymin = ymin - dy
909 0 : ymax = ymax + dy
910 : end if
911 :
912 0 : if (reversed) then
913 0 : ytop = ymin; ybot = ymax
914 : else
915 0 : ytop = ymax; ybot = ymin
916 : end if
917 :
918 0 : end subroutine set_ytop_ybot
919 :
920 :
921 0 : subroutine do1_pgmtxt(side, disp, coord, fjust, label, ch, lw)
922 : character (len = *), intent(in) :: side, label
923 : real, intent(in) :: disp, coord, fjust, ch
924 : integer, intent(in) :: lw
925 : real :: sav_ch
926 : integer :: sav_lw
927 0 : call pgqch(sav_ch)
928 0 : call pgqlw(sav_lw)
929 0 : call pgslw(lw)
930 0 : call pgsch(ch)
931 0 : call pgmtxt(side, disp, coord, fjust, label)
932 0 : call pgslw(sav_lw)
933 0 : call pgsch(sav_ch)
934 0 : end subroutine do1_pgmtxt
935 :
936 :
937 0 : subroutine show_annotations(s, show_annotation1, show_annotation2, show_annotation3)
938 : type (star_info), pointer :: s
939 : logical, intent(in) :: show_annotation1, show_annotation2, show_annotation3
940 0 : if (show_annotation1 .and. len_trim(s% pg% annotation1_text) > 0) then
941 0 : call pgsci(s% pg% annotation1_ci)
942 0 : call pgscf(s% pg% annotation1_cf)
943 : call do1_pgmtxt(s% pg% annotation1_side, s% pg% annotation1_disp, &
944 : s% pg% annotation1_coord, s% pg% annotation1_fjust, s% pg% annotation1_text, &
945 0 : s% pg% annotation1_ch, s% pg% annotation1_lw)
946 : end if
947 0 : if (show_annotation2 .and. len_trim(s% pg% annotation2_text) > 0) then
948 0 : call pgsci(s% pg% annotation2_ci)
949 0 : call pgscf(s% pg% annotation2_cf)
950 : call do1_pgmtxt(s% pg% annotation2_side, s% pg% annotation2_disp, &
951 : s% pg% annotation2_coord, s% pg% annotation2_fjust, s% pg% annotation2_text, &
952 0 : s% pg% annotation2_ch, s% pg% annotation2_lw)
953 : end if
954 0 : if (show_annotation3 .and. len_trim(s% pg% annotation3_text) > 0) then
955 0 : call pgsci(s% pg% annotation3_ci)
956 0 : call pgscf(s% pg% annotation3_cf)
957 : call do1_pgmtxt(s% pg% annotation3_side, s% pg% annotation3_disp, &
958 : s% pg% annotation3_coord, s% pg% annotation3_fjust, s% pg% annotation3_text, &
959 0 : s% pg% annotation3_ch, s% pg% annotation3_lw)
960 : end if
961 0 : end subroutine show_annotations
962 :
963 0 : subroutine set_xaxis_bounds(&
964 : s, xaxis_by, win_xmin_in, win_xmax_in, xaxis_reversed, xmargin, &
965 : xvec, xmin, xmax, xleft, xright, dx, &
966 : grid_min, grid_max, npts, ierr)
967 : use profile_getval, only : get_profile_id, get_profile_val
968 : type (star_info), pointer :: s
969 : character (len = *), intent(in) :: xaxis_by
970 : real, intent(in) :: win_xmin_in, win_xmax_in, xmargin
971 : logical, intent(in) :: xaxis_reversed
972 : real, allocatable, dimension(:) :: xvec
973 : real, intent(out) :: xmin, xmax, xleft, xright, dx
974 : integer, intent(out) :: grid_min, grid_max, npts
975 : integer, intent(out) :: ierr
976 :
977 : integer :: k, nz, xaxis_id
978 : real :: win_xmin, win_xmax
979 :
980 : include 'formats'
981 :
982 0 : ierr = 0
983 0 : win_xmin = win_xmin_in
984 0 : win_xmax = win_xmax_in
985 0 : nz = s% nz
986 :
987 0 : xaxis_id = get_profile_id(s, xaxis_by)
988 0 : if (xaxis_id <= 0) then
989 : write(*, '(a)') &
990 0 : 'pgstar inlist problem: bad value for xaxis_by: <' // trim(xaxis_by) // '>'
991 0 : ierr = -1
992 0 : return
993 : end if
994 0 : do k = 1, nz
995 0 : xvec(k) = get_profile_val(s, xaxis_id, k)
996 : end do
997 :
998 : call set_grid_minmax(&
999 : nz, xvec, xmin, xmax, xleft, xright, xaxis_by, &
1000 0 : win_xmin, win_xmax, xmargin, xaxis_reversed, grid_min, grid_max, 0.0)
1001 0 : dx = xmax - xmin
1002 0 : npts = grid_max - grid_min + 1
1003 0 : if (npts <= 0) then
1004 0 : write(*, *) 'invalid x axis bounds for xaxis_by = ' // trim(xaxis_by)
1005 0 : write(*, 1) 'xmax', xmax
1006 0 : write(*, 1) 'xmin', xmin
1007 0 : write(*, 1) 'dx', dx
1008 0 : write(*, 1) 'xleft', xleft
1009 0 : write(*, 1) 'xright', xright
1010 0 : write(*, 1) 'win_xmin', win_xmin
1011 0 : write(*, 1) 'win_xmax', win_xmax
1012 0 : write(*, 1) 'xmargin', xmargin
1013 0 : write(*, 2) 'grid_min', grid_min
1014 0 : write(*, 2) 'grid_max', grid_max
1015 0 : write(*, 2) 'npts', npts
1016 0 : write(*, 2) 'nz', nz
1017 0 : write(*, 1) 'maxval(xvec(1:nz))', maxval(xvec(1:nz))
1018 0 : write(*, 1) 'minval(xvec(1:nz))', minval(xvec(1:nz))
1019 0 : ierr = -1
1020 : end if
1021 :
1022 : end subroutine set_xaxis_bounds
1023 :
1024 :
1025 0 : subroutine show_xaxis_name(s, name, ierr)
1026 : type (star_info), pointer :: s
1027 : character (len = *), intent(in) :: name
1028 : integer, intent(out) :: ierr
1029 0 : ierr = 0
1030 0 : if (name == 'mass') then
1031 0 : call show_xaxis_label_pgstar(s, 'm/M\d\(2281)')
1032 0 : else if (name == 'grid') then
1033 0 : call show_xaxis_label_pgstar(s, 'grid')
1034 0 : else if (name == 'radius') then
1035 0 : call show_xaxis_label_pgstar(s, 'r/R\d\(2281)')
1036 0 : else if (name == 'logR') then
1037 0 : call show_xaxis_label_pgstar(s, 'log r/R\d\(2281)')
1038 0 : else if (name == 'logT') then
1039 0 : call show_xaxis_label_pgstar(s, 'log T')
1040 0 : else if (name == 'logP') then
1041 0 : call show_xaxis_label_pgstar(s, 'log P')
1042 0 : else if (name == 'logxq') then
1043 0 : if (s% M_center == 0) then
1044 0 : call show_xaxis_label_pgstar(s, 'log(1-q) q=fraction of total mass')
1045 : else
1046 0 : call show_xaxis_label_pgstar(s, 'log(1-q) q=fraction of envelope mass')
1047 : end if
1048 0 : else if (name == 'logxm') then
1049 0 : call show_xaxis_label_pgstar(s, 'log((Mstar-m)/M\d\(2281)\u)')
1050 0 : else if (name == 'r_div_R') then
1051 0 : call show_xaxis_label_pgstar(s, 'r/R')
1052 0 : else if (name == 'log_column_depth') then
1053 0 : call show_xaxis_label_pgstar(s, 'log column depth (g cm\u-2\d)')
1054 : else
1055 0 : call show_xaxis_label_pgstar(s, name)
1056 : end if
1057 0 : end subroutine show_xaxis_name
1058 :
1059 :
1060 0 : subroutine show_mix_regions_on_xaxis(s, ybot_in, ytop, grid_min, grid_max, xvec)
1061 : type (star_info), pointer :: s
1062 : real, intent(in) :: ybot_in, ytop
1063 : integer, intent(in) :: grid_min, grid_max
1064 : real, allocatable, dimension(:) :: xvec
1065 : real :: ybot
1066 0 : ybot = ybot_in + 0.001 * (ytop - ybot_in)
1067 0 : call show_no_mixing_section(s, ybot, grid_min, grid_max, xvec)
1068 0 : call show_convective_section(s, ybot, grid_min, grid_max, xvec)
1069 0 : call show_leftover_convective_section(s, ybot, grid_min, grid_max, xvec)
1070 0 : call show_semiconvective_section(s, ybot, grid_min, grid_max, xvec)
1071 0 : call show_thermohaline_section(s, ybot, grid_min, grid_max, xvec)
1072 0 : call show_rotation_section(s, ybot, grid_min, grid_max, xvec)
1073 0 : call show_overshoot_section(s, ybot, grid_min, grid_max, xvec)
1074 0 : end subroutine show_mix_regions_on_xaxis
1075 :
1076 :
1077 0 : subroutine show_no_mixing_section(s, ybot, grid_min, grid_max, xvec)
1078 : use pgstar_colors
1079 : type (star_info), pointer :: s
1080 : real, intent(in) :: ybot
1081 : integer, intent(in) :: grid_min, grid_max
1082 : real, allocatable, dimension(:) :: xvec
1083 0 : call show_mixing_section(s, ybot, grid_min, grid_max, xvec, no_mixing, clr_no_mixing)
1084 0 : end subroutine show_no_mixing_section
1085 :
1086 :
1087 0 : subroutine show_convective_section(s, ybot, grid_min, grid_max, xvec)
1088 : use pgstar_colors
1089 : type (star_info), pointer :: s
1090 : real, intent(in) :: ybot
1091 : integer, intent(in) :: grid_min, grid_max
1092 : real, allocatable, dimension(:) :: xvec
1093 : call show_mixing_section(&
1094 0 : s, ybot, grid_min, grid_max, xvec, convective_mixing, clr_convection)
1095 0 : end subroutine show_convective_section
1096 :
1097 :
1098 0 : subroutine show_leftover_convective_section(s, ybot, grid_min, grid_max, xvec)
1099 : use pgstar_colors
1100 : type (star_info), pointer :: s
1101 : real, intent(in) :: ybot
1102 : integer, intent(in) :: grid_min, grid_max
1103 : real, allocatable, dimension(:) :: xvec
1104 : call show_mixing_section(&
1105 0 : s, ybot, grid_min, grid_max, xvec, leftover_convective_mixing, clr_leftover_convection)
1106 0 : end subroutine show_leftover_convective_section
1107 :
1108 :
1109 0 : subroutine show_semiconvective_section(s, ybot, grid_min, grid_max, xvec)
1110 : use pgstar_colors
1111 : type (star_info), pointer :: s
1112 : real, intent(in) :: ybot
1113 : integer, intent(in) :: grid_min, grid_max
1114 : real, allocatable, dimension(:) :: xvec
1115 : call show_mixing_section(&
1116 0 : s, ybot, grid_min, grid_max, xvec, semiconvective_mixing, clr_semiconvection)
1117 0 : end subroutine show_semiconvective_section
1118 :
1119 :
1120 0 : subroutine show_thermohaline_section(s, ybot, grid_min, grid_max, xvec)
1121 : use pgstar_colors
1122 : type (star_info), pointer :: s
1123 : real, intent(in) :: ybot
1124 : integer, intent(in) :: grid_min, grid_max
1125 : real, allocatable, dimension(:) :: xvec
1126 : call show_mixing_section(&
1127 0 : s, ybot, grid_min, grid_max, xvec, thermohaline_mixing, clr_thermohaline)
1128 0 : end subroutine show_thermohaline_section
1129 :
1130 :
1131 0 : subroutine show_rotation_section(s, ybot, grid_min, grid_max, xvec)
1132 : use pgstar_colors
1133 : type (star_info), pointer :: s
1134 : real, intent(in) :: ybot
1135 : integer, intent(in) :: grid_min, grid_max
1136 : real, allocatable, dimension(:) :: xvec
1137 : call show_mixing_section(&
1138 0 : s, ybot, grid_min, grid_max, xvec, rotation_mixing, clr_rotation)
1139 0 : end subroutine show_rotation_section
1140 :
1141 :
1142 0 : subroutine show_overshoot_section(s, ybot, grid_min, grid_max, xvec)
1143 : use pgstar_colors
1144 : type (star_info), pointer :: s
1145 : real, intent(in) :: ybot
1146 : integer, intent(in) :: grid_min, grid_max
1147 : real, allocatable, dimension(:) :: xvec
1148 : call show_mixing_section(&
1149 0 : s, ybot, grid_min, grid_max, xvec, overshoot_mixing, clr_overshoot)
1150 0 : end subroutine show_overshoot_section
1151 :
1152 :
1153 0 : subroutine show_mixing_section(s, ybot, grid_min, grid_max, xvec, mixing_type, clr)
1154 : type (star_info), pointer :: s
1155 : real, intent(in) :: ybot
1156 : real, allocatable, dimension(:) :: xvec
1157 : integer, intent(in) :: mixing_type, clr, grid_min, grid_max
1158 :
1159 : integer :: k, first, last
1160 : logical :: inside
1161 : include 'formats'
1162 0 : inside = (s% mixing_type(grid_min) == mixing_type)
1163 0 : first = grid_min
1164 0 : call pgsci(clr)
1165 0 : do k = grid_min, grid_max ! 2,s% nz
1166 0 : if (.not. inside) then
1167 0 : if (s% mixing_type(k) == mixing_type) then ! starting
1168 0 : inside = .true.
1169 0 : first = k
1170 : end if
1171 : else ! inside
1172 0 : if (s% mixing_type(k) /= mixing_type) then ! ending
1173 0 : last = k - 1
1174 0 : call pgmove(xvec(first), ybot)
1175 0 : call pgdraw(xvec(last), ybot)
1176 0 : inside = .false.
1177 : end if
1178 : end if
1179 : end do
1180 0 : if (inside) then
1181 0 : last = grid_max
1182 0 : call pgmove(xvec(first), ybot)
1183 0 : call pgdraw(xvec(last), ybot)
1184 : end if
1185 0 : end subroutine show_mixing_section
1186 :
1187 :
1188 0 : subroutine show_profile_line(&
1189 0 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1190 : show_legend, legend_coord, legend_disp1, legend_del_disp, legend_fjust, &
1191 : show_mass_pts)
1192 : use pgstar_colors
1193 : type (star_info), pointer :: s
1194 : real, intent(in) :: xvec(:), yvec(:), txt_scale, xmin, xmax, ymin, ymax, &
1195 : legend_coord, legend_disp1, legend_del_disp, legend_fjust
1196 : logical, intent(in) :: show_legend, show_mass_pts
1197 :
1198 : real :: disp
1199 : integer :: nz
1200 : logical :: has_convection, has_leftover_convection, has_overshoot, &
1201 : has_semiconvection, has_thermohaline, has_rotation
1202 :
1203 : include 'formats'
1204 :
1205 0 : call pgsave
1206 :
1207 0 : nz = s% nz
1208 0 : call pgsch(s% pg% TRho_Profile_legend_txt_scale * txt_scale)
1209 :
1210 0 : call pgsci(clr_Gold)
1211 0 : call pgslw(14)
1212 0 : call do_show_eps_nuc_section(1d0)
1213 0 : call pgslw(1)
1214 0 : disp = legend_disp1 + 2 * legend_del_disp
1215 0 : if (show_legend) &
1216 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 1 erg g\u-1\d s\u-1\d')
1217 :
1218 0 : call pgsci(clr_Coral)
1219 0 : call pgslw(18)
1220 0 : call do_show_eps_nuc_section(1d3)
1221 0 : call pgslw(1)
1222 0 : disp = legend_disp1 + legend_del_disp
1223 0 : if (show_legend) &
1224 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 1000 erg g\u-1\d s\u-1\d')
1225 :
1226 0 : call pgsci(clr_Crimson)
1227 0 : call pgslw(20)
1228 0 : call do_show_eps_nuc_section(1d7)
1229 0 : call pgslw(1)
1230 0 : disp = legend_disp1
1231 0 : if (show_legend) &
1232 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 10\u7\d erg g\u-1\d s\u-1\d')
1233 :
1234 0 : disp = legend_disp1 + 2 * legend_del_disp
1235 :
1236 0 : call pgsci(clr_no_mixing)
1237 0 : call pgslw(10)
1238 0 : call pgline(nz, xvec, yvec)
1239 : has_convection = do_show_convective_section()
1240 : has_leftover_convection = do_show_leftover_convective_section()
1241 : has_overshoot = do_show_overshoot_section()
1242 : has_semiconvection = do_show_semiconvective_section()
1243 : has_thermohaline = do_show_thermohaline_section()
1244 0 : if (s% rotation_flag) then
1245 : has_rotation = do_show_rotation_section()
1246 : else
1247 : has_rotation = .false.
1248 : end if
1249 0 : call pgslw(1)
1250 0 : if (show_legend) then
1251 0 : call pgslw(1)
1252 0 : disp = disp + legend_del_disp
1253 0 : call show_legend_text(clr_no_mixing, 'no mixing')
1254 0 : if (s% rotation_flag) then
1255 0 : disp = disp + legend_del_disp
1256 0 : call show_legend_text(clr_rotation, 'rotation')
1257 : end if
1258 0 : disp = disp + legend_del_disp
1259 0 : call show_legend_text(clr_convection, 'convection')
1260 0 : disp = disp + legend_del_disp
1261 0 : call show_legend_text(clr_overshoot, 'overshoot')
1262 0 : disp = disp + legend_del_disp
1263 0 : call show_legend_text(clr_semiconvection, 'semiconvection')
1264 0 : disp = disp + legend_del_disp
1265 0 : call show_legend_text(clr_thermohaline, 'thermohaline')
1266 0 : if(s% pg% show_TRho_accretion_mesh_borders) then
1267 0 : disp = disp + legend_del_disp
1268 0 : call show_legend_text(clr_RoyalPurple, 'Lagrangian Outer Border')
1269 0 : disp = disp + legend_del_disp
1270 0 : call show_legend_text(clr_RoyalBlue, 'Homologous Inner Boundary')
1271 0 : disp = disp + legend_del_disp
1272 0 : call show_legend_text(clr_Tan, 'Mass Added This Step')
1273 : end if
1274 0 : call pgslw(10)
1275 : end if
1276 :
1277 0 : if (show_mass_pts) &
1278 0 : call show_mass_points(s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax)
1279 :
1280 0 : call pgunsa
1281 :
1282 :
1283 : contains
1284 :
1285 0 : subroutine show_legend_text(clr, txt)
1286 : integer, intent(in) :: clr
1287 : character (len = *), intent(in) :: txt
1288 0 : call pgsci(clr)
1289 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, txt)
1290 0 : end subroutine show_legend_text
1291 :
1292 :
1293 0 : subroutine do_show_eps_nuc_section(eps)
1294 : real(dp), intent(in) :: eps
1295 : integer :: k, first, last
1296 : logical :: inside
1297 0 : inside = (s% eps_nuc(1) > eps)
1298 : if (inside) first = 1
1299 0 : do k = 2, s% nz
1300 0 : if (.not. inside) then
1301 0 : if (s% eps_nuc(k) > eps) then ! starting
1302 0 : inside = .true.
1303 0 : first = k
1304 : end if
1305 : else ! inside
1306 0 : if (s% eps_nuc(k) <= eps) then ! ending
1307 0 : last = k - 1
1308 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1309 0 : inside = .false.
1310 : end if
1311 : end if
1312 : end do
1313 0 : if (inside) then
1314 0 : last = nz
1315 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1316 : end if
1317 0 : end subroutine do_show_eps_nuc_section
1318 :
1319 :
1320 0 : logical function do_show_convective_section()
1321 0 : do_show_convective_section = do_show_mixing_section(convective_mixing, clr_convection)
1322 : end function do_show_convective_section
1323 :
1324 :
1325 0 : logical function do_show_leftover_convective_section()
1326 0 : do_show_leftover_convective_section = do_show_mixing_section(leftover_convective_mixing, clr_leftover_convection)
1327 : end function do_show_leftover_convective_section
1328 :
1329 :
1330 0 : logical function do_show_semiconvective_section()
1331 0 : do_show_semiconvective_section = do_show_mixing_section(semiconvective_mixing, clr_semiconvection)
1332 : end function do_show_semiconvective_section
1333 :
1334 :
1335 0 : logical function do_show_thermohaline_section()
1336 0 : do_show_thermohaline_section = do_show_mixing_section(thermohaline_mixing, clr_thermohaline)
1337 : end function do_show_thermohaline_section
1338 :
1339 :
1340 0 : logical function do_show_rotation_section()
1341 0 : do_show_rotation_section = do_show_mixing_section(rotation_mixing, clr_rotation)
1342 : end function do_show_rotation_section
1343 :
1344 :
1345 0 : logical function do_show_overshoot_section()
1346 0 : do_show_overshoot_section = do_show_mixing_section(overshoot_mixing, clr_overshoot)
1347 : end function do_show_overshoot_section
1348 :
1349 :
1350 0 : logical function do_show_mixing_section(mixing_type, clr)
1351 : integer, intent(in) :: mixing_type, clr
1352 : integer :: k, first, last
1353 : logical :: inside
1354 : include 'formats'
1355 0 : call pgsave
1356 0 : call pgsci(clr)
1357 0 : inside = (s% mixing_type(1) == mixing_type)
1358 : if (inside) first = 1
1359 0 : do_show_mixing_section = .false.
1360 0 : do k = 2, s% nz
1361 0 : if (.not. inside) then
1362 0 : if (s% mixing_type(k) == mixing_type) then ! starting
1363 0 : inside = .true.
1364 0 : first = k
1365 : end if
1366 : else ! inside
1367 0 : if (s% mixing_type(k) /= mixing_type) then ! ending
1368 0 : last = k - 1
1369 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1370 0 : do_show_mixing_section = .true.
1371 0 : inside = .false.
1372 : end if
1373 : end if
1374 : end do
1375 0 : if (inside) then
1376 0 : last = nz
1377 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1378 0 : do_show_mixing_section = .true.
1379 : end if
1380 0 : call pgunsa
1381 0 : end function do_show_mixing_section
1382 :
1383 :
1384 : end subroutine show_profile_line
1385 :
1386 :
1387 0 : subroutine show_mass_points(s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax)
1388 : type (star_info), pointer :: s
1389 : real, intent(in) :: xvec(:), yvec(:), txt_scale, xmin, xmax, ymin, ymax
1390 : integer :: i
1391 0 : do i = 1, s% pg% num_profile_mass_points
1392 : call show_mass_point(&
1393 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1394 : s% pg% profile_mass_point_q(i), &
1395 : s% pg% profile_mass_point_color_index(i), &
1396 : s% pg% profile_mass_point_symbol(i), &
1397 : s% pg% profile_mass_point_symbol_scale(i), &
1398 : s% pg% profile_mass_point_str(i), &
1399 : s% pg% profile_mass_point_str_clr(i), &
1400 0 : s% pg% profile_mass_point_str_scale(i))
1401 : end do
1402 0 : end subroutine show_mass_points
1403 :
1404 :
1405 0 : subroutine show_mass_point(&
1406 0 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1407 : q_in, clr_index, symbol, symbol_scale, str, str_clr, str_scale)
1408 : type (star_info), pointer :: s
1409 : real, intent(in) :: xvec(:), yvec(:), txt_scale, q_in, &
1410 : xmin, xmax, ymin, ymax, symbol_scale, str_scale
1411 : integer, intent(in) :: clr_index, symbol, str_clr
1412 : character (len = *), intent(in) :: str
1413 : real :: q, q0, q1, x, y, dy
1414 : integer :: nz, i, j, k
1415 : include 'formats'
1416 0 : q = max(0.0, min(1.0, q_in))
1417 0 : nz = s% nz
1418 0 : i = nz
1419 0 : dy = ymax - ymin
1420 0 : do k = 1, s% nz - 1
1421 0 : if (s% q(k) >= q .and. q > s% q(k + 1)) then
1422 : i = k; exit
1423 : end if
1424 : end do
1425 0 : j = i + 1
1426 0 : if (j >= nz) j = i
1427 0 : q0 = s% q(i)
1428 0 : q1 = s% q(j)
1429 0 : if ((q0 - q) * (q - q1) < 0) then
1430 0 : j = i - 1
1431 0 : q1 = s% q(j)
1432 : end if
1433 0 : x = find0(xvec(i), q0 - q, xvec(j), q1 - q)
1434 0 : if (x > xmax .or. x < xmin) return
1435 0 : y = find0(yvec(i), q0 - q, yvec(j), q1 - q)
1436 0 : if (y > ymax .or. y < ymin) return
1437 0 : call pgsave
1438 0 : call pgscf(1)
1439 0 : call pgslw(1)
1440 0 : call pgsci(clr_index)
1441 0 : call pgsch(symbol_scale * txt_scale)
1442 0 : call pgpt(1, x, y, symbol)
1443 0 : call pgsci(str_clr)
1444 0 : call pgsch(str_scale * txt_scale)
1445 0 : call pgptxt(x, y - 0.015 * dy, 0.0, 0.0, trim(str))
1446 0 : call pgunsa
1447 : end subroutine show_mass_point
1448 :
1449 :
1450 0 : real function find0(xx1, yy1, xx2, yy2)
1451 : real :: xx1, yy1, xx2, yy2
1452 : real :: a, b, xz
1453 : ! returns x where y is 0 on line connecting the points (xx1,yy1) and (xx2,yy2)
1454 0 : a = (xx1 * yy2) - (xx2 * yy1)
1455 0 : b = yy2 - yy1
1456 0 : if ((abs(a) >= abs(b) * 1e30) .and. ((yy1 >= 0 .and. yy2 <= 0) &
1457 : .or. (yy1 <= 0 .and. yy2 > 0))) then
1458 0 : xz = 0.5 * (xx1 + xx2)
1459 : else
1460 0 : xz = a / b
1461 : end if
1462 0 : find0 = xz
1463 0 : end function find0
1464 :
1465 :
1466 0 : subroutine show_box_pgstar(s, str1, str2)
1467 : type (star_info), pointer :: s
1468 : character (len = *), intent(in) :: str1, str2
1469 : real :: ch
1470 : integer :: lw
1471 0 : call pgqch(ch)
1472 0 : call pgqlw(lw)
1473 0 : call pgsch(s% pg% pgstar_num_scale * ch)
1474 0 : call pgslw(s% pg% pgstar_box_lw)
1475 0 : call pgbox(str1, 0.0, 0, str2, 0.0, 0)
1476 0 : call pgsch(ch)
1477 0 : call pgslw(lw)
1478 0 : end subroutine show_box_pgstar
1479 :
1480 :
1481 0 : subroutine show_grid_title_pgstar(s, title, pad)
1482 : type (star_info), pointer :: s
1483 : character (len = *), intent(in) :: title
1484 : real, intent(in) :: pad
1485 : optional pad
1486 : real :: ch, disp
1487 0 : if (.not. s% pg% pgstar_grid_show_title) return
1488 0 : if (len_trim(title) == 0) return
1489 0 : call pgqch(ch)
1490 0 : disp = s% pg% pgstar_grid_title_disp
1491 0 : if (present(pad)) disp = disp + pad
1492 : call do1_pgmtxt('T', disp, &
1493 : s% pg% pgstar_grid_title_coord, s% pg% pgstar_grid_title_fjust, title, &
1494 0 : s% pg% pgstar_grid_title_scale * ch, s% pg% pgstar_grid_title_lw)
1495 : end subroutine show_grid_title_pgstar
1496 :
1497 :
1498 0 : subroutine show_title_pgstar(s, title, pad)
1499 : type (star_info), pointer :: s
1500 : character (len = *), intent(in) :: title
1501 : real, intent(in) :: pad
1502 : optional pad
1503 : real :: ch, disp
1504 0 : if (.not. s% pg% pgstar_show_title) return
1505 0 : if (len_trim(title) == 0) return
1506 0 : call pgqch(ch)
1507 0 : disp = s% pg% pgstar_title_disp
1508 0 : if (present(pad)) disp = disp + pad
1509 : call do1_pgmtxt('T', disp, &
1510 : s% pg% pgstar_title_coord, s% pg% pgstar_title_fjust, title, &
1511 0 : s% pg% pgstar_title_scale * ch, s% pg% pgstar_title_lw)
1512 : end subroutine show_title_pgstar
1513 :
1514 :
1515 0 : subroutine show_title_label_pgmtxt_pgstar(&
1516 : s, coord, fjust, label, pad)
1517 : type (star_info), pointer :: s
1518 : character (len = *), intent(in) :: label
1519 : real, intent(in) :: pad, coord, fjust
1520 : optional pad
1521 : real :: disp
1522 0 : disp = s% pg% pgstar_title_disp
1523 0 : if (present(pad)) disp = disp + pad
1524 0 : call pgmtxt('T', disp, coord, fjust, label)
1525 0 : end subroutine show_title_label_pgmtxt_pgstar
1526 :
1527 :
1528 0 : subroutine show_xaxis_label_pgstar(s, label, pad)
1529 : type (star_info), pointer :: s
1530 : character (len = *), intent(in) :: label
1531 : real, intent(in) :: pad
1532 : optional pad
1533 : real :: ch, disp
1534 0 : call pgqch(ch)
1535 0 : disp = s% pg% pgstar_xaxis_label_disp
1536 0 : if (present(pad)) disp = disp + pad
1537 : call do1_pgmtxt('B', disp, 0.5, 0.5, label, &
1538 0 : s% pg% pgstar_xaxis_label_scale * ch, s% pg% pgstar_xaxis_label_lw)
1539 0 : end subroutine show_xaxis_label_pgstar
1540 :
1541 :
1542 0 : subroutine show_xaxis_label_pgmtxt_pgstar(&
1543 : s, coord, fjust, label, pad)
1544 : type (star_info), pointer :: s
1545 : character (len = *), intent(in) :: label
1546 : real, intent(in) :: pad, coord, fjust
1547 : optional pad
1548 : real :: disp
1549 0 : disp = s% pg% pgstar_xaxis_label_disp
1550 0 : if (present(pad)) disp = disp + pad
1551 0 : call pgmtxt('B', disp, coord, fjust, label)
1552 0 : end subroutine show_xaxis_label_pgmtxt_pgstar
1553 :
1554 :
1555 0 : subroutine show_left_yaxis_label_pgstar(s, label, pad)
1556 : type (star_info), pointer :: s
1557 : character (len = *), intent(in) :: label
1558 : real, intent(in) :: pad
1559 : optional pad
1560 : real :: ch, disp
1561 0 : call pgqch(ch)
1562 0 : disp = s% pg% pgstar_left_yaxis_label_disp
1563 0 : if (present(pad)) disp = disp + pad
1564 : call do1_pgmtxt('L', disp, 0.5, 0.5, label, &
1565 0 : s% pg% pgstar_left_yaxis_label_scale * ch, s% pg% pgstar_left_yaxis_label_lw)
1566 0 : end subroutine show_left_yaxis_label_pgstar
1567 :
1568 :
1569 0 : subroutine show_right_yaxis_label_pgstar(s, label, pad)
1570 : type (star_info), pointer :: s
1571 : character (len = *), intent(in) :: label
1572 : real, intent(in) :: pad
1573 : optional pad
1574 : real :: ch, disp
1575 0 : call pgqch(ch)
1576 0 : disp = s% pg% pgstar_right_yaxis_label_disp
1577 0 : if (present(pad)) disp = disp + pad
1578 : call do1_pgmtxt('R', disp, 0.5, 0.5, label, &
1579 0 : s% pg% pgstar_right_yaxis_label_scale * ch, s% pg% pgstar_right_yaxis_label_lw)
1580 0 : end subroutine show_right_yaxis_label_pgstar
1581 :
1582 :
1583 0 : subroutine show_left_yaxis_label_pgmtxt_pgstar(&
1584 : s, coord, fjust, label, pad)
1585 : type (star_info), pointer :: s
1586 : character (len = *), intent(in) :: label
1587 : real, intent(in) :: pad, coord, fjust
1588 : optional pad
1589 : real :: ch, disp
1590 0 : call pgqch(ch)
1591 0 : call pgsch(1.1 * ch)
1592 0 : disp = s% pg% pgstar_left_yaxis_label_disp
1593 0 : if (present(pad)) disp = disp + pad
1594 0 : call pgmtxt('L', disp, coord, fjust, label)
1595 0 : call pgsch(ch)
1596 0 : end subroutine show_left_yaxis_label_pgmtxt_pgstar
1597 :
1598 :
1599 0 : subroutine show_right_yaxis_label_pgmtxt_pgstar(&
1600 : s, coord, fjust, label, pad)
1601 : type (star_info), pointer :: s
1602 : character (len = *), intent(in) :: label
1603 : real, intent(in) :: pad, coord, fjust
1604 : optional pad
1605 : real :: ch, disp
1606 0 : call pgqch(ch)
1607 0 : call pgsch(1.1 * ch)
1608 0 : disp = s% pg% pgstar_right_yaxis_label_disp
1609 0 : if (present(pad)) disp = disp + pad
1610 0 : call pgmtxt('R', disp, coord, fjust, label)
1611 0 : call pgsch(ch)
1612 0 : end subroutine show_right_yaxis_label_pgmtxt_pgstar
1613 :
1614 :
1615 0 : subroutine show_model_number_pgstar(s)
1616 : type (star_info), pointer :: s
1617 : character (len = 32) :: str
1618 : real :: ch
1619 0 : if (.not. s% pg% pgstar_show_model_number) return
1620 0 : write(str, '(i9)') s% model_number
1621 0 : str = 'model ' // trim(adjustl(str))
1622 0 : call pgqch(ch)
1623 : call do1_pgmtxt('T', &
1624 : s% pg% pgstar_model_disp, s% pg% pgstar_model_coord, &
1625 : s% pg% pgstar_model_fjust, str, &
1626 0 : s% pg% pgstar_model_scale * ch, s% pg% pgstar_model_lw)
1627 : end subroutine show_model_number_pgstar
1628 :
1629 :
1630 0 : subroutine show_age_pgstar(s)
1631 : type (star_info), pointer :: s
1632 : character (len = 32) :: age_str, units_str
1633 : real(dp) :: age
1634 : real :: ch
1635 : integer :: len, i, j, iE
1636 0 : if (.not. s% pg% pgstar_show_age) return
1637 0 : age = s% star_age
1638 0 : if (s% pg% pgstar_show_age_in_seconds) then
1639 0 : age = age * secyer
1640 0 : units_str = 'secs'
1641 0 : else if (s% pg% pgstar_show_age_in_minutes) then
1642 0 : age = age * secyer / 60
1643 0 : units_str = 'mins'
1644 0 : else if (s% pg% pgstar_show_age_in_hours) then
1645 0 : age = age * secyer / (60 * 60)
1646 0 : units_str = 'hrs'
1647 0 : else if (s% pg% pgstar_show_age_in_days) then
1648 0 : age = age * secyer / secday
1649 0 : units_str = 'days'
1650 0 : else if (s% pg% pgstar_show_age_in_years) then
1651 : !age = age
1652 0 : units_str = 'yrs'
1653 0 : else if (s% pg% pgstar_show_log_age_in_years) then
1654 0 : age = log10(max(1d-99, age))
1655 0 : units_str = 'log yrs'
1656 0 : else if (age * secyer < 60) then
1657 0 : age = age * secyer
1658 0 : units_str = 'secs'
1659 0 : else if (age * secyer < 60 * 60) then
1660 0 : age = age * secyer / 60
1661 0 : units_str = 'mins'
1662 0 : else if (age * secyer < secday) then
1663 0 : age = age * secyer / (60 * 60)
1664 0 : units_str = 'hrs'
1665 0 : else if (age * secyer < secday * 500) then
1666 0 : age = age * secyer / secday
1667 0 : units_str = 'days'
1668 : else
1669 : !age = age
1670 0 : units_str = 'yrs'
1671 : end if
1672 0 : if (abs(age) > 1e-3 .and. abs(age) < 1e3) then
1673 0 : write(age_str, '(f14.6)') age
1674 : else
1675 0 : write(age_str, '(1pe14.6)') age
1676 0 : len = len_trim(age_str)
1677 0 : iE = 0
1678 0 : do i = 1, len
1679 0 : if (age_str(i:i) == 'E') then
1680 0 : iE = i
1681 0 : age_str(i:i) = 'e'
1682 0 : exit
1683 : end if
1684 : end do
1685 0 : if (iE > 0) then
1686 0 : i = iE + 1
1687 0 : if (age_str(i:i) == '+') then
1688 0 : do j = i, len - 1
1689 0 : age_str(j:j) = age_str(j + 1:j + 1)
1690 : end do
1691 0 : age_str(len:len) = ' '
1692 0 : len = len - 1
1693 : else
1694 0 : i = i + 1
1695 : end if
1696 0 : if (age_str(i:i) == '0') then
1697 0 : do j = i, len - 1
1698 0 : age_str(j:j) = age_str(j + 1:j + 1)
1699 : end do
1700 0 : age_str(len:len) = ' '
1701 0 : len = len - 1
1702 : end if
1703 : end if
1704 : end if
1705 0 : age_str = adjustl(age_str)
1706 0 : age_str = 'age ' // trim(age_str) // ' ' // trim(units_str)
1707 0 : call pgqch(ch)
1708 : call do1_pgmtxt('T', &
1709 : s% pg% pgstar_age_disp, s% pg% pgstar_age_coord, &
1710 : s% pg% pgstar_age_fjust, age_str, &
1711 0 : s% pg% pgstar_age_scale * ch, s% pg% pgstar_age_lw)
1712 : end subroutine show_age_pgstar
1713 :
1714 :
1715 0 : logical function read_values_from_file(fname, x_data, y_data, data_len)
1716 : character(len = *), intent(in) :: fname
1717 : real, allocatable, dimension(:) :: x_data, y_data
1718 : integer, intent(out) :: data_len
1719 : integer :: iounit, ierr, i
1720 : include 'formats'
1721 0 : read_values_from_file = .false.
1722 0 : ierr = 0
1723 0 : open(newunit = iounit, file = trim(fname), action = 'read', status = 'old', iostat = ierr)
1724 0 : if (ierr /= 0) then
1725 : !write(*, *) 'failed to open ' // trim(fname)
1726 : return
1727 : end if
1728 0 : read(iounit, *, iostat = ierr) data_len
1729 0 : if (ierr /= 0) then
1730 0 : write(*, *) 'failed to read num points on 1st line ' // trim(fname)
1731 0 : return
1732 : end if
1733 : !write(*,2) trim(fname) // ' data_len', data_len
1734 0 : allocate(x_data(data_len), y_data(data_len))
1735 0 : do i = 1, data_len
1736 0 : read(iounit, *, iostat = ierr) x_data(i), y_data(i)
1737 0 : if (ierr /= 0) then
1738 0 : write(*, *) 'failed to read data ' // trim(fname)
1739 0 : deallocate(x_data, y_data)
1740 0 : return
1741 : end if
1742 : end do
1743 0 : close(iounit)
1744 0 : read_values_from_file = .true.
1745 0 : end function read_values_from_file
1746 :
1747 :
1748 0 : subroutine show_pgstar_decorator(id, use_flag, pgstar_decorator, plot_num, ierr)
1749 : logical, intent(in) :: use_flag
1750 : real :: xmin, xmax, ymin, ymax
1751 : integer, intent(in) :: id, plot_num
1752 : integer, intent(inout) :: ierr
1753 : procedure(pgstar_decorator_interface), pointer :: pgstar_decorator
1754 :
1755 0 : if(use_flag)then
1756 0 : if(associated(pgstar_decorator))then
1757 0 : call pgsave
1758 0 : call PGQWIN(xmin, xmax, ymin, ymax)
1759 0 : call pgstar_decorator(id, xmin, xmax, ymin, ymax, plot_num, ierr)
1760 0 : call pgunsa
1761 0 : if(ierr/=0)then
1762 0 : write(*, *) "Error in pgstar_decorator"
1763 : end if
1764 : end if
1765 : end if
1766 :
1767 0 : end subroutine show_pgstar_decorator
1768 :
1769 : end module pgstar_support
1770 :
|
