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 :
|