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_abundance
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 : contains
30 :
31 0 : subroutine abundance_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 :
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_abundance_plot(s, id, device_id, &
46 : s% pg% Abundance_xleft, s% pg% Abundance_xright, &
47 : s% pg% Abundance_ybot, s% pg% Abundance_ytop, .false., &
48 0 : s% pg% Abundance_title, s% pg% Abundance_txt_scale, ierr)
49 :
50 0 : call pgebuf()
51 :
52 : end subroutine abundance_plot
53 :
54 :
55 0 : subroutine do_abundance_plot(s, id, device_id, &
56 : winxmin, winxmax, winymin, winymax, subplot, &
57 : title, txt_scale, ierr)
58 : type (star_info), pointer :: s
59 : integer, intent(in) :: id, device_id
60 : real, intent(in) :: winxmin, winxmax, winymin, winymax
61 : logical, intent(in) :: subplot
62 : character (len=*), intent(in) :: title
63 : real, intent(in) :: txt_scale
64 : integer, intent(out) :: ierr
65 : call do_abundance_panel(s, id, device_id, &
66 : winxmin, winxmax, winymin, winymax, subplot, &
67 : title, txt_scale, s% pg% Abundance_xaxis_name, &
68 : s% pg% Abundance_xmin, s% pg% Abundance_xmax, s% pg% Abundance_xaxis_reversed, &
69 : s% pg% abundance_log_mass_frac_min, s% pg% abundance_log_mass_frac_max, &
70 0 : .false., .true., ierr)
71 0 : end subroutine do_abundance_plot
72 :
73 :
74 0 : subroutine do_abundance_panel(s, id, device_id, &
75 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, &
76 : xaxis_name, xaxis_min, xaxis_max, xaxis_reversed, ymin_in, ymax_in, &
77 : panel_flag, xaxis_numeric_labels_flag, ierr)
78 : use utils_lib
79 : use chem_def
80 : use net_def
81 : use const_def, only: Msun, Rsun
82 : use pgstar_colors
83 :
84 : type (star_info), pointer :: s
85 : integer, intent(in) :: id, device_id
86 : real, intent(in) :: &
87 : winxmin, winxmax, winymin, winymax, xaxis_min, xaxis_max, ymin_in, ymax_in
88 : character (len=*), intent(in) :: title, xaxis_name
89 : real, intent(in) :: txt_scale
90 : logical, intent(in) :: subplot, &
91 : xaxis_reversed, panel_flag, xaxis_numeric_labels_flag
92 : integer, intent(out) :: ierr
93 :
94 : real, allocatable, dimension(:) :: xvec, yvec
95 0 : real :: xmin, xmax, xleft, xright, dx, dylbl, chScale, windy, xmargin, &
96 : ymin, ymax, legend_xmin, legend_xmax, legend_ymin, legend_ymax
97 : integer :: lw, lw_sav, grid_min, grid_max, npts, nz
98 : integer, parameter :: num_colors = 14
99 : integer :: colors(num_colors)
100 : integer, parameter :: max_num_labels = 30
101 : integer :: num_labels, iloc_abundance_label(max_num_labels)
102 0 : real :: xloc_abundance_label(max_num_labels)
103 :
104 : include 'formats'
105 0 : ierr = 0
106 0 : nz = s% nz
107 :
108 : colors(:) = [ &
109 : clr_Gold, clr_LightSkyBlue, clr_Crimson, clr_Goldenrod, clr_MediumSlateBlue, &
110 : clr_Coral, clr_LightSkyGreen, clr_DarkGray, clr_Lilac, &
111 0 : clr_Tan, clr_IndianRed, clr_Teal, clr_Silver, clr_BrightBlue ]
112 :
113 0 : chScale = txt_scale
114 :
115 0 : ymin = ymin_in
116 0 : ymax = ymax_in
117 0 : if (abs(ymin+101.0) < 1e-6) ymin = s% pg% abundance_log_mass_frac_min
118 0 : if (abs(ymax+101.0) < 1e-6) ymax = s% pg% abundance_log_mass_frac_max
119 :
120 0 : windy = winymax - winymin
121 :
122 0 : legend_xmin = winxmax
123 0 : legend_xmax = min(1.0, winxmax + (winxmax - winxmin)/5)
124 :
125 0 : legend_ymin = winymin
126 0 : legend_ymax = winymax
127 :
128 0 : allocate(xvec(nz), yvec(nz))
129 :
130 0 : xmargin = 0
131 : call set_xaxis_bounds( &
132 : s, xaxis_name, xaxis_min, xaxis_max, xaxis_reversed, &
133 : xmargin, xvec, xmin, xmax, xleft, xright, dx, &
134 0 : grid_min, grid_max, npts, ierr)
135 0 : if (ierr == 0) call plot(ierr)
136 :
137 0 : deallocate(xvec, yvec)
138 :
139 : contains
140 :
141 :
142 0 : subroutine plot(ierr)
143 0 : use rates_def
144 : integer, intent(out) :: ierr
145 :
146 : integer :: i, k
147 : logical, parameter :: dbg = .false.
148 0 : real(dp) :: photosphere_logxm
149 0 : real :: lgz, x, ybot
150 :
151 : include 'formats'
152 0 : ierr = 0
153 :
154 0 : if (ymin > 0) then
155 0 : ymin = -5.1
156 0 : lgz = log10(s% initial_z + 1e-9) - 1
157 0 : if (lgz-1 < ymin) ymin = lgz
158 : end if
159 0 : if (ymax >= 100) then
160 0 : if (ymin < -8) then
161 0 : ymax = 0.5
162 : else
163 0 : ymax = 0.25
164 : end if
165 : end if
166 :
167 0 : num_labels = max(0,min(max_num_labels, s% pg% num_abundance_line_labels))
168 :
169 0 : iloc_abundance_label = -HUGE(grid_min)
170 0 : xloc_abundance_label = -HUGE(grid_min)
171 0 : do i=1,num_labels
172 0 : x = xmin + (i-0.5)*dx/num_labels
173 0 : do k=2,nz
174 0 : if ((xvec(k-1)-x)*(x-xvec(k)) >= 0) then
175 0 : iloc_abundance_label(i) = k
176 0 : xloc_abundance_label(i) = x
177 0 : exit
178 : end if
179 : end do
180 : end do
181 :
182 0 : dylbl = (ymax - ymin)*0.015
183 :
184 0 : lw = s% pg% pgstar_lw
185 0 : call pgqlw(lw_sav)
186 :
187 0 : call pgsave
188 0 : call pgsch(txt_scale)
189 0 : call pgsvp(legend_xmin, legend_xmax, legend_ymin, legend_ymax)
190 0 : call pgswin(0.0, 1.0, ymin, ymax)
191 0 : call do_all(.true.)
192 0 : call pgunsa
193 :
194 0 : call pgsave
195 0 : call pgsch(txt_scale)
196 :
197 0 : call pgsvp(winxmin, winxmax, winymin, winymax)
198 0 : if (.not. panel_flag) then
199 0 : if (.not. subplot) then
200 0 : call show_model_number_pgstar(s)
201 0 : call show_age_pgstar(s)
202 : end if
203 0 : call show_title_pgstar(s, title)
204 0 : call pgsci(clr_Foreground)
205 0 : call show_xaxis_name(s,xaxis_name,ierr)
206 0 : if (ierr /= 0) return
207 : end if
208 :
209 0 : ybot = -0.05
210 0 : call pgswin(xleft, xright, ymin+ybot, ymax)
211 0 : call pgscf(1)
212 0 : call pgsci(clr_Foreground)
213 0 : if (xaxis_numeric_labels_flag) then
214 0 : call show_box_pgstar(s,'BCNST','BCNSTV')
215 : else
216 0 : call show_box_pgstar(s,'BCST','BCNSTV')
217 : end if
218 0 : call show_left_yaxis_label_pgstar(s, 'log mass fraction')
219 :
220 0 : call pgsave
221 0 : call pgsch(txt_scale*1.05)
222 0 : call do_all(.false.)
223 0 : call pgunsa
224 :
225 0 : if (.not. panel_flag) then ! show mix regions at bottom of plot
226 0 : call pgslw(10)
227 : call show_mix_regions_on_xaxis( &
228 0 : s,ymin+ybot,ymax,grid_min,grid_max,xvec)
229 : end if
230 :
231 0 : call pgunsa
232 :
233 0 : if (s% pg% Abundance_show_photosphere_location .and. &
234 : (xaxis_name == 'mass' .or. &
235 : xaxis_name == 'logxm' .or. &
236 : xaxis_name == 'radius')) then
237 0 : call pgsave
238 0 : call pgsvp(winxmin, winxmax, winymin, winymax)
239 0 : call pgswin(0.0, 1.0, 0.0, 1.0)
240 0 : call pgsci(clr_Gray)
241 0 : call pgsls(Line_Type_Dash)
242 0 : call pgslw(5)
243 0 : if (xaxis_name == 'radius') then
244 0 : dx = (s% photosphere_r - xmin)/(xmax - xmin)
245 0 : else if (xaxis_name == 'radius_cm') then
246 0 : dx = (s% photosphere_r*Rsun - xmin)/(xmax - xmin)
247 0 : else if (xaxis_name == 'mass') then
248 0 : dx = (s% photosphere_m - xmin)/(xmax - xmin)
249 : else
250 0 : if (s% star_mass > s% photosphere_m) then
251 0 : photosphere_logxm = log10(s% star_mass - s% photosphere_m)
252 : else
253 0 : photosphere_logxm = -99
254 : end if
255 : !dx = (xmax - photosphere_logxm)/(xmin - xmax)
256 0 : dx = 1 - (photosphere_logxm - xmin)/(xmax - xmin)
257 0 : write(*,2) 'photosphere_xm xmin photo_x xmax dx', s% model_number, &
258 0 : s% star_mass - s% photosphere_m, &
259 0 : xmin, photosphere_logxm, xmax, dx
260 : end if
261 0 : call pgmove(dx, 0.0)
262 0 : call pgdraw(dx, 1.0)
263 0 : call pgunsa
264 : end if
265 :
266 0 : call show_pgstar_decorator(s%id,s% pg% Abundance_use_decorator,s% pg% Abundance_pgstar_decorator,0, ierr)
267 :
268 0 : end subroutine plot
269 :
270 :
271 0 : subroutine do_all(legend_flag)
272 : logical, intent(in) :: legend_flag
273 : integer :: cnt, num_to_show, i, j, jmax
274 0 : real(dp) :: max_abund(s% species)
275 : include 'formats'
276 0 : cnt = 0
277 0 : num_to_show = s% pg% Abundance_num_isos_to_show
278 0 : do j=1, s% species
279 0 : max_abund(j) = maxval(s% xa(j,grid_min:grid_max))
280 : end do
281 0 : if (num_to_show < 0) then ! show as many as fit
282 0 : if (legend_flag) then
283 0 : jmax = min(s% species, max_num_labels)
284 : else
285 : jmax = s% species
286 : end if
287 0 : do j = 1, jmax
288 0 : i = maxloc(max_abund(:),dim=1)
289 0 : cnt = do1(cnt, chem_isos% name(s% chem_id(i)), legend_flag)
290 0 : max_abund(i) = -1
291 : end do
292 : else
293 0 : do i = 1, num_to_show
294 0 : cnt = do1(cnt, s% pg% Abundance_which_isos_to_show(i), legend_flag)
295 : end do
296 : end if
297 0 : end subroutine do_all
298 :
299 :
300 0 : integer function do1(cnt, str, legend_flag)
301 : use chem_lib, only: chem_get_iso_id
302 : integer, intent(in) :: cnt
303 : character (len=*), intent(in) :: str
304 : logical, intent(in) :: legend_flag
305 : integer :: i, cid, k
306 : include 'formats'
307 0 : do1 = cnt
308 0 : if (len_trim(str) == 0) return
309 0 : cid = chem_get_iso_id(str)
310 0 : if (cid <= 0) return
311 0 : i = s% net_iso(cid)
312 0 : if (i == 0) return
313 0 : do k=1,nz
314 0 : yvec(k) = safe_log10(s% xa(i,k))
315 : end do
316 0 : if (s% pg% abundance_log_mass_frac_min < 0) then
317 0 : if (maxval(yvec(grid_min:grid_max)) < &
318 : s% pg% abundance_log_mass_frac_min) return
319 : end if
320 0 : if (legend_flag) then
321 0 : do1 = abundance_line_legend(cnt, str)
322 : else
323 0 : do1 = abundance_line(cnt, i, str)
324 : end if
325 0 : end function do1
326 :
327 :
328 0 : subroutine set_line_style(cnt)
329 : integer, intent(in) :: cnt
330 : integer :: iclr, itype
331 0 : iclr = cnt - num_colors*(cnt/num_colors) + 1
332 0 : call pgsci(colors(iclr))
333 0 : if (cnt >= num_colors) then
334 0 : itype = Line_Type_Dot
335 : else
336 0 : itype = Line_Type_Solid
337 : end if
338 0 : call pgsls(itype)
339 0 : end subroutine set_line_style
340 :
341 :
342 0 : integer function abundance_line(cnt, j, str)
343 : integer, intent(in) :: cnt, j
344 : character (len=*), intent(in) :: str
345 : integer :: i, ii
346 0 : real :: x, frac, y, dx
347 : include 'formats'
348 0 : call set_line_style(cnt)
349 0 : call pgslw(lw)
350 0 : call pgline(npts, xvec(grid_min:grid_max), yvec(grid_min:grid_max))
351 0 : call pgslw(lw_sav)
352 0 : call pgsch(txt_scale*s% pg% Abundance_line_txt_scale_factor)
353 0 : abundance_line = cnt + 1
354 0 : do i=1,num_labels
355 0 : ii = iloc_abundance_label(i)
356 0 : if (ii > grid_max .or. ii < grid_min) cycle
357 0 : x = xloc_abundance_label(i)
358 0 : if (ii > grid_min) then
359 0 : dx = xvec(ii)-xvec(ii-1)
360 0 : if (abs(dx) > 1e-20) then
361 0 : frac = (x-xvec(ii-1))/dx
362 : else
363 : frac = 1.0
364 : end if
365 0 : y = frac*yvec(ii) + (1.0-frac)*yvec(ii-1) + dylbl
366 : else
367 0 : y = yvec(ii) + dylbl
368 : end if
369 0 : if (y < ymin .or. y > ymax) cycle
370 0 : call pgptxt(x, y, 0.0, 0.5, str)
371 : end do
372 0 : end function abundance_line
373 :
374 :
375 0 : integer function abundance_line_legend(cnt, str)
376 : integer, intent(in) :: cnt
377 : character (len=*), intent(in) :: str
378 0 : real :: dx, dyline, ypos, xpts(2), ypts(2)
379 : integer :: max_cnt
380 0 : max_cnt = min(max_num_labels, s% pg% Abundance_legend_max_cnt)
381 0 : if (cnt >= max_cnt) then
382 0 : abundance_line_legend = cnt
383 : return
384 : end if
385 0 : call set_line_style(cnt)
386 0 : dx = 0.1
387 0 : dyline = (ymax-ymin)/max_cnt
388 0 : ypos = ymax - (cnt+0.5)*dyline
389 0 : xpts(1) = 1.3*dx
390 0 : xpts(2) = xpts(1) + 2.3*dx
391 0 : ypts = ypos + dyline*0.1
392 0 : call pgslw(lw)
393 0 : call pgline(2, xpts, ypts)
394 0 : call pgslw(lw_sav)
395 0 : call pgsci(clr_Foreground)
396 0 : call pgsch(txt_scale*s% pg% Abundance_legend_txt_scale_factor)
397 0 : call pgptxt(xpts(2) + dx, ypos, 0.0, 0.0, trim(str))
398 0 : abundance_line_legend = cnt + 1
399 0 : end function abundance_line_legend
400 :
401 : end subroutine do_abundance_panel
402 :
403 : end module pgstar_abundance
|
