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 import h5py
5 import os
6 
7 # This script will plot the different contributions to the energy generation
8 # from the (cold) CNO-cycles, pp-chain, and Ne-Na, as well as Mg-Al cycles.
9 # Note that given temperatures should not be smaller than 1d-2GK
10 # which is the lower limit of the reaclib fits.
11 # The energy generation of the reactions is not written to the output individually
12 # and, as a consequence, the plot is only an approximation to some degree.
13 
14 # Get names of the runs, ignore hidden files
15 paths = [i for i in os.listdir(".") if i[0] != "."]
16 
17 # Empty list to store various quantites
18 
19 # Loop through all folders
20 for p in paths:
21  # Skip non folders
22  if not os.path.isdir(p):
23  continue
24  if p.strip() == "network_data":
25  continue
26  if not os.path.isfile(os.path.join(p,"WinNet_data.h5")):
27  raise Exception("Could not find WinNet_data.h5 in '"+str(p)+"'. "\
28  "Are you sure you enabled hdf5 output in the makefile "\
29  "before running?")
30 
31  # Open the file
32  f = h5py.File(os.path.join(p,"WinNet_data.h5"),"r")
33 
34  # A bit handwaving formula to get the energy plateau at the end.
35  # This plateau can be estimated in a more physical motivated way rather
36  # then fitting it by hand as done here. To estimate it, one can
37  # consider when the pariticipating nuclei are in equilibrium and take
38  # this time. A detailed derivation is given in
39  # e.g., Nuclear physics of stars, C. Iliadis, Sect. 5.1.
40  time_plateau = 10**(-4/2*float(p)+17.5)
41 
42  # Total energy
43  time = f["energy/time"][:]
44  idx_plateau = np.argmin(abs(time-time_plateau))
45 
46  # Get nuclear properties
47  A = f["energy/A"][:]
48  Z = f["energy/Z"][:]
49  N = f["energy/N"][:]
50 
51  # In the following we separate the energy generation from
52  # the different cycles. The reactions and cycles can e.g.,
53  # be seen at https://cococubed.com/code_pages/burn_hydrogen.shtml .
54 
55  # CNO-cycles and pp-chain
56  idx_p = (A == 1) & (Z ==1)
57  idx_d = (A == 2) & (Z ==1)
58  idx_he3 = (A == 3) & (Z ==2)
59  idx_he4 = (A == 4) & (Z ==2)
60  idx_be7 = (A == 7) & (Z ==4)
61  idx_be8 = (A == 8) & (Z ==4)
62  idx_b8 = (A == 8) & (Z ==5)
63  idx_li7 = (A == 7) & (Z ==3)
64  idx_c12 = (A == 12) & (Z ==6)
65  idx_c13 = (A == 13) & (Z ==6)
66  idx_n13 = (A == 13) & (Z ==7)
67  idx_n14 = (A == 14) & (Z ==7)
68  idx_n15 = (A == 15) & (Z ==7)
69  idx_o15 = (A == 15) & (Z ==8)
70  idx_o16 = (A == 16) & (Z ==8)
71  idx_o17 = (A == 17) & (Z ==8)
72  idx_o18 = (A == 18) & (Z ==8)
73  idx_f17 = (A == 17) & (Z ==9)
74  idx_f18 = (A == 18) & (Z ==9)
75  idx_f19 = (A == 19) & (Z ==9)
76 
77  e_detailed_pg = f["energy/detailed (p,g)"][idx_plateau,:]
78  e_detailed_bet = f["energy/detailed decay"][idx_plateau,:]
79  e_detailed_ap = f["energy/detailed (a,p)"][idx_plateau,:]
80  e_detailed_other = f["energy/detailed other"][idx_plateau,:]
81  e_detailed_ag = f["energy/detailed (a,g)"][idx_plateau,:]
82 
83  #############
84  # CNO cycle #
85  #############
86 
87  # CNO cycle 1
88  # (p,g) reactions
89  c13pg = e_detailed_pg[idx_c13]
90  n14pg = e_detailed_pg[idx_n14]
91  c12pg = e_detailed_pg[idx_c12]
92  # Decays
93  n13d = e_detailed_bet[idx_n13]
94  o15d = e_detailed_bet[idx_o15]
95  # (p,a) reactions
96  n15pa = -e_detailed_ap[idx_n15]
97  # Total CNO 1
98  e_cno = c13pg+n14pg+c12pg+n13d+o15d+n15pa
99 
100  # CNO cycle 2
101  # (p,g) reactions
102  n15pg = e_detailed_pg[idx_n15]
103  o16pg = e_detailed_pg[idx_o16]
104  # Decays
105  f17d = e_detailed_bet[idx_f17]
106  # (p,a) reactions
107  o17pa = -e_detailed_ap[idx_o17]
108  # Total CNO 1 + CNO 2
109  e_cno = e_cno+n15pg+o16pg+f17d+o17pa
110 
111  # CNO cycle 3
112  # (p,g) reactions
113  o17pg = e_detailed_pg[idx_o17]
114  # Decays
115  f18pg = e_detailed_bet[idx_f17]
116  # (p,a) reactions
117  o18pa = -e_detailed_ap[idx_o18]
118  # Total CNO 1 + CNO 2 + CNO 3
119  e_cno = e_cno+o17pg+f18pg+o18pa
120 
121  # CNO cycle 4
122  # (p,g) reactions
123  o18pg = e_detailed_pg[idx_o18]
124  # (p,a) reactions
125  f19pa = -e_detailed_ap[idx_f19]
126  # Total CNO 1 + CNO 2 + CNO 3
127  e_cno = e_cno+o18pg+f19pa
128 
129  #############
130  # PP-Chains #
131  #############
132  ppg = e_detailed_pg[idx_p]
133  dpg = e_detailed_pg[idx_d]
134  pd = e_detailed_bet[idx_p]
135  # Branch 1
136  he3he3 = e_detailed_other[idx_he3]
137  # Branch 2 and 3
138  he3a = e_detailed_ag[idx_he3]
139  be7d = e_detailed_bet[idx_be7]
140  li7p = e_detailed_other[idx_li7]
141  be7pg = e_detailed_pg[idx_be7]
142  b8d = e_detailed_bet[idx_b8]
143  be8o = e_detailed_other[idx_be8]
144 
145  # Total energy generation of pp-chain
146  e_pp = ppg+dpg+he3he3+pd+he3a+li7p+be7d+b8d+be8o
147 
148  #########################
149  # Ne-Na and Mg-Al chain #
150  #########################
151  idx_ne20 = (A == 20) & (Z ==10)
152  idx_ne21 = (A == 21) & (Z ==10)
153  idx_ne22 = (A == 22) & (Z ==10)
154  idx_na21 = (A == 21) & (Z ==11)
155  idx_na22 = (A == 22) & (Z ==11)
156  idx_na23 = (A == 23) & (Z ==11)
157  idx_mg22 = (A == 22) & (Z ==12)
158  idx_mg23 = (A == 23) & (Z ==12)
159  idx_mg24 = (A == 24) & (Z ==12)
160  idx_mg25 = (A == 25) & (Z ==12)
161  idx_mg26 = (A == 26) & (Z ==12)
162  idx_al25 = (A == 25) & (Z ==13)
163  idx_al26 = (A == 26) & (Z ==13)
164  idx_al27 = (A == 27) & (Z ==13)
165  idx_si26 = (A == 26) & (Z ==14)
166  idx_si27 = (A == 27) & (Z ==14)
167 
168  # Ne-Na chain
169  #(p,g)
170  ne20pg = e_detailed_pg[idx_ne20]
171  ne21pg = e_detailed_pg[idx_ne21]
172  na21pg = e_detailed_pg[idx_na21]
173  na22pg = e_detailed_pg[idx_na22]
174  ne22pg = e_detailed_pg[idx_ne22]
175  #decay
176  na21d = e_detailed_bet[idx_na21]
177  mg22d = e_detailed_bet[idx_mg22]
178  na22d = e_detailed_bet[idx_na22]
179  mg23d = e_detailed_bet[idx_mg23]
180  #(p,a)
181  na23pa = -e_detailed_ap[idx_na23]
182 
183  e_nena = ne20pg+ne21pg+na21pg+na22pg+ne22pg+na21d+mg22d+na22d+mg23d+na23pa
184 
185  # Mg-al chain
186  #(p,g)
187  na23pg = e_detailed_pg[idx_na23]
188  mg24pg = e_detailed_pg[idx_mg24]
189  al25pg = e_detailed_pg[idx_al25]
190  mg25pg = e_detailed_pg[idx_mg25]
191  al26pg = np.sum(e_detailed_pg[idx_al26])
192  mg26pg = e_detailed_pg[idx_mg26]
193  #decay
194  al25d = e_detailed_bet[idx_al25]
195  si26d = e_detailed_bet[idx_si26]
196  al26d = np.sum(e_detailed_bet[idx_al26])
197  si27d = e_detailed_bet[idx_si27]
198  #(p,a)
199  al27pa = -e_detailed_ap[idx_al27]
200 
201  e_mgal = na23pg+mg24pg+al25pg+mg25pg+al26pg+mg26pg+al25d+si26d+al26d+si27d+al27pa
202 
203  # Plot Ne-Na cycle
204  plt.scatter(float(p),e_nena,color="tab:orange")
205  # Plot Mg-Al cycle
206  plt.scatter(float(p),e_mgal,color="tab:green")
207  # Plot PP-Chains
208  plt.scatter(float(p),e_pp,color="tab:blue")
209  # Plot CNO cycles
210  plt.scatter(float(p),e_cno,color="tab:red")
211 
212 
213 # Indicate the suns core temperature
214 plt.axvline(1.56,lw=2,ls="--",color="lightgrey")
215 plt.text(1.56,5e6,r"Temperature of the sun",rotation=90,va="top",ha="right",fontsize=14,color="lightgrey")
216 
217 # Create labels for the legend
218 plt.scatter(np.nan,np.nan,color="tab:red" ,label="CNO-cycle")
219 plt.scatter(np.nan,np.nan,color="tab:blue" ,label="pp-chain")
220 plt.scatter(np.nan,np.nan,color="tab:orange",label="Ne-Na cycle")
221 plt.scatter(np.nan,np.nan,color="tab:green",label="Mg-Al cycle")
222 plt.legend()
223 
224 # Set the limits of the plot
225 plt.ylim(3e0,1e7)
226 plt.xlim(1,4)
227 plt.yscale("log")
228 plt.xlabel(r"Temperature [10$^7$ K]")
229 plt.ylabel(r"Energy [erg g$^{-1}$ s$^{-1}$]")
230 plt.title(r"$\rho$ = 100 g cm$^{-3}$")
231 plt.savefig("energy_generation.pdf",bbox_inches="tight")
232 plt.show()