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 0 : 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 0 : real :: width, ratio
248 :
249 0 : if (is_file) then
250 0 : dir = p% file_dir
251 0 : white_on_black_flag = s% pg% file_white_on_black_flag
252 : else
253 0 : 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 0 : 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 0 : 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 0 : index = i
614 0 : return
615 : end if
616 : end do
617 0 : 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 0 : 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 0 : 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)) + 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))
750 : else
751 0 : xmin = minval(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)) + given_xmax
757 0 : else if (use_given_xmax) then
758 0 : xmax = given_xmax
759 : else
760 0 : xmax = maxval(xvec(1:nz))
761 : end if
762 0 : dx = xmax - xmin
763 :
764 0 : if (.not. use_given_xmin) xmin = xmin - margin * dx
765 0 : if (.not. use_given_xmax) xmax = xmax + margin * dx
766 :
767 0 : dx = xmax - xmin
768 0 : if (dx < dxmin) then
769 0 : dx = dxmin
770 0 : xmax = (xmax + xmin) / 2 + dx / 2
771 0 : xmin = xmax - dx
772 : end if
773 :
774 0 : if (xmin == xmax) then
775 0 : xmin = xmin - margin / 2
776 0 : xmax = xmax + margin / 2
777 : end if
778 :
779 0 : if (reversed) then
780 0 : xright = xmin; xleft = xmax
781 : else
782 0 : xright = xmax; xleft = xmin
783 : end if
784 :
785 0 : if (xvec(1) < xvec(nz)) then ! increasing xs
786 0 : grid_max = nz
787 0 : do k = nz - 1, 1, -1 ! in decreasing order
788 0 : if (xvec(k) < xmax) then ! this is the first one < xmax
789 0 : grid_max = k + 1
790 0 : exit
791 : end if
792 : end do
793 0 : grid_min = 1
794 0 : do k = grid_max, 1, -1
795 0 : if (xvec(k) <= xmin) then ! this is the first one <= xmin
796 0 : grid_min = k
797 0 : exit
798 : end if
799 : end do
800 : else ! decreasing
801 0 : grid_min = 1
802 0 : do k = 2, nz ! in decreasing order
803 0 : if (xvec(k) < xmax) then ! this is the first one < xmax
804 0 : grid_min = k - 1
805 0 : exit
806 : end if
807 : end do
808 0 : grid_max = nz
809 0 : do k = grid_min, nz
810 0 : if (xvec(k) <= xmin) then ! this is the first one <= xmin
811 0 : grid_max = k
812 0 : exit
813 : end if
814 : end do
815 : end if
816 :
817 0 : end subroutine set_grid_minmax
818 :
819 :
820 0 : subroutine set_xleft_xright(&
821 0 : npts, xvec, given_xmin, given_xmax, xmargin_in, &
822 : reversed, dxmin, xleft, xright)
823 : integer, intent(in) :: npts
824 : real, intent(in), dimension(:) :: xvec
825 : real, intent(in) :: given_xmin, given_xmax, xmargin_in
826 : logical, intent(in) :: reversed
827 : real, intent(in) :: dxmin
828 : real, intent(out) :: xleft, xright
829 : call set_ytop_ybot(&
830 : npts, xvec, given_xmin, given_xmax, -101.0, xmargin_in, &
831 0 : reversed, dxmin, xleft, xright)
832 0 : end subroutine set_xleft_xright
833 :
834 :
835 0 : subroutine set_ytop_ybot(&
836 0 : npts, yvec, given_ymin, given_ymax, given_ycenter, &
837 : ymargin_in, reversed, dymin, ybot, ytop)
838 : integer, intent(in) :: npts
839 : real, intent(in), dimension(:) :: yvec
840 : real, intent(in) :: given_ymin, given_ymax, given_ycenter, ymargin_in
841 : logical, intent(in) :: reversed
842 : real, intent(in) :: dymin
843 : real, intent(out) :: ybot, ytop
844 :
845 0 : real :: dy, ymax, ymin, ymargin
846 : logical :: use_given_ymin, use_given_ymax
847 : real, parameter :: dymin_min = 1d-34
848 :
849 : include 'formats'
850 :
851 0 : ymargin = max(0.0, ymargin_in)
852 :
853 0 : use_given_ymin = abs(given_ymin + 101.0) > 1e-6
854 0 : if (use_given_ymin) then
855 : ymin = given_ymin
856 : else
857 0 : ymin = minval(yvec(1:npts))
858 : end if
859 :
860 0 : use_given_ymax = abs(given_ymax + 101.0) > 1e-6
861 0 : if (use_given_ymax) then
862 : ymax = given_ymax
863 : else
864 0 : ymax = maxval(yvec(1:npts))
865 : end if
866 0 : dy = ymax - ymin
867 :
868 0 : if (.not. use_given_ymin) ymin = ymin - ymargin * dy
869 0 : if (.not. use_given_ymax) ymax = ymax + ymargin * dy
870 :
871 : if (abs(given_ycenter + 101.0) > 1e-6 .and. &
872 0 : ymax > given_ycenter .and. given_ycenter > ymin) then
873 0 : dy = 2d0 * max(given_ycenter - ymin, ymax - given_ycenter)
874 0 : ymax = given_ycenter + 0.5d0 * dy
875 0 : ymin = given_ycenter - 0.5d0 * dy
876 : end if
877 :
878 0 : dy = ymax - ymin
879 0 : if (dy == 0.0) then
880 0 : ymax = ymax + dy
881 0 : ymin = ymin - dy
882 0 : else if (dy < max(dymin, dymin_min)) then
883 : !dymin_min prevents graphical glitches and segmentation faults when dy is too small
884 0 : dy = max(dymin, dymin_min)
885 0 : ymax = (ymax + ymin) / 2 + dy / 2
886 0 : ymin = ymax - dy
887 : end if
888 :
889 0 : if (ymin == ymax) then
890 0 : ymin = ymin - max(1.0, ymargin) / 2
891 0 : ymax = ymax + max(1.0, ymargin) / 2
892 : end if
893 :
894 0 : if (ymin == ymax) then ! round off problems
895 0 : dy = 1e-6 * abs(ymax)
896 0 : ymin = ymin - dy
897 0 : ymax = ymax + dy
898 : end if
899 :
900 0 : if (reversed) then
901 0 : ytop = ymin; ybot = ymax
902 : else
903 0 : ytop = ymax; ybot = ymin
904 : end if
905 :
906 0 : end subroutine set_ytop_ybot
907 :
908 :
909 0 : subroutine do1_pgmtxt(side, disp, coord, fjust, label, ch, lw)
910 : character (len = *), intent(in) :: side, label
911 : real, intent(in) :: disp, coord, fjust, ch
912 : integer, intent(in) :: lw
913 0 : real :: sav_ch
914 : integer :: sav_lw
915 0 : call pgqch(sav_ch)
916 0 : call pgqlw(sav_lw)
917 0 : call pgslw(lw)
918 0 : call pgsch(ch)
919 0 : call pgmtxt(side, disp, coord, fjust, label)
920 0 : call pgslw(sav_lw)
921 0 : call pgsch(sav_ch)
922 0 : end subroutine do1_pgmtxt
923 :
924 :
925 0 : subroutine show_annotations(s, show_annotation1, show_annotation2, show_annotation3)
926 : type (star_info), pointer :: s
927 : logical, intent(in) :: show_annotation1, show_annotation2, show_annotation3
928 0 : if (show_annotation1 .and. len_trim(s% pg% annotation1_text) > 0) then
929 0 : call pgsci(s% pg% annotation1_ci)
930 0 : call pgscf(s% pg% annotation1_cf)
931 : call do1_pgmtxt(s% pg% annotation1_side, s% pg% annotation1_disp, &
932 : s% pg% annotation1_coord, s% pg% annotation1_fjust, s% pg% annotation1_text, &
933 0 : s% pg% annotation1_ch, s% pg% annotation1_lw)
934 : end if
935 0 : if (show_annotation2 .and. len_trim(s% pg% annotation2_text) > 0) then
936 0 : call pgsci(s% pg% annotation2_ci)
937 0 : call pgscf(s% pg% annotation2_cf)
938 : call do1_pgmtxt(s% pg% annotation2_side, s% pg% annotation2_disp, &
939 : s% pg% annotation2_coord, s% pg% annotation2_fjust, s% pg% annotation2_text, &
940 0 : s% pg% annotation2_ch, s% pg% annotation2_lw)
941 : end if
942 0 : if (show_annotation3 .and. len_trim(s% pg% annotation3_text) > 0) then
943 0 : call pgsci(s% pg% annotation3_ci)
944 0 : call pgscf(s% pg% annotation3_cf)
945 : call do1_pgmtxt(s% pg% annotation3_side, s% pg% annotation3_disp, &
946 : s% pg% annotation3_coord, s% pg% annotation3_fjust, s% pg% annotation3_text, &
947 0 : s% pg% annotation3_ch, s% pg% annotation3_lw)
948 : end if
949 0 : end subroutine show_annotations
950 :
951 0 : subroutine set_xaxis_bounds(&
952 : s, xaxis_by, win_xmin_in, win_xmax_in, xaxis_reversed, xmargin, &
953 : xvec, xmin, xmax, xleft, xright, dx, &
954 : grid_min, grid_max, npts, ierr)
955 : use profile_getval, only : get_profile_id, get_profile_val
956 : type (star_info), pointer :: s
957 : character (len = *), intent(in) :: xaxis_by
958 : real, intent(in) :: win_xmin_in, win_xmax_in, xmargin
959 : logical, intent(in) :: xaxis_reversed
960 : real, allocatable, dimension(:) :: xvec
961 : real, intent(out) :: xmin, xmax, xleft, xright, dx
962 : integer, intent(out) :: grid_min, grid_max, npts
963 : integer, intent(out) :: ierr
964 :
965 : integer :: k, nz, xaxis_id
966 : real :: win_xmin, win_xmax
967 :
968 : include 'formats'
969 :
970 0 : ierr = 0
971 0 : win_xmin = win_xmin_in
972 0 : win_xmax = win_xmax_in
973 0 : nz = s% nz
974 :
975 0 : xaxis_id = get_profile_id(s, xaxis_by)
976 0 : if (xaxis_id <= 0) then
977 : write(*, '(a)') &
978 0 : 'pgstar inlist problem: bad value for xaxis_by: <' // trim(xaxis_by) // '>'
979 0 : ierr = -1
980 : return
981 : end if
982 0 : do k = 1, nz
983 0 : xvec(k) = get_profile_val(s, xaxis_id, k)
984 : end do
985 :
986 : call set_grid_minmax(&
987 : nz, xvec, xmin, xmax, xleft, xright, xaxis_by, &
988 0 : win_xmin, win_xmax, xmargin, xaxis_reversed, grid_min, grid_max, 0.0)
989 0 : dx = xmax - xmin
990 0 : npts = grid_max - grid_min + 1
991 0 : if (npts <= 0) then
992 0 : write(*, *) 'invalid x axis bounds for xaxis_by = ' // trim(xaxis_by)
993 0 : write(*, 1) 'xmax', xmax
994 0 : write(*, 1) 'xmin', xmin
995 0 : write(*, 1) 'dx', dx
996 0 : write(*, 1) 'xleft', xleft
997 0 : write(*, 1) 'xright', xright
998 0 : write(*, 1) 'win_xmin', win_xmin
999 0 : write(*, 1) 'win_xmax', win_xmax
1000 0 : write(*, 1) 'xmargin', xmargin
1001 0 : write(*, 2) 'grid_min', grid_min
1002 0 : write(*, 2) 'grid_max', grid_max
1003 0 : write(*, 2) 'npts', npts
1004 0 : write(*, 2) 'nz', nz
1005 0 : write(*, 1) 'maxval(xvec(1:nz))', maxval(xvec(1:nz))
1006 0 : write(*, 1) 'minval(xvec(1:nz))', minval(xvec(1:nz))
1007 0 : ierr = -1
1008 : end if
1009 :
1010 0 : end subroutine set_xaxis_bounds
1011 :
1012 :
1013 0 : subroutine show_xaxis_name(s, name, ierr)
1014 : type (star_info), pointer :: s
1015 : character (len = *), intent(in) :: name
1016 : integer, intent(out) :: ierr
1017 0 : ierr = 0
1018 0 : if (name == 'mass') then
1019 0 : call show_xaxis_label_pgstar(s, 'm/M\d\(2281)')
1020 0 : else if (name == 'grid') then
1021 0 : call show_xaxis_label_pgstar(s, 'grid')
1022 0 : else if (name == 'radius') then
1023 0 : call show_xaxis_label_pgstar(s, 'r/R\d\(2281)')
1024 0 : else if (name == 'logR') then
1025 0 : call show_xaxis_label_pgstar(s, 'log r/R\d\(2281)')
1026 0 : else if (name == 'logT') then
1027 0 : call show_xaxis_label_pgstar(s, 'log T')
1028 0 : else if (name == 'logP') then
1029 0 : call show_xaxis_label_pgstar(s, 'log P')
1030 0 : else if (name == 'logxq') then
1031 0 : if (s% M_center == 0) then
1032 0 : call show_xaxis_label_pgstar(s, 'log(1-q) q=fraction of total mass')
1033 : else
1034 0 : call show_xaxis_label_pgstar(s, 'log(1-q) q=fraction of envelope mass')
1035 : end if
1036 0 : else if (name == 'logxm') then
1037 0 : call show_xaxis_label_pgstar(s, 'log((Mstar-m)/M\d\(2281)\u)')
1038 0 : else if (name == 'r_div_R') then
1039 0 : call show_xaxis_label_pgstar(s, 'r/R')
1040 0 : else if (name == 'log_column_depth') then
1041 0 : call show_xaxis_label_pgstar(s, 'log column depth (g cm\u-2\d)')
1042 : else
1043 0 : call show_xaxis_label_pgstar(s, name)
1044 : end if
1045 0 : end subroutine show_xaxis_name
1046 :
1047 :
1048 0 : subroutine show_mix_regions_on_xaxis(s, ybot_in, ytop, grid_min, grid_max, xvec)
1049 : type (star_info), pointer :: s
1050 : real, intent(in) :: ybot_in, ytop
1051 : integer, intent(in) :: grid_min, grid_max
1052 : real, allocatable, dimension(:) :: xvec
1053 : real :: ybot
1054 0 : ybot = ybot_in + 0.001 * (ytop - ybot_in)
1055 0 : call show_no_mixing_section(s, ybot, grid_min, grid_max, xvec)
1056 0 : call show_convective_section(s, ybot, grid_min, grid_max, xvec)
1057 0 : call show_leftover_convective_section(s, ybot, grid_min, grid_max, xvec)
1058 0 : call show_semiconvective_section(s, ybot, grid_min, grid_max, xvec)
1059 0 : call show_thermohaline_section(s, ybot, grid_min, grid_max, xvec)
1060 0 : call show_rotation_section(s, ybot, grid_min, grid_max, xvec)
1061 0 : call show_overshoot_section(s, ybot, grid_min, grid_max, xvec)
1062 0 : end subroutine show_mix_regions_on_xaxis
1063 :
1064 :
1065 0 : subroutine show_no_mixing_section(s, ybot, grid_min, grid_max, xvec)
1066 : use pgstar_colors
1067 : type (star_info), pointer :: s
1068 : real, intent(in) :: ybot
1069 : integer, intent(in) :: grid_min, grid_max
1070 : real, allocatable, dimension(:) :: xvec
1071 0 : call show_mixing_section(s, ybot, grid_min, grid_max, xvec, no_mixing, clr_no_mixing)
1072 0 : end subroutine show_no_mixing_section
1073 :
1074 :
1075 0 : subroutine show_convective_section(s, ybot, grid_min, grid_max, xvec)
1076 : use pgstar_colors
1077 : type (star_info), pointer :: s
1078 : real, intent(in) :: ybot
1079 : integer, intent(in) :: grid_min, grid_max
1080 : real, allocatable, dimension(:) :: xvec
1081 : call show_mixing_section(&
1082 0 : s, ybot, grid_min, grid_max, xvec, convective_mixing, clr_convection)
1083 0 : end subroutine show_convective_section
1084 :
1085 :
1086 0 : subroutine show_leftover_convective_section(s, ybot, grid_min, grid_max, xvec)
1087 : use pgstar_colors
1088 : type (star_info), pointer :: s
1089 : real, intent(in) :: ybot
1090 : integer, intent(in) :: grid_min, grid_max
1091 : real, allocatable, dimension(:) :: xvec
1092 : call show_mixing_section(&
1093 0 : s, ybot, grid_min, grid_max, xvec, leftover_convective_mixing, clr_leftover_convection)
1094 0 : end subroutine show_leftover_convective_section
1095 :
1096 :
1097 0 : subroutine show_semiconvective_section(s, ybot, grid_min, grid_max, xvec)
1098 : use pgstar_colors
1099 : type (star_info), pointer :: s
1100 : real, intent(in) :: ybot
1101 : integer, intent(in) :: grid_min, grid_max
1102 : real, allocatable, dimension(:) :: xvec
1103 : call show_mixing_section(&
1104 0 : s, ybot, grid_min, grid_max, xvec, semiconvective_mixing, clr_semiconvection)
1105 0 : end subroutine show_semiconvective_section
1106 :
1107 :
1108 0 : subroutine show_thermohaline_section(s, ybot, grid_min, grid_max, xvec)
1109 : use pgstar_colors
1110 : type (star_info), pointer :: s
1111 : real, intent(in) :: ybot
1112 : integer, intent(in) :: grid_min, grid_max
1113 : real, allocatable, dimension(:) :: xvec
1114 : call show_mixing_section(&
1115 0 : s, ybot, grid_min, grid_max, xvec, thermohaline_mixing, clr_thermohaline)
1116 0 : end subroutine show_thermohaline_section
1117 :
1118 :
1119 0 : subroutine show_rotation_section(s, ybot, grid_min, grid_max, xvec)
1120 : use pgstar_colors
1121 : type (star_info), pointer :: s
1122 : real, intent(in) :: ybot
1123 : integer, intent(in) :: grid_min, grid_max
1124 : real, allocatable, dimension(:) :: xvec
1125 : call show_mixing_section(&
1126 0 : s, ybot, grid_min, grid_max, xvec, rotation_mixing, clr_rotation)
1127 0 : end subroutine show_rotation_section
1128 :
1129 :
1130 0 : subroutine show_overshoot_section(s, ybot, grid_min, grid_max, xvec)
1131 : use pgstar_colors
1132 : type (star_info), pointer :: s
1133 : real, intent(in) :: ybot
1134 : integer, intent(in) :: grid_min, grid_max
1135 : real, allocatable, dimension(:) :: xvec
1136 : call show_mixing_section(&
1137 0 : s, ybot, grid_min, grid_max, xvec, overshoot_mixing, clr_overshoot)
1138 0 : end subroutine show_overshoot_section
1139 :
1140 :
1141 0 : subroutine show_mixing_section(s, ybot, grid_min, grid_max, xvec, mixing_type, clr)
1142 : type (star_info), pointer :: s
1143 : real, intent(in) :: ybot
1144 : real, allocatable, dimension(:) :: xvec
1145 : integer, intent(in) :: mixing_type, clr, grid_min, grid_max
1146 :
1147 : integer :: k, first, last
1148 : logical :: inside
1149 : include 'formats'
1150 0 : inside = (s% mixing_type(grid_min) == mixing_type)
1151 0 : first = grid_min
1152 0 : call pgsci(clr)
1153 0 : do k = grid_min, grid_max ! 2,s% nz
1154 0 : if (.not. inside) then
1155 0 : if (s% mixing_type(k) == mixing_type) then ! starting
1156 0 : inside = .true.
1157 0 : first = k
1158 : end if
1159 : else ! inside
1160 0 : if (s% mixing_type(k) /= mixing_type) then ! ending
1161 0 : last = k - 1
1162 0 : call pgmove(xvec(first), ybot)
1163 0 : call pgdraw(xvec(last), ybot)
1164 0 : inside = .false.
1165 : end if
1166 : end if
1167 : end do
1168 0 : if (inside) then
1169 0 : last = grid_max
1170 0 : call pgmove(xvec(first), ybot)
1171 0 : call pgdraw(xvec(last), ybot)
1172 : end if
1173 0 : end subroutine show_mixing_section
1174 :
1175 :
1176 0 : subroutine show_profile_line(&
1177 0 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1178 : show_legend, legend_coord, legend_disp1, legend_del_disp, legend_fjust, &
1179 : show_mass_pts)
1180 : use pgstar_colors
1181 : type (star_info), pointer :: s
1182 : real, intent(in) :: xvec(:), yvec(:), txt_scale, xmin, xmax, ymin, ymax, &
1183 : legend_coord, legend_disp1, legend_del_disp, legend_fjust
1184 : logical, intent(in) :: show_legend, show_mass_pts
1185 :
1186 0 : real :: disp
1187 : integer :: nz
1188 : logical :: has_convection, has_leftover_convection, has_overshoot, &
1189 : has_semiconvection, has_thermohaline, has_rotation
1190 :
1191 : include 'formats'
1192 :
1193 0 : call pgsave
1194 :
1195 0 : nz = s% nz
1196 0 : call pgsch(s% pg% TRho_Profile_legend_txt_scale * txt_scale)
1197 :
1198 0 : call pgsci(clr_Gold)
1199 0 : call pgslw(14)
1200 0 : call do_show_eps_nuc_section(1d0)
1201 0 : call pgslw(1)
1202 0 : disp = legend_disp1 + 2 * legend_del_disp
1203 0 : if (show_legend) &
1204 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 1 erg g\u-1\d s\u-1\d')
1205 :
1206 0 : call pgsci(clr_Coral)
1207 0 : call pgslw(18)
1208 0 : call do_show_eps_nuc_section(1d3)
1209 0 : call pgslw(1)
1210 0 : disp = legend_disp1 + legend_del_disp
1211 0 : if (show_legend) &
1212 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 1000 erg g\u-1\d s\u-1\d')
1213 :
1214 0 : call pgsci(clr_Crimson)
1215 0 : call pgslw(20)
1216 0 : call do_show_eps_nuc_section(1d7)
1217 0 : call pgslw(1)
1218 0 : disp = legend_disp1
1219 0 : if (show_legend) &
1220 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, '> 10\u7\d erg g\u-1\d s\u-1\d')
1221 :
1222 0 : disp = legend_disp1 + 2 * legend_del_disp
1223 :
1224 0 : call pgsci(clr_no_mixing)
1225 0 : call pgslw(10)
1226 0 : call pgline(nz, xvec, yvec)
1227 : has_convection = do_show_convective_section()
1228 : has_leftover_convection = do_show_leftover_convective_section()
1229 : has_overshoot = do_show_overshoot_section()
1230 : has_semiconvection = do_show_semiconvective_section()
1231 : has_thermohaline = do_show_thermohaline_section()
1232 0 : if (s% rotation_flag) then
1233 : has_rotation = do_show_rotation_section()
1234 : else
1235 : has_rotation = .false.
1236 : end if
1237 0 : call pgslw(1)
1238 0 : if (show_legend) then
1239 0 : call pgslw(1)
1240 0 : disp = disp + legend_del_disp
1241 0 : call show_legend_text(clr_no_mixing, 'no mixing')
1242 0 : if (s% rotation_flag) then
1243 0 : disp = disp + legend_del_disp
1244 0 : call show_legend_text(clr_rotation, 'rotation')
1245 : end if
1246 0 : disp = disp + legend_del_disp
1247 0 : call show_legend_text(clr_convection, 'convection')
1248 0 : disp = disp + legend_del_disp
1249 0 : call show_legend_text(clr_overshoot, 'overshoot')
1250 0 : disp = disp + legend_del_disp
1251 0 : call show_legend_text(clr_semiconvection, 'semiconvection')
1252 0 : disp = disp + legend_del_disp
1253 0 : call show_legend_text(clr_thermohaline, 'thermohaline')
1254 0 : if(s% pg% show_TRho_accretion_mesh_borders) then
1255 0 : disp = disp + legend_del_disp
1256 0 : call show_legend_text(clr_RoyalPurple, 'Lagrangian Outer Border')
1257 0 : disp = disp + legend_del_disp
1258 0 : call show_legend_text(clr_RoyalBlue, 'Homologous Inner Boundary')
1259 0 : disp = disp + legend_del_disp
1260 0 : call show_legend_text(clr_Tan, 'Mass Added This Step')
1261 : end if
1262 0 : call pgslw(10)
1263 : end if
1264 :
1265 0 : if (show_mass_pts) &
1266 0 : call show_mass_points(s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax)
1267 :
1268 0 : call pgunsa
1269 :
1270 :
1271 : contains
1272 :
1273 0 : subroutine show_legend_text(clr, txt)
1274 : integer, intent(in) :: clr
1275 : character (len = *), intent(in) :: txt
1276 0 : call pgsci(clr)
1277 0 : call pgmtxt('T', disp, legend_coord, legend_fjust, txt)
1278 0 : end subroutine show_legend_text
1279 :
1280 :
1281 0 : subroutine do_show_eps_nuc_section(eps)
1282 : real(dp), intent(in) :: eps
1283 : integer :: k, first, last
1284 : logical :: inside
1285 0 : inside = (s% eps_nuc(1) > eps)
1286 : if (inside) first = 1
1287 0 : do k = 2, s% nz
1288 0 : if (.not. inside) then
1289 0 : if (s% eps_nuc(k) > eps) then ! starting
1290 0 : inside = .true.
1291 0 : first = k
1292 : end if
1293 : else ! inside
1294 0 : if (s% eps_nuc(k) <= eps) then ! ending
1295 0 : last = k - 1
1296 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1297 0 : inside = .false.
1298 : end if
1299 : end if
1300 : end do
1301 0 : if (inside) then
1302 0 : last = nz
1303 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1304 : end if
1305 0 : end subroutine do_show_eps_nuc_section
1306 :
1307 :
1308 0 : logical function do_show_convective_section()
1309 0 : do_show_convective_section = do_show_mixing_section(convective_mixing, clr_convection)
1310 : end function do_show_convective_section
1311 :
1312 :
1313 0 : logical function do_show_leftover_convective_section()
1314 0 : do_show_leftover_convective_section = do_show_mixing_section(leftover_convective_mixing, clr_leftover_convection)
1315 : end function do_show_leftover_convective_section
1316 :
1317 :
1318 0 : logical function do_show_semiconvective_section()
1319 0 : do_show_semiconvective_section = do_show_mixing_section(semiconvective_mixing, clr_semiconvection)
1320 : end function do_show_semiconvective_section
1321 :
1322 :
1323 0 : logical function do_show_thermohaline_section()
1324 0 : do_show_thermohaline_section = do_show_mixing_section(thermohaline_mixing, clr_thermohaline)
1325 : end function do_show_thermohaline_section
1326 :
1327 :
1328 0 : logical function do_show_rotation_section()
1329 0 : do_show_rotation_section = do_show_mixing_section(rotation_mixing, clr_rotation)
1330 : end function do_show_rotation_section
1331 :
1332 :
1333 0 : logical function do_show_overshoot_section()
1334 0 : do_show_overshoot_section = do_show_mixing_section(overshoot_mixing, clr_overshoot)
1335 : end function do_show_overshoot_section
1336 :
1337 :
1338 0 : logical function do_show_mixing_section(mixing_type, clr)
1339 : integer, intent(in) :: mixing_type, clr
1340 : integer :: k, first, last
1341 : logical :: inside
1342 : include 'formats'
1343 0 : call pgsave
1344 0 : call pgsci(clr)
1345 0 : inside = (s% mixing_type(1) == mixing_type)
1346 : if (inside) first = 1
1347 0 : do_show_mixing_section = .false.
1348 0 : do k = 2, s% nz
1349 0 : if (.not. inside) then
1350 0 : if (s% mixing_type(k) == mixing_type) then ! starting
1351 0 : inside = .true.
1352 0 : first = k
1353 : end if
1354 : else ! inside
1355 0 : if (s% mixing_type(k) /= mixing_type) then ! ending
1356 0 : last = k - 1
1357 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1358 0 : do_show_mixing_section = .true.
1359 0 : inside = .false.
1360 : end if
1361 : end if
1362 : end do
1363 0 : if (inside) then
1364 0 : last = nz
1365 0 : call pgline(k - first, xvec(first:last), yvec(first:last))
1366 0 : do_show_mixing_section = .true.
1367 : end if
1368 0 : call pgunsa
1369 0 : end function do_show_mixing_section
1370 :
1371 :
1372 : end subroutine show_profile_line
1373 :
1374 :
1375 0 : subroutine show_mass_points(s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax)
1376 : type (star_info), pointer :: s
1377 : real, intent(in) :: xvec(:), yvec(:), txt_scale, xmin, xmax, ymin, ymax
1378 : integer :: i
1379 0 : do i = 1, s% pg% num_profile_mass_points
1380 : call show_mass_point(&
1381 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1382 : s% pg% profile_mass_point_q(i), &
1383 : s% pg% profile_mass_point_color_index(i), &
1384 : s% pg% profile_mass_point_symbol(i), &
1385 : s% pg% profile_mass_point_symbol_scale(i), &
1386 : s% pg% profile_mass_point_str(i), &
1387 : s% pg% profile_mass_point_str_clr(i), &
1388 0 : s% pg% profile_mass_point_str_scale(i))
1389 : end do
1390 0 : end subroutine show_mass_points
1391 :
1392 :
1393 0 : subroutine show_mass_point(&
1394 0 : s, xvec, yvec, txt_scale, xmin, xmax, ymin, ymax, &
1395 : q_in, clr_index, symbol, symbol_scale, str, str_clr, str_scale)
1396 : type (star_info), pointer :: s
1397 : real, intent(in) :: xvec(:), yvec(:), txt_scale, q_in, &
1398 : xmin, xmax, ymin, ymax, symbol_scale, str_scale
1399 : integer, intent(in) :: clr_index, symbol, str_clr
1400 : character (len = *), intent(in) :: str
1401 0 : real :: q, q0, q1, x, y, dy
1402 : integer :: nz, i, j, k
1403 : include 'formats'
1404 0 : q = max(0.0, min(1.0, q_in))
1405 0 : nz = s% nz
1406 0 : i = nz
1407 0 : dy = ymax - ymin
1408 0 : do k = 1, s% nz - 1
1409 0 : if (s% q(k) >= q .and. q > s% q(k + 1)) then
1410 : i = k; exit
1411 : end if
1412 : end do
1413 0 : j = i + 1
1414 0 : if (j >= nz) j = i
1415 0 : q0 = s% q(i)
1416 0 : q1 = s% q(j)
1417 0 : if ((q0 - q) * (q - q1) < 0) then
1418 0 : j = i - 1
1419 0 : q1 = s% q(j)
1420 : end if
1421 0 : x = find0(xvec(i), q0 - q, xvec(j), q1 - q)
1422 0 : if (x > xmax .or. x < xmin) return
1423 0 : y = find0(yvec(i), q0 - q, yvec(j), q1 - q)
1424 0 : if (y > ymax .or. y < ymin) return
1425 0 : call pgsave
1426 0 : call pgscf(1)
1427 0 : call pgslw(1)
1428 0 : call pgsci(clr_index)
1429 0 : call pgsch(symbol_scale * txt_scale)
1430 0 : call pgpt(1, x, y, symbol)
1431 0 : call pgsci(str_clr)
1432 0 : call pgsch(str_scale * txt_scale)
1433 0 : call pgptxt(x, y - 0.015 * dy, 0.0, 0.0, trim(str))
1434 0 : call pgunsa
1435 : end subroutine show_mass_point
1436 :
1437 :
1438 0 : real function find0(xx1, yy1, xx2, yy2)
1439 : real :: xx1, yy1, xx2, yy2
1440 0 : real :: a, b, xz
1441 : ! returns x where y is 0 on line connecting the points (xx1,yy1) and (xx2,yy2)
1442 0 : a = (xx1 * yy2) - (xx2 * yy1)
1443 0 : b = yy2 - yy1
1444 0 : if ((abs(a) >= abs(b) * 1e30) .and. ((yy1 >= 0 .and. yy2 <= 0) &
1445 : .or. (yy1 <= 0 .and. yy2 > 0))) then
1446 0 : xz = 0.5 * (xx1 + xx2)
1447 : else
1448 0 : xz = a / b
1449 : end if
1450 0 : find0 = xz
1451 0 : end function find0
1452 :
1453 :
1454 0 : subroutine show_box_pgstar(s, str1, str2)
1455 : type (star_info), pointer :: s
1456 : character (len = *), intent(in) :: str1, str2
1457 0 : real :: ch
1458 : integer :: lw
1459 0 : call pgqch(ch)
1460 0 : call pgqlw(lw)
1461 0 : call pgsch(s% pg% pgstar_num_scale * ch)
1462 0 : call pgslw(s% pg% pgstar_box_lw)
1463 0 : call pgbox(str1, 0.0, 0, str2, 0.0, 0)
1464 0 : call pgsch(ch)
1465 0 : call pgslw(lw)
1466 0 : end subroutine show_box_pgstar
1467 :
1468 :
1469 0 : subroutine show_grid_title_pgstar(s, title, pad)
1470 : type (star_info), pointer :: s
1471 : character (len = *), intent(in) :: title
1472 : real, intent(in) :: pad
1473 : optional pad
1474 0 : real :: ch, disp
1475 0 : if (.not. s% pg% pgstar_grid_show_title) return
1476 0 : if (len_trim(title) == 0) return
1477 0 : call pgqch(ch)
1478 0 : disp = s% pg% pgstar_grid_title_disp
1479 0 : if (present(pad)) disp = disp + pad
1480 : call do1_pgmtxt('T', disp, &
1481 : s% pg% pgstar_grid_title_coord, s% pg% pgstar_grid_title_fjust, title, &
1482 0 : s% pg% pgstar_grid_title_scale * ch, s% pg% pgstar_grid_title_lw)
1483 : end subroutine show_grid_title_pgstar
1484 :
1485 :
1486 0 : subroutine show_title_pgstar(s, title, pad)
1487 : type (star_info), pointer :: s
1488 : character (len = *), intent(in) :: title
1489 : real, intent(in) :: pad
1490 : optional pad
1491 0 : real :: ch, disp
1492 0 : if (.not. s% pg% pgstar_show_title) return
1493 0 : if (len_trim(title) == 0) return
1494 0 : call pgqch(ch)
1495 0 : disp = s% pg% pgstar_title_disp
1496 0 : if (present(pad)) disp = disp + pad
1497 : call do1_pgmtxt('T', disp, &
1498 : s% pg% pgstar_title_coord, s% pg% pgstar_title_fjust, title, &
1499 0 : s% pg% pgstar_title_scale * ch, s% pg% pgstar_title_lw)
1500 : end subroutine show_title_pgstar
1501 :
1502 :
1503 0 : subroutine show_title_label_pgmtxt_pgstar(&
1504 : s, coord, fjust, label, pad)
1505 : type (star_info), pointer :: s
1506 : character (len = *), intent(in) :: label
1507 : real, intent(in) :: pad, coord, fjust
1508 : optional pad
1509 : real :: disp
1510 0 : disp = s% pg% pgstar_title_disp
1511 0 : if (present(pad)) disp = disp + pad
1512 0 : call pgmtxt('T', disp, coord, fjust, label)
1513 0 : end subroutine show_title_label_pgmtxt_pgstar
1514 :
1515 :
1516 0 : subroutine show_xaxis_label_pgstar(s, label, pad)
1517 : type (star_info), pointer :: s
1518 : character (len = *), intent(in) :: label
1519 : real, intent(in) :: pad
1520 : optional pad
1521 0 : real :: ch, disp
1522 0 : call pgqch(ch)
1523 0 : disp = s% pg% pgstar_xaxis_label_disp
1524 0 : if (present(pad)) disp = disp + pad
1525 : call do1_pgmtxt('B', disp, 0.5, 0.5, label, &
1526 0 : s% pg% pgstar_xaxis_label_scale * ch, s% pg% pgstar_xaxis_label_lw)
1527 0 : end subroutine show_xaxis_label_pgstar
1528 :
1529 :
1530 0 : subroutine show_xaxis_label_pgmtxt_pgstar(&
1531 : s, coord, fjust, label, pad)
1532 : type (star_info), pointer :: s
1533 : character (len = *), intent(in) :: label
1534 : real, intent(in) :: pad, coord, fjust
1535 : optional pad
1536 : real :: disp
1537 0 : disp = s% pg% pgstar_xaxis_label_disp
1538 0 : if (present(pad)) disp = disp + pad
1539 0 : call pgmtxt('B', disp, coord, fjust, label)
1540 0 : end subroutine show_xaxis_label_pgmtxt_pgstar
1541 :
1542 :
1543 0 : subroutine show_left_yaxis_label_pgstar(s, label, pad)
1544 : type (star_info), pointer :: s
1545 : character (len = *), intent(in) :: label
1546 : real, intent(in) :: pad
1547 : optional pad
1548 0 : real :: ch, disp
1549 0 : call pgqch(ch)
1550 0 : disp = s% pg% pgstar_left_yaxis_label_disp
1551 0 : if (present(pad)) disp = disp + pad
1552 : call do1_pgmtxt('L', disp, 0.5, 0.5, label, &
1553 0 : s% pg% pgstar_left_yaxis_label_scale * ch, s% pg% pgstar_left_yaxis_label_lw)
1554 0 : end subroutine show_left_yaxis_label_pgstar
1555 :
1556 :
1557 0 : subroutine show_right_yaxis_label_pgstar(s, label, pad)
1558 : type (star_info), pointer :: s
1559 : character (len = *), intent(in) :: label
1560 : real, intent(in) :: pad
1561 : optional pad
1562 0 : real :: ch, disp
1563 0 : call pgqch(ch)
1564 0 : disp = s% pg% pgstar_right_yaxis_label_disp
1565 0 : if (present(pad)) disp = disp + pad
1566 : call do1_pgmtxt('R', disp, 0.5, 0.5, label, &
1567 0 : s% pg% pgstar_right_yaxis_label_scale * ch, s% pg% pgstar_right_yaxis_label_lw)
1568 0 : end subroutine show_right_yaxis_label_pgstar
1569 :
1570 :
1571 0 : subroutine show_left_yaxis_label_pgmtxt_pgstar(&
1572 : s, coord, fjust, label, pad)
1573 : type (star_info), pointer :: s
1574 : character (len = *), intent(in) :: label
1575 : real, intent(in) :: pad, coord, fjust
1576 : optional pad
1577 0 : real :: ch, disp
1578 0 : call pgqch(ch)
1579 0 : call pgsch(1.1 * ch)
1580 0 : disp = s% pg% pgstar_left_yaxis_label_disp
1581 0 : if (present(pad)) disp = disp + pad
1582 0 : call pgmtxt('L', disp, coord, fjust, label)
1583 0 : call pgsch(ch)
1584 0 : end subroutine show_left_yaxis_label_pgmtxt_pgstar
1585 :
1586 :
1587 0 : subroutine show_right_yaxis_label_pgmtxt_pgstar(&
1588 : s, coord, fjust, label, pad)
1589 : type (star_info), pointer :: s
1590 : character (len = *), intent(in) :: label
1591 : real, intent(in) :: pad, coord, fjust
1592 : optional pad
1593 0 : real :: ch, disp
1594 0 : call pgqch(ch)
1595 0 : call pgsch(1.1 * ch)
1596 0 : disp = s% pg% pgstar_right_yaxis_label_disp
1597 0 : if (present(pad)) disp = disp + pad
1598 0 : call pgmtxt('R', disp, coord, fjust, label)
1599 0 : call pgsch(ch)
1600 0 : end subroutine show_right_yaxis_label_pgmtxt_pgstar
1601 :
1602 :
1603 0 : subroutine show_model_number_pgstar(s)
1604 : type (star_info), pointer :: s
1605 : character (len = 32) :: str
1606 0 : real :: ch
1607 0 : if (.not. s% pg% pgstar_show_model_number) return
1608 0 : write(str, '(i9)') s% model_number
1609 0 : str = 'model ' // trim(adjustl(str))
1610 0 : call pgqch(ch)
1611 : call do1_pgmtxt('T', &
1612 : s% pg% pgstar_model_disp, s% pg% pgstar_model_coord, &
1613 : s% pg% pgstar_model_fjust, str, &
1614 0 : s% pg% pgstar_model_scale * ch, s% pg% pgstar_model_lw)
1615 : end subroutine show_model_number_pgstar
1616 :
1617 :
1618 0 : subroutine show_age_pgstar(s)
1619 : type (star_info), pointer :: s
1620 : character (len = 32) :: age_str, units_str
1621 : real(dp) :: age
1622 0 : real :: ch
1623 : integer :: len, i, j, iE
1624 0 : if (.not. s% pg% pgstar_show_age) return
1625 0 : age = s% star_age
1626 0 : if (s% pg% pgstar_show_age_in_seconds) then
1627 0 : age = age * secyer
1628 0 : units_str = 'secs'
1629 0 : else if (s% pg% pgstar_show_age_in_minutes) then
1630 0 : age = age * secyer / 60
1631 0 : units_str = 'mins'
1632 0 : else if (s% pg% pgstar_show_age_in_hours) then
1633 0 : age = age * secyer / (60 * 60)
1634 0 : units_str = 'hrs'
1635 0 : else if (s% pg% pgstar_show_age_in_days) then
1636 0 : age = age * secyer / secday
1637 0 : units_str = 'days'
1638 0 : else if (s% pg% pgstar_show_age_in_years) then
1639 : !age = age
1640 0 : units_str = 'yrs'
1641 0 : else if (s% pg% pgstar_show_log_age_in_years) then
1642 0 : age = log10(max(1d-99, age))
1643 0 : units_str = 'log yrs'
1644 0 : else if (age * secyer < 60) then
1645 0 : age = age * secyer
1646 0 : units_str = 'secs'
1647 0 : else if (age * secyer < 60 * 60) then
1648 0 : age = age * secyer / 60
1649 0 : units_str = 'mins'
1650 0 : else if (age * secyer < secday) then
1651 0 : age = age * secyer / (60 * 60)
1652 0 : units_str = 'hrs'
1653 0 : else if (age * secyer < secday * 500) then
1654 0 : age = age * secyer / secday
1655 0 : units_str = 'days'
1656 : else
1657 : !age = age
1658 0 : units_str = 'yrs'
1659 : end if
1660 0 : if (abs(age) > 1e-3 .and. abs(age) < 1e3) then
1661 0 : write(age_str, '(f14.6)') age
1662 : else
1663 0 : write(age_str, '(1pe14.6)') age
1664 0 : len = len_trim(age_str)
1665 0 : iE = 0
1666 0 : do i = 1, len
1667 0 : if (age_str(i:i) == 'E') then
1668 0 : iE = i
1669 0 : age_str(i:i) = 'e'
1670 0 : exit
1671 : end if
1672 : end do
1673 0 : if (iE > 0) then
1674 0 : i = iE + 1
1675 0 : if (age_str(i:i) == '+') then
1676 0 : do j = i, len - 1
1677 0 : age_str(j:j) = age_str(j + 1:j + 1)
1678 : end do
1679 0 : age_str(len:len) = ' '
1680 0 : len = len - 1
1681 : else
1682 0 : i = i + 1
1683 : end if
1684 0 : if (age_str(i:i) == '0') then
1685 0 : do j = i, len - 1
1686 0 : age_str(j:j) = age_str(j + 1:j + 1)
1687 : end do
1688 0 : age_str(len:len) = ' '
1689 0 : len = len - 1
1690 : end if
1691 : end if
1692 : end if
1693 0 : age_str = adjustl(age_str)
1694 0 : age_str = 'age ' // trim(age_str) // ' ' // trim(units_str)
1695 0 : call pgqch(ch)
1696 : call do1_pgmtxt('T', &
1697 : s% pg% pgstar_age_disp, s% pg% pgstar_age_coord, &
1698 : s% pg% pgstar_age_fjust, age_str, &
1699 0 : s% pg% pgstar_age_scale * ch, s% pg% pgstar_age_lw)
1700 : end subroutine show_age_pgstar
1701 :
1702 :
1703 0 : logical function read_values_from_file(fname, x_data, y_data, data_len)
1704 : character(len = *), intent(in) :: fname
1705 : real, allocatable, dimension(:) :: x_data, y_data
1706 : integer, intent(out) :: data_len
1707 : integer :: iounit, ierr, i
1708 : include 'formats'
1709 0 : read_values_from_file = .false.
1710 0 : ierr = 0
1711 0 : open(newunit = iounit, file = trim(fname), action = 'read', status = 'old', iostat = ierr)
1712 0 : if (ierr /= 0) then
1713 : !write(*, *) 'failed to open ' // trim(fname)
1714 : return
1715 : end if
1716 0 : read(iounit, *, iostat = ierr) data_len
1717 0 : if (ierr /= 0) then
1718 0 : write(*, *) 'failed to read num points on 1st line ' // trim(fname)
1719 0 : return
1720 : end if
1721 : !write(*,2) trim(fname) // ' data_len', data_len
1722 0 : allocate(x_data(data_len), y_data(data_len))
1723 0 : do i = 1, data_len
1724 0 : read(iounit, *, iostat = ierr) x_data(i), y_data(i)
1725 0 : if (ierr /= 0) then
1726 0 : write(*, *) 'failed to read data ' // trim(fname)
1727 0 : deallocate(x_data, y_data)
1728 0 : return
1729 : end if
1730 : end do
1731 0 : close(iounit)
1732 0 : read_values_from_file = .true.
1733 0 : end function read_values_from_file
1734 :
1735 :
1736 0 : subroutine show_pgstar_decorator(id, use_flag, pgstar_decorator, plot_num, ierr)
1737 : logical, intent(in) :: use_flag
1738 0 : real :: xmin, xmax, ymin, ymax
1739 : integer, intent(in) :: id, plot_num
1740 : integer, intent(inout) :: ierr
1741 : procedure(pgstar_decorator_interface), pointer :: pgstar_decorator
1742 :
1743 0 : if(use_flag)then
1744 0 : if(associated(pgstar_decorator))then
1745 0 : call pgsave
1746 0 : call PGQWIN(xmin, xmax, ymin, ymax)
1747 0 : call pgstar_decorator(id, xmin, xmax, ymin, ymax, plot_num, ierr)
1748 0 : call pgunsa
1749 0 : if(ierr/=0)then
1750 0 : write(*, *) "Error in pgstar_decorator"
1751 : end if
1752 : end if
1753 : end if
1754 :
1755 0 : end subroutine show_pgstar_decorator
1756 :
1757 : end module pgstar_support
1758 :
|
