Plot_me.py
Go to the documentation of this file.
1 # Author: M. Reichert
2 import numpy as np
3 import matplotlib.pyplot as plt
4 from matplotlib import cm
5 # Import the WinNet python plot-routine
6 # Note that you can also just append the path to your .bashrc
7 import sys
8 sys.path.append('../../bin')
9 from class_files.winnet_class import winnet
10 
11 
12 # First, plot the final abundances
13 # Create figures to plot at
14 fig_finab = plt.figure(figsize=(5,3))
15 ax_finab = fig_finab.gca()
16 
17 # Plot the result with default reaclib decays
18 # The table has been created out of the reaclib data
19 path = "beta_decay_reaclib.dat/finabsum.dat"
20 A,Y,X = np.loadtxt(path,unpack=True)
21 ax_finab.plot(A,X,label='REACLIB')
22 
23 # Plot the Moeller beta decay data
24 # The decay data was accessed at
25 # https://www.matthewmumpower.com/publications/paper/2019/moller/nuclear-properties-for-astrophysical-and-radioactive-ion-beam-applications-ii
26 # and converted to the WinNet format
27 # The publication of Moeller et al. 2019 to these rates can be found here:
28 # https://www.sciencedirect.com/science/article/pii/S0092640X18300263?via%3Dihub
29 path = "beta_decay_moeller.dat/finabsum.dat"
30 A,Y,X = np.loadtxt(path,unpack=True)
31 ax_finab.plot(A,X,label='Möller et al. 2019')
32 
33 # Plot the result with Marketin beta decay data
34 # The data and paper of Marketin et al. 2016 can be
35 # accessed at:
36 # https://journals.aps.org/prc/abstract/10.1103/PhysRevC.93.025805
37 path = "beta_decay_marketin.dat/finabsum.dat"
38 A,Y,X = np.loadtxt(path,unpack=True)
39 ax_finab.plot(A,X,label='Marketin et al. 2016')
40 
41 # Make the plot look nice
42 ax_finab.legend(fontsize=7)
43 ax_finab.set_yscale('log')
44 ax_finab.set_ylim(1e-5,1)
45 ax_finab.set_xlim(55,240)
46 ax_finab.set_ylabel("Mass fraction")
47 ax_finab.set_xlabel("Mass number")
48 
49 # Also plot the mass fractions in the nuclear chart at freezeout
50 fig_chart,ax_chart = plt.subplots(1,3,figsize=(15,4),sharex=True,sharey=True)
51 plt.subplots_adjust(hspace=0,wspace=0)
52 
53 
54 # Setup for nuclear chart plot
55 lw = 0.1 # use a small linewidth for the nuclei-rectangle
56 min_X = 1e-12 # Min of mass fraction
57 max_X = 1e-4 # Max of mass fraction
58 cmap = cm.magma_r # Use the jet colormap
59 
60 # Reaclib
61 w = winnet('./beta_decay_reaclib.dat/')
62 w.read_mainout()
63 t,yn,yh = w.get_mainout_time(), w.get_mainout_yn(), w.get_mainout_yheavy()
64 # Find the neutron freezeout time (yn/yh = 1)
65 t_freezeout = t[np.argmin(np.abs(yn/yh-1))]
66 
67 # Plot the path in the nuclear chart
68 w.plot_nuclear_chart_at(t_freezeout,figure=ax_chart[0],axes_label=False,element_labels=False,fig_is_ax=True,
69  colorbar=False,min_X=min_X,max_X=max_X,cmap= cmap,nuclei_linewidths=lw)
70 # Label the plot
71 ax_chart[0].text(0.01,0.98,"Reaclib",transform=ax_chart[0].transAxes,va='top',ha='left',fontsize=13,color="grey")
72 
73 
74 # Moeller
75 w = winnet('./beta_decay_moeller.dat/')
76 w.read_mainout()
77 t,yn,yh = w.get_mainout_time(), w.get_mainout_yn(), w.get_mainout_yheavy()
78 # Find the neutron freezeout time (yn/yh = 1)
79 t_freezeout = t[np.argmin(np.abs(yn/yh-1))]
80 # Plot the path in the nuclear chart
81 w.plot_nuclear_chart_at(t_freezeout,figure=ax_chart[1],axes_label=False,element_labels=False,fig_is_ax=True,
82  colorbar=False,min_X=min_X,max_X=max_X,cmap= cmap,nuclei_linewidths=lw)
83 # Label the plot
84 ax_chart[1].text(0.01,0.98,"Möller et al. 2019",transform=ax_chart[1].transAxes,va='top',ha='left',fontsize=13,color="grey")
85 
86 
87 # Marketin
88 w = winnet('./beta_decay_marketin.dat/')
89 w.read_mainout()
90 t,yn,yh = w.get_mainout_time(), w.get_mainout_yn(), w.get_mainout_yheavy()
91 # Find the neutron freezeout time (yn/yh = 1)
92 t_freezeout = t[np.argmin(np.abs(yn/yh-1))]
93 # Plot the path in the nuclear chart
94 w.plot_nuclear_chart_at(t_freezeout,figure=ax_chart[2],axes_label=False,element_labels=False,fig_is_ax=True,colorbar=True,
95  colorbar_position=[0.27, 0.85, 0.5, 0.025],colorbar_inset=True,min_X=min_X,max_X=max_X
96  ,cmap= cmap,nuclei_linewidths=lw)
97 # Label the plot
98 ax_chart[2].text(0.01,0.98,"Marketin et al. 2016",transform=ax_chart[2].transAxes,va='top',ha='left',fontsize=13,color="grey")
99 
100 
101 # Set the labels
102 ax_chart[0].set_ylabel("Proton number")
103 ax_chart[1].set_xlabel("Neutron number")
104 # Set the correct limits
105 [a.set_xlim(50,150) for a in ax_chart]
106 [a.set_ylim(30,85) for a in ax_chart]
107 
108 # Save and show the figures
109 fig_finab.savefig("final_massfractions.pdf",bbox_inches='tight')
110 fig_chart.savefig("nuclear_chart_beta_decay.pdf",bbox_inches='tight')
111 plt.show()