Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-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_summary_burn
21 :
22 : use star_private_def
23 : use const_def, only: dp, ln10
24 : use pgstar_support
25 : use star_pgstar
26 :
27 : implicit none
28 :
29 : contains
30 :
31 0 : subroutine summary_burn_plot(id, device_id, ierr)
32 : integer, intent(in) :: id, device_id
33 : integer, intent(out) :: ierr
34 :
35 : type (star_info), pointer :: s
36 : ierr = 0
37 0 : call get_star_ptr(id, s, ierr)
38 0 : if (ierr /= 0) return
39 :
40 0 : call pgslct(device_id)
41 0 : call pgbbuf()
42 0 : call pgeras()
43 :
44 : call do_summary_burn_plot(s, id, device_id, &
45 : s% pg% Summary_Burn_xleft, s% pg% Summary_Burn_xright, &
46 : s% pg% Summary_Burn_ybot, s% pg% Summary_Burn_ytop, .false., &
47 0 : s% pg% Summary_Burn_title, s% pg% Summary_Burn_txt_scale, ierr)
48 :
49 0 : call pgebuf()
50 :
51 : end subroutine summary_burn_plot
52 :
53 :
54 0 : subroutine do_summary_burn_plot(s, id, device_id, &
55 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, ierr)
56 : use utils_lib
57 : use chem_def
58 : use net_def
59 : use const_def, only: Msun, Rsun
60 : type (star_info), pointer :: s
61 : integer, intent(in) :: id, device_id
62 : real, intent(in) :: winxmin, winxmax, winymin, winymax
63 : logical, intent(in) :: subplot
64 : character (len=*), intent(in) :: title
65 : real, intent(in) :: txt_scale
66 : integer, intent(out) :: ierr
67 :
68 : character (len=strlen) :: xaxis_name, str
69 : logical :: xaxis_reversed
70 0 : real, allocatable, dimension(:) :: xvec, yvec, yvec2, yvec3
71 0 : real :: xmin, xmax, xleft, xright, dx, windy, dy, &
72 0 : ymin, ymax, xaxis_min, xaxis_max, xmargin, &
73 0 : legend_xmin, legend_xmax, legend_ymin, legend_ymax
74 : integer :: lw, lw_sav, grid_min, grid_max, npts, nz
75 : integer :: docat(num_categories), eps_k(num_categories), num_cat
76 0 : real :: xpos_nuc(num_categories)
77 0 : real :: xnuc_cat(num_categories), coord
78 0 : real(dp) :: eps_max, eps, maxv
79 :
80 : include 'formats'
81 :
82 0 : ierr = 0
83 0 : xaxis_name = s% pg% Summary_Burn_xaxis_name
84 0 : xaxis_min = s% pg% Summary_Burn_xmin
85 0 : xaxis_max = s% pg% Summary_Burn_xmax
86 0 : xaxis_reversed = s% pg% Summary_Burn_xaxis_reversed
87 :
88 0 : nz = s% nz
89 0 : windy = winymax - winymin
90 :
91 0 : legend_xmin = winxmax - 0.01
92 0 : legend_xmax = 0.99
93 0 : legend_ymin = winymin
94 0 : legend_ymax = winymax
95 :
96 0 : allocate (xvec(nz), yvec(nz), yvec2(nz), yvec3(nz))
97 :
98 0 : xmargin = 0
99 : call set_xaxis_bounds( &
100 : s, xaxis_name, xaxis_min, xaxis_max, xaxis_reversed, xmargin, &
101 : xvec, xmin, xmax, xleft, xright, dx, &
102 0 : grid_min, grid_max, npts, ierr)
103 :
104 0 : if (ierr == 0) then
105 0 : call pgsave
106 0 : call pgsch(txt_scale)
107 0 : call plot(ierr)
108 0 : call pgunsa
109 : end if
110 :
111 0 : deallocate(xvec, yvec, yvec2, yvec3)
112 :
113 :
114 : contains
115 :
116 :
117 0 : subroutine plot(ierr)
118 0 : use rates_def
119 : use pgstar_colors
120 : integer, intent(out) :: ierr
121 :
122 : integer :: j, ii, jj, i, cnt, k
123 : logical, parameter :: dbg = .false.
124 0 : real :: ybot
125 :
126 : include 'formats'
127 :
128 0 : ymax = 1.02
129 0 : ymin = 0.0
130 0 : ybot = -0.02
131 :
132 0 : lw = s% pg% pgstar_lw
133 0 : call pgqlw(lw_sav)
134 :
135 0 : call pgsvp(winxmin, winxmax, winymin, winymax)
136 0 : if (.not. subplot) then
137 0 : call show_model_number_pgstar(s)
138 0 : call show_age_pgstar(s)
139 : end if
140 0 : call show_title_pgstar(s, title, s% pg% Summary_Burn_title_shift)
141 :
142 : ! logT
143 0 : do k=grid_min,grid_max
144 0 : yvec(k) = s% lnT(k)/ln10
145 : end do
146 0 : ymax = maxval(yvec(grid_min:grid_max))
147 0 : ymin = minval(yvec(grid_min:grid_max))
148 0 : dy = ymax-ymin
149 0 : ymax = ymax + 0.1*dy
150 0 : ymin = ymin - 0.1*dy
151 :
152 0 : call pgswin(xleft, xright, ymin+ybot, ymax)
153 :
154 0 : call pgscf(1)
155 0 : call pgsci(clr_Foreground)
156 0 : call pgsch(txt_scale)
157 0 : call pgbox('BCNST',0.0,0,'CMSTV',0.0,0)
158 0 : call pgsci(clr_Goldenrod)
159 0 : call pgmtxt('R',5.0,0.5,0.5,'log T')
160 0 : call pgslw(lw)
161 0 : call pgline(npts, xvec(grid_min:grid_max), yvec(grid_min:grid_max))
162 0 : call pgslw(lw_sav)
163 :
164 0 : do k=grid_min,grid_max
165 0 : yvec(k) = safe_log10(s% non_nuc_neu(k) + s% eps_nuc_neu_total(k))
166 0 : yvec2(k) = s% lnd(k)/ln10
167 0 : yvec3(k) = safe_log10(abs(s% eps_nuc(k)))
168 : end do
169 :
170 0 : ymax = maxval(yvec(grid_min:grid_max))
171 0 : ymax = max(ymax,maxval(yvec2(grid_min:grid_max)))
172 0 : ymax = max(ymax,maxval(yvec3(grid_min:grid_max)))
173 0 : if (ymax <= 10) then
174 0 : ymax = 11.1
175 : else
176 0 : ymax = ymax*1.1
177 : end if
178 :
179 0 : ymin = minval(yvec(grid_min:grid_min))
180 0 : ymin = min(ymin,minval(yvec2(grid_min:grid_min)))
181 0 : ymin = min(ymin,minval(yvec3(grid_min:grid_min)))
182 0 : ymin = max(-6.6, ymin)
183 :
184 0 : ymax = max(ymax, ymin+1)
185 0 : dy = ymax-ymin
186 0 : ymin = ymin-dy*0.1
187 :
188 0 : call pgswin(xleft, xright, ymin+ybot, ymax)
189 0 : call pgscf(1)
190 0 : call pgsci(clr_Foreground)
191 0 : call pgsch(txt_scale)
192 0 : call pgbox('',0.0,0,'BNSTV',0.0,0)
193 :
194 0 : call pgsci(clr_Lilac)
195 0 : call pgmtxt('L',4.9,0.5,0.5,'log \gr')
196 0 : call pgslw(lw)
197 0 : call pgline(npts, xvec(grid_min:grid_max), yvec2(grid_min:grid_max))
198 0 : call pgslw(lw_sav)
199 :
200 0 : call pgsci(clr_LightSteelBlue)
201 0 : call pgsch(txt_scale)
202 0 : call pgmtxt('L',3.7,0.5,0.5,'log \ge\dnuc\u')
203 0 : call pgqlw(lw)
204 0 : call pgslw(8)
205 0 : call pgline(npts, xvec(grid_min:grid_max), yvec3(grid_min:grid_max))
206 0 : call pgslw(lw_sav)
207 :
208 0 : call pgsci(clr_LightSkyGreen)
209 0 : call pgsch(txt_scale)
210 0 : call pgmtxt('L',2.5,0.5,0.5,'log \ge\d\gn\u burn+thermal')
211 0 : call pgqlw(lw)
212 0 : call pgslw(8)
213 0 : call pgline(npts, xvec(grid_min:grid_max), yvec(grid_min:grid_max))
214 0 : call pgslw(lw_sav)
215 :
216 : ! label peaks of eps_nuc
217 : ! (stack the labels if they are too close in the x-direction)
218 0 : num_cat = 0
219 0 : do i = 1, num_categories
220 0 : if (i == iphoto) cycle
221 0 : maxv = safe_log10(maxval(s% eps_nuc_categories(i,1:nz)))
222 0 : if (maxv > -1) then
223 0 : num_cat = num_cat + 1
224 0 : docat(num_cat) = i
225 0 : xnuc_cat(num_cat) = maxv
226 : end if
227 : end do
228 0 : call pgsci(clr_Foreground)
229 0 : call pgsch(txt_scale*0.8)
230 0 : do ii = 1, num_cat
231 0 : eps_max = -100; i = 0
232 0 : do jj = 1, num_cat
233 0 : if (xnuc_cat(jj) < -1) cycle
234 0 : if (xnuc_cat(jj) > eps_max) then
235 0 : eps_max = xnuc_cat(jj)
236 0 : i = jj
237 : end if
238 : end do
239 0 : if (i == 0) exit
240 : if (dbg) write(*,2) 'place ' // category_name(docat(i)), i, eps_max
241 0 : xnuc_cat(i) = -1e10 ! mark as done
242 0 : eps_max = -100
243 0 : eps_k(i) = 0
244 0 : do k = 1, nz ! if limit this to grid_min:grid_max, locations jump around too much
245 0 : eps = s% eps_nuc_categories(docat(i),k)
246 0 : if(eps > eps_max) then
247 0 : eps_max = eps
248 0 : eps_k(i) = k
249 : end if
250 : end do
251 0 : if(eps_k(i) > 0) then
252 0 : k = eps_k(i)
253 0 : xpos_nuc(i) = xvec(k)
254 0 : if (xpos_nuc(i) < xmin .or. xpos_nuc(i) > xmax) cycle
255 : cnt = 0
256 0 : do j = 1, num_cat ! compare location to ones already placed
257 0 : if (j == i) cycle
258 0 : if (xnuc_cat(j) > -1) cycle ! haven't done this one yet
259 0 : if (abs(xpos_nuc(i) - xpos_nuc(j)) < 0.1*dx) then
260 0 : cnt = cnt + 1
261 : if (dbg) write(*,*) 'conflicts with ' // category_name(docat(j))
262 : end if
263 : end do
264 0 : if (cnt < 3) then ! only show 3 max
265 0 : str = category_name(docat(i))
266 0 : if (str(1:5) == 'burn_') then
267 0 : str = str(6:len_trim(str))
268 0 : else if (str == 'tri_alpha') then
269 0 : str = '3a'
270 0 : else if (str == 'c12_c12') then
271 0 : str = 'c+c'
272 0 : else if (str == 'c12_o16') then
273 0 : str = 'c+o'
274 0 : else if (str == 'o16_o16') then
275 0 : str = 'o+o'
276 : end if
277 :
278 0 : coord = (xpos_nuc(i) - xmin)/(xmax - xmin)
279 0 : if (xaxis_reversed) coord = 1.0 - coord
280 : call show_title_label_pgmtxt_pgstar( &
281 0 : s, coord, 0.5, trim(str), 1.03*(cnt+1))
282 : !call pgptxt(xpos_nuc(i), ylb, 0.0, 0.5, trim(str))
283 : end if
284 : end if
285 : end do
286 :
287 0 : call pgsci(clr_Foreground)
288 : ierr = 0
289 0 : call show_xaxis_name(s,xaxis_name,ierr)
290 0 : if (ierr == 0) then ! show mix regions at bottom of plot
291 0 : call pgslw(10)
292 : call show_mix_regions_on_xaxis( &
293 0 : s,ymin+ybot,ymax,grid_min,grid_max,xvec)
294 : end if
295 :
296 : call show_pgstar_decorator(s%id, s% pg% summary_burn_use_decorator, &
297 0 : s% pg% summary_burn_pgstar_decorator, 0, ierr)
298 :
299 0 : end subroutine plot
300 :
301 : end subroutine do_summary_burn_plot
302 :
303 : end module pgstar_summary_burn
|
