LCOV - code coverage report
Current view: top level - star/private - pgstar_hr.f90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 0.0 % 111 0
Test Date: 2026-11-11 22:04:34 Functions: 0.0 % 4 0

            Line data    Source code
       1              : ! ***********************************************************************
       2              : !
       3              : !   Copyright (C) 2010  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_HR
      21              : 
      22              :       use star_private_def
      23              :       use const_def, only: dp
      24              :       use pgstar_support
      25              :       use star_pgstar
      26              : 
      27              :       implicit none
      28              : 
      29              : 
      30              :       contains
      31              : 
      32              : 
      33            0 :       subroutine HR_Plot(id, device_id, ierr)
      34              :          integer, intent(in) :: id, device_id
      35              :          integer, intent(out) :: ierr
      36              :          type (star_info), pointer :: s
      37              :          ierr = 0
      38            0 :          call get_star_ptr(id, s, ierr)
      39            0 :          if (ierr /= 0) return
      40              : 
      41            0 :          call pgslct(device_id)
      42            0 :          call pgbbuf()
      43            0 :          call pgeras()
      44              : 
      45              :          call do_HR_Plot(s, id, device_id, &
      46              :             s% pg% HR_xleft, s% pg% HR_xright, &
      47              :             s% pg% HR_ybot, s% pg% HR_ytop, .false., &
      48            0 :             s% pg% HR_title, s% pg% HR_txt_scale, ierr)
      49            0 :          if (ierr /= 0) return
      50              : 
      51            0 :          call pgebuf()
      52              : 
      53              :       end subroutine HR_Plot
      54              : 
      55              : 
      56            0 :       subroutine HR_decorate(id, ierr)
      57              :          use pgstar_colors
      58              :          integer, intent(in) :: id
      59              :          integer, intent(out) :: ierr
      60              :          type (star_info), pointer :: s
      61              :          real :: logT1, logT2, logL1, logL2
      62              :          ierr = 0
      63            0 :          call get_star_ptr(id, s, ierr)
      64            0 :          if (ierr /= 0) return
      65            0 :          if (s% pg% show_HR_radius_lines) call HR_radius_lines()
      66            0 :          if (s% pg% show_HR_Mira_instability_region) then
      67            0 :             logL1 = 2.8
      68            0 :             logL2 = 4.0
      69            0 :             logT2 = 3.7
      70            0 :             logT1 = 3.45
      71              : 
      72            0 :             call pgsls(Line_Type_Solid)
      73            0 :             call pgslw(s% pg% pgstar_lw)
      74            0 :             call pgsci(clr_Blue)
      75            0 :             call pgmove(logT1, logL1)
      76            0 :             call pgdraw(logT1, logL2)
      77            0 :             call pgdraw(logT2, logL2)
      78            0 :             call pgdraw(logT2, logL1)
      79            0 :             call pgdraw(logT1, logL1)
      80              :          end if
      81            0 :          if (s% pg% show_HR_classical_instability_strip) then
      82              : 
      83            0 :             call pgsls(Line_Type_Solid)
      84            0 :             call pgslw(s% pg% pgstar_lw)
      85              : 
      86              :             ! approximate edges
      87              : 
      88              :             ! blue edge
      89            0 :             logT1 = 3.70
      90            0 :             logL1 = 5.5
      91              : 
      92            0 :             logT2 = 3.93
      93            0 :             logL2 = 1.0
      94              : 
      95            0 :             call pgsci(clr_Blue)
      96            0 :             call pgmove(logT1, logL1)
      97            0 :             call pgdraw(logT2, logL2)
      98              : 
      99              :             ! red edge
     100            0 :             logT1 = 3.60
     101            0 :             logL1 = 5.5
     102              : 
     103            0 :             logT2 = 3.83
     104            0 :             logL2 = 1.0
     105              : 
     106            0 :             call pgsci(clr_FireBrick)
     107            0 :             call pgmove(logT1, logL1)
     108            0 :             call pgdraw(logT2, logL2)
     109              :          end if
     110            0 :          if (s% pg% show_HR_WD_instabilities) then
     111              :             ! from Winget & Kepler, Annu. Rev. Astron. Astrophys., 2008, Fig 3.
     112            0 :             call pgsave
     113            0 :             call pgsls(Line_Type_Solid)
     114            0 :             call pgslw(s% pg% pgstar_lw)
     115            0 :             call pgsci(clr_Silver)
     116            0 :             call pgmove(5.1, 4.3)  ! DOV
     117            0 :             call pgdraw(4.8, 0.3)
     118            0 :             call pgdraw(4.95, 0.9)
     119            0 :             call pgdraw(5.38, 4.2)
     120            0 :             call pgdraw(5.1, 4.3)
     121            0 :             call pgmove(4.42, -0.7)  ! DBV
     122            0 :             call pgdraw(4.34, -1.2)
     123            0 :             call pgdraw(4.38, -1.7)
     124            0 :             call pgdraw(4.45, -1.4)
     125            0 :             call pgdraw(4.42, -0.7)
     126            0 :             call pgmove(4.03, -2.4)  ! DAV
     127            0 :             call pgdraw(4.03, -3.4)
     128            0 :             call pgdraw(4.1, -3.2)
     129            0 :             call pgdraw(4.1, -2.2)
     130            0 :             call pgdraw(4.03, -2.4)
     131            0 :             call pgsch(s% pg% HR_txt_scale*0.7)
     132            0 :             call pgslw(1)
     133            0 :             call pgptxt(5.1 - 0.05, 4.3, 0.0, 0.0, 'DOV')
     134            0 :             call pgptxt(4.42 - 0.05, -0.7, 0.0, 0.0, 'DBV')
     135            0 :             call pgptxt(4.03 - 0.05, -2.4, 0.0, 0.0, 'DAV')
     136            0 :             call pgunsa
     137              :          end if
     138              :       end subroutine HR_decorate
     139              : 
     140              : 
     141            0 :       subroutine HR_radius_lines()
     142              :          use const_def, only: pi, boltz_sigma, Rsun, Lsun
     143              :          use math_lib, only: log10
     144              :          use pgstar_colors, only: clr_Gray
     145              :          character(len=48) :: label
     146              :          character(len=16) :: radius
     147              :          real :: xmin, xmax, ymin, ymax
     148              :          real :: x(2), y(2), xlo, xhi, ylo, yhi, dx, dy
     149              :          real :: intercept, label_lo, label_hi, label_x, char_height
     150              :          real :: solar_intercept
     151              :          integer :: i, imin, imax
     152              : 
     153            0 :          call pgqwin(xmin, xmax, ymin, ymax)
     154              :          ! HR temperature axes run from hot to cool; handle either ordering.
     155            0 :          xlo = min(xmin, xmax)
     156            0 :          xhi = max(xmin, xmax)
     157            0 :          ylo = min(ymin, ymax)
     158            0 :          yhi = max(ymin, ymax)
     159            0 :          dx = xhi - xlo
     160            0 :          dy = yhi - ylo
     161            0 :          if (dx <= 0.0 .or. dy <= 0.0) return
     162              : 
     163              :          ! save the current drawing settings -- we'll restore them later
     164              :          ! so other decorations don't inherit this style
     165            0 :          call pgsave
     166              : 
     167            0 :          call pgsci(clr_Gray)
     168            0 :          call pgsls(Line_Type_Dash)
     169            0 :          call pgslw(1)
     170            0 :          call pgqch(char_height)
     171              :          ! do_Hist_Track leaves the multiplied the character size by 2.8, so undo
     172              :          ! that here.
     173            0 :          call pgsch(0.7*char_height/2.8)
     174              : 
     175              :          ! log10(L/Lsun) = 4*log10(Teff) + 2*log10(R/Rsun) + solar_intercept.
     176            0 :          solar_intercept = real(log10(4.0_dp*pi*boltz_sigma*pow2(Rsun)/Lsun))  ! cast to real b/c pgstar uses single precision
     177              : 
     178              :          ! figure out which radii (powers of 10) are visible
     179              :          ! smallest visible line is at highest Teff and lowest L (xhi, ylo)
     180              :          ! largest visible line is at lowest Teff and highest L (xlo, yhi)
     181            0 :          imin = ceiling(0.5*(ylo - 4.0*xhi - solar_intercept))
     182            0 :          imax = floor(0.5*(yhi - 4.0*xlo - solar_intercept))
     183              : 
     184            0 :          do i = imin, imax
     185              :             ! get the intercept for the current radius
     186            0 :             intercept = solar_intercept + 2.0*i
     187              : 
     188              :             ! Clip the line to the current axes before drawing it.
     189            0 :             x(1) = max(xlo, (ylo-intercept)/4.0)
     190            0 :             x(2) = min(xhi, (yhi-intercept)/4.0)
     191            0 :             if (x(1) >= x(2)) cycle
     192            0 :             y = 4.0*x + intercept
     193            0 :             call pgline(2, x, y)
     194              : 
     195              :             ! Label an interior segment so labels follow changing axis limits.
     196            0 :             label_lo = max(xlo + 0.12*dx, (ylo + 0.08*dy-intercept)/4.0)
     197            0 :             label_hi = min(xhi - 0.12*dx, (yhi - 0.12*dy-intercept)/4.0)
     198              : 
     199            0 :             if (label_lo >= label_hi) cycle
     200              : 
     201            0 :             label_x = 0.5*(label_lo + label_hi)
     202              : 
     203            0 :             if (i >= 0 .and. i <= 3) then
     204            0 :                write(radius, '(i0)') 10**i                    ! 1, 10, 100, 1000
     205            0 :                label = trim(radius) // ' R\d\(2281)\u'
     206            0 :             else if (i < 0 .and. i >= -3) then
     207            0 :                radius = '0.' // repeat('0', -i-1) // '1'      ! 0.1, 0.01, 0.001
     208            0 :                label = trim(radius) // ' R\d\(2281)\u'
     209              :             else
     210            0 :                write(radius, '(i0)') i                        ! i0 includes the minus sign
     211            0 :                label = '10\u' // trim(radius) // '\d R\d\(2281)\u'
     212              :             end if
     213              : 
     214              :             call pgptxt(label_x, 4.0*label_x + intercept + 0.02*dy, &
     215            0 :                0.0, 0.5, trim(label))
     216              :          end do
     217              : 
     218              :          ! restore the drawing settings
     219            0 :          call pgunsa
     220              : 
     221              :       end subroutine HR_radius_lines
     222              : 
     223              : 
     224            0 :       subroutine do_HR_Plot(s, id, device_id, &
     225              :             xleft, xright, ybot, ytop, subplot, &
     226              :             title, txt_scale, ierr)
     227              :          use pgstar_hist_track, only: do_Hist_Track
     228              :          type (star_info), pointer :: s
     229              :          integer, intent(in) :: id, device_id
     230              :          real, intent(in) :: xleft, xright, ybot, ytop, txt_scale
     231              :          logical, intent(in) :: subplot
     232              :          character (len=*), intent(in) :: title
     233              :          integer, intent(out) :: ierr
     234              :          logical, parameter :: &
     235              :             reverse_xaxis = .true., reverse_yaxis = .false.
     236              :          ierr = 0
     237              :          call do_Hist_Track(s, id, device_id, &
     238              :             xleft, xright, ybot, ytop, subplot, title, txt_scale, &
     239              :             'log_Teff', 'log_L', &
     240              :             'log Teff', 'log L/L\d\(2281)', &
     241              :             s% pg% HR_logT_min, s% pg% HR_logT_max, &
     242              :             s% pg% HR_logT_margin, s% pg% HR_dlogT_min, &
     243              :             s% pg% HR_logL_min, s% pg% HR_logL_max, &
     244              :             s% pg% HR_logL_margin, s% pg% HR_dlogL_min, &
     245              :             s% pg% HR_step_min, s% pg% HR_step_max, &
     246              :             reverse_xaxis, reverse_yaxis, .false., .false., &
     247              :             s% pg% show_HR_target_box, s% pg% HR_target_n_sigma, &
     248              :             s% pg% HR_target_logT, s% pg% HR_target_logL, &
     249              :             s% pg% HR_target_logT_sigma, s% pg% HR_target_logL_sigma, &
     250              :             s% pg% show_HR_annotation1, &
     251              :             s% pg% show_HR_annotation2, &
     252              :             s% pg% show_HR_annotation3, &
     253              :             s% pg% HR_fname, &
     254              :             s% pg% HR_use_decorator, &
     255              :             s% pg% HR_pgstar_decorator, &
     256            0 :             HR_decorate, ierr)
     257            0 :       end subroutine do_HR_Plot
     258              : 
     259              : 
     260              :       end module pgstar_HR
     261              : 
        

Generated by: LCOV version 2.0-1