Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2015-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_production
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 production_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_Production_plot(s, id, device_id, &
46 : s% pg% Production_xleft, s% pg% Production_xright, &
47 : s% pg% Production_ybot, s% pg% Production_ytop, .false., &
48 0 : s% pg% Production_title, s% pg% Production_txt_scale, ierr)
49 :
50 0 : call pgebuf()
51 :
52 : end subroutine production_plot
53 :
54 :
55 0 : subroutine do_production_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_production_panel(s, id, device_id, &
66 : winxmin, winxmax, winymin, winymax, subplot, &
67 0 : title, txt_scale, ierr)
68 0 : end subroutine do_production_plot
69 :
70 :
71 0 : subroutine do_production_panel(s, id, device_id, &
72 : winxmin, winxmax, winymin, winymax, subplot, title, txt_scale, &
73 : ierr)
74 : use utils_lib
75 : use chem_def
76 : use net_def
77 : use const_def, only: Msun
78 : use pgstar_colors
79 :
80 : type (star_info), pointer :: s
81 : integer, intent(in) :: id, device_id
82 : real, intent(in) :: winxmin, winxmax, winymin, winymax
83 : character (len=*), intent(in) :: title
84 : real, intent(in) :: txt_scale
85 : logical, intent(in) :: subplot
86 : integer, intent(out) :: ierr
87 :
88 0 : real :: xleft, xright, chScale, xmargin
89 : integer, parameter :: num_colors = 14
90 : integer :: colors(num_colors)
91 :
92 : include 'formats'
93 : ierr = 0
94 :
95 : colors(:) = [ &
96 : clr_Gold, clr_LightSkyBlue, clr_Crimson, clr_Goldenrod, clr_MediumSlateBlue, &
97 : clr_Coral, clr_LightSkyGreen, clr_DarkGray, clr_Lilac, &
98 0 : clr_Tan, clr_IndianRed, clr_Teal, clr_Silver, clr_BrightBlue ]
99 :
100 0 : chScale = txt_scale
101 :
102 0 : xmargin = 0
103 0 : call plot(ierr)
104 :
105 :
106 : contains
107 :
108 0 : subroutine plot(ierr)
109 0 : use chem_def
110 : use chem_lib
111 : use adjust_xyz, only: get_xa_for_standard_metals
112 : integer, intent(out) :: ierr
113 :
114 : integer :: i, j
115 :
116 : integer :: amin,amax,z,n,a,zmin,zmax
117 : integer :: min_zone,max_zone,alternate
118 0 : real :: extra_pad
119 0 : real :: min_mass,max_mass,yloc
120 : real,parameter :: point_size=0.1
121 0 : real :: ymin, ymax
122 : real,parameter :: pad=1.0
123 0 : real :: last_x,last_y
124 0 : real(dp),dimension(1:solsiz) :: scaled_abun,scaled_abun_init
125 0 : real(dp),dimension(:),allocatable :: init_comp,abun
126 :
127 0 : real(dp) :: initial_z,initial_y,initial_h1,initial_h2
128 0 : real(dp) :: initial_he3,initial_he4,xsol_he3,xsol_he4,la,lac
129 :
130 : include 'formats'
131 0 : ierr = 0
132 :
133 0 : call pgsave
134 0 : call pgsch(txt_scale)
135 0 : call pgsvp(winxmin, winxmax, winymin, winymax)
136 :
137 0 : amax=0
138 0 : amin=0.0
139 0 : zmax=0
140 0 : zmin=HUGE(zmin)
141 0 : ymax=-HUGE(ymax)
142 0 : ymin=HUGE(ymin)
143 :
144 :
145 0 : if (s% pg% production_min_mass > 0.0) then
146 0 : min_mass = s% pg% production_min_mass
147 : else
148 : min_mass = 0.d0
149 : end if
150 :
151 0 : if (s% pg% production_max_mass > -100) then
152 : max_mass = s% pg% production_max_mass
153 : else
154 0 : max_mass = s% mstar/msun
155 : end if
156 :
157 0 : min_zone=1
158 0 : max_zone=s%nz
159 :
160 0 : do i=1,s%nz
161 0 : if(s%m(i)<=max_mass) then
162 0 : min_zone=i-1
163 0 : exit
164 : end if
165 : end do
166 :
167 0 : do i=min_zone,s%nz
168 0 : if(s%m(i)<=min_mass) then
169 0 : max_zone=i-1
170 0 : exit
171 : end if
172 : end do
173 :
174 0 : allocate(abun(1:s%species),init_comp(1:s%species))
175 0 : abun=0.d0
176 :
177 : !Get initial composition
178 :
179 : !Stolen from star/private/create_initial_model
180 0 : initial_z = s% initial_z
181 0 : initial_y = s% initial_y
182 0 : initial_h1 = max(0d0, min(1d0, 1d0 - (initial_z + initial_y)))
183 0 : initial_h2 = chem_Xsol('h2')
184 0 : xsol_he3 = chem_Xsol('he3')
185 0 : xsol_he4 = chem_Xsol('he4')
186 0 : initial_he3 = initial_y*xsol_he3/(xsol_he3 + xsol_he4)
187 0 : initial_he4 = initial_y*xsol_he4/(xsol_he3 + xsol_he4)
188 : !
189 :
190 : call get_xa_for_standard_metals( &
191 : s, s%species, s%chem_id, s%net_iso,&
192 : initial_h1,initial_h2,initial_he3,initial_he4, &
193 0 : s% job% initial_zfracs, s% job% dump_missing_metals_into_heaviest, init_comp, ierr)
194 :
195 : !compute abundances
196 0 : do i=1,s%species
197 :
198 0 : Z=chem_isos%Z(s%chem_id(i))
199 0 : N=chem_isos%N(s%chem_id(i))
200 :
201 0 : zmin=min(Z,zmin)
202 0 : zmax=max(Z,zmax)
203 :
204 : abun(i)=(dot_product(s%xa(i,min_zone:max_zone),s%dm(min_zone:max_zone))/msun)/&
205 0 : ((max_mass-min_mass)-(s%m_center/msun))
206 :
207 : end do
208 :
209 : !Get stable isotope abundances
210 0 : call get_stable_mass_frac(s%chem_id,s%species,dble(abun),scaled_abun)
211 0 : call get_stable_mass_frac(s%chem_id,s%species,init_comp,scaled_abun_init)
212 :
213 0 : do i=1,solsiz
214 :
215 0 : la=safe_log10(scaled_abun(i))
216 0 : lac=safe_log10(scaled_abun_init(i))
217 :
218 : !Remove low abundance isotopes, low in star and low in solar can lead to large production factor
219 0 : if(la <s% pg%production_min_mass_frac .or. lac <s% pg%production_min_mass_frac) then
220 0 : scaled_abun(i)=-HUGE(ymin)
221 : else
222 0 : scaled_abun(i)=real(la-lac)
223 0 : ymax=max(ymax,real(scaled_abun(i),kind=kind(ymax)))
224 0 : ymin=min(ymin,real(scaled_abun(i),kind=kind(ymin)))
225 : end if
226 :
227 0 : if(zmax==izsol(i)) amax=int(iasol(i))
228 :
229 : end do
230 :
231 0 : if (s% pg% production_ymax > -100) then
232 : ymax = s% pg% production_ymax
233 : else
234 0 : ymax = ymax+0.01
235 : end if
236 :
237 0 : if (s% pg% production_ymin > -100) then
238 0 : ymin = s% pg% production_ymin
239 : else
240 : ymin = ymin
241 : end if
242 :
243 0 : if (s% pg% production_amax > -100) then
244 0 : xright = s% pg% production_amax
245 : else
246 0 : xright= amax
247 : end if
248 :
249 0 : if (s% pg% production_amin > -100) then
250 0 : xleft = s% pg% production_amin
251 : else
252 0 : xleft = amin
253 : end if
254 :
255 0 : if (s% pg% Production_show_element_names) THEN
256 : extra_pad=2.5
257 : else
258 0 : extra_pad=1.0
259 : end if
260 :
261 : !Set xaxis and yaxis bounds
262 0 : call pgswin(xleft-1,xright+1,ymin-pad,ymax+extra_pad*pad)
263 : !Create a box with ticks
264 0 : call show_box_pgstar(s,'BCNST','BCNSTV')
265 : !Labels
266 0 : call show_xaxis_name(s,'A',ierr)
267 0 : call show_left_yaxis_label_pgstar(s,'Log Production Factor',0.0)
268 :
269 0 : if (.not. subplot) then
270 0 : call show_model_number_pgstar(s)
271 0 : call show_age_pgstar(s)
272 : end if
273 0 : call show_title_pgstar(s, title)
274 :
275 0 : alternate=-1
276 0 : i=1
277 : outer: do
278 0 : if(i>solsiz) exit outer
279 :
280 : ! Z is greater than zmax
281 0 : if(izsol(i)>zmax) exit outer
282 :
283 : !Sets color
284 0 : call set_line_style(izsol(i))
285 :
286 : !Shows element name, alternates between two levels to spread them out
287 : if (s% pg% Production_show_element_names &
288 0 : .and.(iasol(i)<=xright).and.(iasol(i)>=xleft)) THEN
289 0 : yloc=(ymax*1.0)+abs(0.75+(alternate)/2.0)
290 0 : call pgtext(iasol(i)*1.0,yloc,el_name(izsol(i)))
291 0 : alternate=alternate*(-1.0)
292 : end if
293 :
294 : !When we have more than one isotope per element we draw a line, this marks the last
295 : ! x cordinate we "saw" for an element
296 0 : last_x=-HUGE(last_x)
297 0 : last_y=-HUGE(last_y)
298 :
299 0 : inner: do j=i,solsiz
300 :
301 0 : if(izsol(j)==izsol(i))then
302 : if((scaled_abun(j)>= ymin) .and. (scaled_abun(j) <= ymax)&
303 0 : .and.(iasol(j)<=xright).and.(iasol(j)>=xleft)) then
304 0 : a=iasol(j)
305 :
306 : !Draw point at values
307 0 : call PGCIRC(A*1.0,real(scaled_abun(j)),point_size)
308 : !Not the first isotope of an element
309 0 : if(last_x>0)then
310 : !Then draw a line between isotopes of same element
311 0 : call pgline(2,[last_x,A*1.0],[last_y,real(scaled_abun(j)*1.0)])
312 : end if
313 :
314 : !Save last x,y pair we saw
315 0 : last_x=A*1.0
316 0 : last_y=scaled_abun(j)
317 : end if
318 : else
319 : exit inner
320 : end if
321 :
322 : end do inner
323 : !Jump to next element
324 0 : i=j
325 : end do outer
326 :
327 0 : call pgunsa
328 0 : deallocate(abun,init_comp)
329 :
330 : call show_pgstar_decorator(id,s% pg% production_use_decorator,&
331 0 : s% pg% production_pgstar_decorator, 0, ierr)
332 :
333 :
334 0 : end subroutine plot
335 :
336 0 : subroutine set_line_style(cnt)
337 : integer, intent(in) :: cnt
338 : integer :: iclr, itype
339 0 : iclr = cnt - num_colors*(cnt/num_colors) + 1
340 0 : call pgsci(colors(iclr))
341 0 : if (cnt >= num_colors) then
342 0 : itype = Line_Type_Dot
343 : else
344 0 : itype = Line_Type_Solid
345 : end if
346 0 : call pgsls(itype)
347 0 : end subroutine set_line_style
348 :
349 : end subroutine do_production_panel
350 :
351 : end module pgstar_production
|
