A He-burning Network#
Here we create a network that can be used as an \(\alpha\)-chain for modeling He detonations, and includes some of the iron group. It is similar in goals to the classic aprox19 network, but it has some additional rates that are important.
import pynucastro as pyna
Assembling our library#
Core nuclei#
We start with a list of nuclei, including all \(\alpha\)-nuclei up to \({}^{56}\mathrm{Ni}\).
We also add the intermediate nuclei that would participate in \((\alpha,p)(p,\gamma)\) reactions,
as well as a few more nuclei, including those identified by Shen & Bildsten 2009 for bypassing the \({}^{12}\mathrm{C}(\alpha,\gamma){}^{16}\mathrm{O}\) rate. These rates are not present in the classic aprox networks, but can
be very important for the ignition of a detonation.
nuclei = ["p",
"he4", "c12", "o16", "ne20", "mg24", "si28", "s32",
"ar36", "ca40", "ti44", "cr48", "fe52", "ni56",
"al27", "p31", "cl35", "k39", "sc43", "v47", "mn51", "co55",
"n13", "n14", "f18", "ne21", "na22", "na23"]
For our initial library, we take all of the ReacLib rates that link these.
reaclib_lib = pyna.ReacLibLibrary()
core_lib = reaclib_lib.linking_nuclei(nuclei)
Since we didn’t include neutrons in our list of nuclei, we are missing some potentially important rates, which we now add manually. However, we do not want to carry neutrons, so we modify the endpoints of these reactions to assume that the original \((X, n)\) reaction is immediately followed by a neutron capture.
other_rates = [("c12(c12,n)mg23", "mg24"),
("o16(o16,n)s31", "s32"),
("o16(c12,n)si27", "si28")]
for r, mp in other_rates:
_r = reaclib_lib.get_rate_by_name(r)
_r.modify_products(mp)
core_lib += pyna.Library(rates=[_r])
Now we remove some rates that traditionally are not present in the aprox-style networks. This includes the
reverse rates of \({}^{12}\mathrm{C} + {}^{12}\mathrm{C}\), \({}^{16}\mathrm{O} + {}^{12}\mathrm{C}\), and \({}^{16}\mathrm{O} + {}^{16}\mathrm{O}\). We also remove some \(\alpha\)-captures and \(p\)-captures on neutron-rich nuclei that will be eliminated when we approximate them out of the network.
for r in core_lib.get_rates():
if sorted(r.products) in [[pyna.Nucleus("c12"), pyna.Nucleus("c12")],
[pyna.Nucleus("c12"), pyna.Nucleus("o16")],
[pyna.Nucleus("o16"), pyna.Nucleus("o16")]]:
core_lib.remove_rate(r)
rates_to_remove = ["p31(p,c12)ne20",
"si28(a,c12)ne20",
"ne20(c12,p)p31",
"ne20(c12,a)si28",
"na23(a,g)al27",
"al27(g,a)na23",
"al27(a,g)p31",
"p31(g,a)al27"]
for r in rates_to_remove:
_r = core_lib.get_rate_by_name(r)
core_lib.remove_rate(_r)
Here’s the current list of rates
core_lib
3 He4 ⟶ C12 + 𝛾 [Q = 7.28 MeV] (3 He4 --> C12 <fy05_reaclib__>)
C12 + p ⟶ N13 + 𝛾 [Q = 1.94 MeV] (C12 + p --> N13 <ls09_reaclib__>)
C12 + He4 ⟶ O16 + 𝛾 [Q = 7.16 MeV] (C12 + He4 --> O16 <nac2_reaclib__>)
C12 + C12 ⟶ He4 + Ne20 [Q = 4.62 MeV] (C12 + C12 --> He4 + Ne20 <cf88_reaclib__>)
C12 + C12 ⟶ p + Na23 [Q = 2.24 MeV] (C12 + C12 --> p + Na23 <cf88_reaclib__>)
C12 + C12 ⟶ Mg24 + 𝛾 [Q = 13.93 MeV] (C12 + C12 --> Mg24 <cf88_reaclib__derived_from_inverse>)
N13 + He4 ⟶ p + O16 [Q = 5.22 MeV] (N13 + He4 --> p + O16 <cf88_reaclib__>)
N14 + He4 ⟶ F18 + 𝛾 [Q = 4.41 MeV] (N14 + He4 --> F18 <il10_reaclib__>)
O16 + He4 ⟶ Ne20 + 𝛾 [Q = 4.73 MeV] (O16 + He4 --> Ne20 <co10_reaclib__>)
O16 + C12 ⟶ He4 + Mg24 [Q = 6.77 MeV] (O16 + C12 --> He4 + Mg24 <cf88_reaclib__>)
O16 + C12 ⟶ p + Al27 [Q = 5.17 MeV] (O16 + C12 --> p + Al27 <cf88_reaclib__>)
O16 + C12 ⟶ Si28 + 𝛾 [Q = 16.76 MeV] (O16 + C12 --> Si28 <cf88_reaclib__derived_from_inverse>)
O16 + O16 ⟶ He4 + Si28 [Q = 9.59 MeV] (O16 + O16 --> He4 + Si28 <cf88_reaclib__>)
O16 + O16 ⟶ p + P31 [Q = 7.68 MeV] (O16 + O16 --> p + P31 <cf88_reaclib__>)
O16 + O16 ⟶ S32 + 𝛾 [Q = 16.54 MeV] (O16 + O16 --> S32 <cf88_reaclib__>)
F18 + He4 ⟶ p + Ne21 [Q = 1.74 MeV] (F18 + He4 --> p + Ne21 <rpsm_reaclib__>)
F18 + He4 ⟶ Na22 + 𝛾 [Q = 8.48 MeV] (F18 + He4 --> Na22 <rpsm_reaclib__>)
Ne20 + He4 ⟶ Mg24 + 𝛾 [Q = 9.32 MeV] (Ne20 + He4 --> Mg24 <il10_reaclib__>)
Ne21 + p ⟶ Na22 + 𝛾 [Q = 6.74 MeV] (Ne21 + p --> Na22 <il10_reaclib__>)
Na23 + p ⟶ He4 + Ne20 [Q = 2.38 MeV] (Na23 + p --> He4 + Ne20 <il10_reaclib__>)
Na23 + p ⟶ Mg24 + 𝛾 [Q = 11.69 MeV] (Na23 + p --> Mg24 <il10_reaclib__>)
Mg24 + He4 ⟶ Si28 + 𝛾 [Q = 9.98 MeV] (Mg24 + He4 --> Si28 <st08_reaclib__>)
Al27 + p ⟶ He4 + Mg24 [Q = 1.60 MeV] (Al27 + p --> He4 + Mg24 <il10_reaclib__>)
Al27 + p ⟶ Si28 + 𝛾 [Q = 11.59 MeV] (Al27 + p --> Si28 <il10_reaclib__>)
Si28 + He4 ⟶ S32 + 𝛾 [Q = 6.95 MeV] (Si28 + He4 --> S32 <ths8_reaclib__>)
P31 + p ⟶ He4 + Si28 [Q = 1.92 MeV] (P31 + p --> He4 + Si28 <il10_reaclib__>)
P31 + p ⟶ S32 + 𝛾 [Q = 8.86 MeV] (P31 + p --> S32 <il10_reaclib__>)
P31 + He4 ⟶ Cl35 + 𝛾 [Q = 7.00 MeV] (P31 + He4 --> Cl35 <ths8_reaclib__>)
S32 + He4 ⟶ Ar36 + 𝛾 [Q = 6.64 MeV] (S32 + He4 --> Ar36 <ths8_reaclib__>)
Cl35 + p ⟶ He4 + S32 [Q = 1.87 MeV] (Cl35 + p --> He4 + S32 <il10_reaclib__>)
Cl35 + p ⟶ Ar36 + 𝛾 [Q = 8.51 MeV] (Cl35 + p --> Ar36 <il10_reaclib__>)
Cl35 + He4 ⟶ K39 + 𝛾 [Q = 7.22 MeV] (Cl35 + He4 --> K39 <ths8_reaclib__>)
Ar36 + He4 ⟶ Ca40 + 𝛾 [Q = 7.04 MeV] (Ar36 + He4 --> Ca40 <ths8_reaclib__>)
K39 + p ⟶ He4 + Ar36 [Q = 1.29 MeV] (K39 + p --> He4 + Ar36 <ths8_reaclib__>)
K39 + p ⟶ Ca40 + 𝛾 [Q = 8.33 MeV] (K39 + p --> Ca40 <lo18_reaclib__>)
K39 + He4 ⟶ Sc43 + 𝛾 [Q = 4.81 MeV] (K39 + He4 --> Sc43 <ths8_reaclib__>)
Ca40 + He4 ⟶ Ti44 + 𝛾 [Q = 5.13 MeV] (Ca40 + He4 --> Ti44 <chw0_reaclib__>)
Sc43 + p ⟶ He4 + Ca40 [Q = 3.52 MeV] (Sc43 + p --> He4 + Ca40 <ths8_reaclib__>)
Sc43 + p ⟶ Ti44 + 𝛾 [Q = 8.65 MeV] (Sc43 + p --> Ti44 <ths8_reaclib__>)
Sc43 + He4 ⟶ V47 + 𝛾 [Q = 8.24 MeV] (Sc43 + He4 --> V47 <ths8_reaclib__>)
Ti44 + He4 ⟶ Cr48 + 𝛾 [Q = 7.70 MeV] (Ti44 + He4 --> Cr48 <ths8_reaclib__>)
V47 + p ⟶ He4 + Ti44 [Q = 0.41 MeV] (V47 + p --> He4 + Ti44 <chw0_reaclib__derived_from_inverse>)
V47 + p ⟶ Cr48 + 𝛾 [Q = 8.11 MeV] (V47 + p --> Cr48 <nfis_reaclib__>)
V47 + He4 ⟶ Mn51 + 𝛾 [Q = 8.66 MeV] (V47 + He4 --> Mn51 <ths8_reaclib__>)
Cr48 + He4 ⟶ p + Mn51 [Q = 0.56 MeV] (Cr48 + He4 --> p + Mn51 <ths8_reaclib__>)
Cr48 + He4 ⟶ Fe52 + 𝛾 [Q = 7.94 MeV] (Cr48 + He4 --> Fe52 <ths8_reaclib__>)
Mn51 + p ⟶ Fe52 + 𝛾 [Q = 7.38 MeV] (Mn51 + p --> Fe52 <ths8_reaclib__>)
Mn51 + He4 ⟶ Co55 + 𝛾 [Q = 8.21 MeV] (Mn51 + He4 --> Co55 <ths8_reaclib__>)
Fe52 + He4 ⟶ p + Co55 [Q = 0.83 MeV] (Fe52 + He4 --> p + Co55 <ths8_reaclib__>)
Fe52 + He4 ⟶ Ni56 + 𝛾 [Q = 8.00 MeV] (Fe52 + He4 --> Ni56 <ths8_reaclib__>)
Co55 + p ⟶ Ni56 + 𝛾 [Q = 7.17 MeV] (Co55 + p --> Ni56 <ths8_reaclib__>)
C12 ⟶ 3 He4 [Q = -7.28 MeV] (C12 --> 3 He4 <fy05_reaclib__derived_from_inverse>)
N13 ⟶ p + C12 [Q = -1.94 MeV] (N13 --> p + C12 <ls09_reaclib__derived_from_inverse>)
O16 + p ⟶ He4 + N13 [Q = -5.22 MeV] (O16 + p --> He4 + N13 <cf88_reaclib__derived_from_inverse>)
O16 ⟶ He4 + C12 [Q = -7.16 MeV] (O16 --> He4 + C12 <nac2_reaclib__derived_from_inverse>)
F18 ⟶ He4 + N14 [Q = -4.41 MeV] (F18 --> He4 + N14 <il10_reaclib__derived_from_inverse>)
Ne20 + He4 ⟶ p + Na23 [Q = -2.38 MeV] (Ne20 + He4 --> p + Na23 <il10_reaclib__derived_from_inverse>)
Ne20 ⟶ He4 + O16 [Q = -4.73 MeV] (Ne20 --> He4 + O16 <co10_reaclib__derived_from_inverse>)
Ne21 + p ⟶ He4 + F18 [Q = -1.74 MeV] (Ne21 + p --> He4 + F18 <rpsm_reaclib__derived_from_inverse>)
Na22 ⟶ He4 + F18 [Q = -8.48 MeV] (Na22 --> He4 + F18 <rpsm_reaclib__derived_from_inverse>)
Na22 ⟶ p + Ne21 [Q = -6.74 MeV] (Na22 --> p + Ne21 <il10_reaclib__derived_from_inverse>)
Mg24 + He4 ⟶ p + Al27 [Q = -1.60 MeV] (Mg24 + He4 --> p + Al27 <il10_reaclib__derived_from_inverse>)
Mg24 ⟶ He4 + Ne20 [Q = -9.32 MeV] (Mg24 --> He4 + Ne20 <il10_reaclib__derived_from_inverse>)
Mg24 ⟶ p + Na23 [Q = -11.69 MeV] (Mg24 --> p + Na23 <il10_reaclib__derived_from_inverse>)
Si28 + He4 ⟶ p + P31 [Q = -1.92 MeV] (Si28 + He4 --> p + P31 <il10_reaclib__derived_from_inverse>)
Si28 ⟶ He4 + Mg24 [Q = -9.98 MeV] (Si28 --> He4 + Mg24 <st08_reaclib__derived_from_inverse>)
Si28 ⟶ p + Al27 [Q = -11.59 MeV] (Si28 --> p + Al27 <il10_reaclib__derived_from_inverse>)
S32 + He4 ⟶ p + Cl35 [Q = -1.87 MeV] (S32 + He4 --> p + Cl35 <il10_reaclib__derived_from_inverse>)
S32 ⟶ He4 + Si28 [Q = -6.95 MeV] (S32 --> He4 + Si28 <ths8_reaclib__derived_from_inverse>)
S32 ⟶ p + P31 [Q = -8.86 MeV] (S32 --> p + P31 <il10_reaclib__derived_from_inverse>)
Cl35 ⟶ He4 + P31 [Q = -7.00 MeV] (Cl35 --> He4 + P31 <ths8_reaclib__derived_from_inverse>)
Ar36 + He4 ⟶ p + K39 [Q = -1.29 MeV] (Ar36 + He4 --> p + K39 <ths8_reaclib__derived_from_inverse>)
Ar36 ⟶ He4 + S32 [Q = -6.64 MeV] (Ar36 --> He4 + S32 <ths8_reaclib__derived_from_inverse>)
Ar36 ⟶ p + Cl35 [Q = -8.51 MeV] (Ar36 --> p + Cl35 <il10_reaclib__derived_from_inverse>)
K39 ⟶ He4 + Cl35 [Q = -7.22 MeV] (K39 --> He4 + Cl35 <ths8_reaclib__derived_from_inverse>)
Ca40 + He4 ⟶ p + Sc43 [Q = -3.52 MeV] (Ca40 + He4 --> p + Sc43 <ths8_reaclib__derived_from_inverse>)
Ca40 ⟶ He4 + Ar36 [Q = -7.04 MeV] (Ca40 --> He4 + Ar36 <ths8_reaclib__derived_from_inverse>)
Ca40 ⟶ p + K39 [Q = -8.33 MeV] (Ca40 --> p + K39 <lo18_reaclib__derived_from_inverse>)
Sc43 ⟶ He4 + K39 [Q = -4.81 MeV] (Sc43 --> He4 + K39 <ths8_reaclib__derived_from_inverse>)
Ti44 + He4 ⟶ p + V47 [Q = -0.41 MeV] (Ti44 + He4 --> p + V47 <chw0_reaclib__>)
Ti44 ⟶ He4 + Ca40 [Q = -5.13 MeV] (Ti44 --> He4 + Ca40 <chw0_reaclib__derived_from_inverse>)
Ti44 ⟶ p + Sc43 [Q = -8.65 MeV] (Ti44 --> p + Sc43 <ths8_reaclib__derived_from_inverse>)
V47 ⟶ He4 + Sc43 [Q = -8.24 MeV] (V47 --> He4 + Sc43 <ths8_reaclib__derived_from_inverse>)
Cr48 ⟶ He4 + Ti44 [Q = -7.70 MeV] (Cr48 --> He4 + Ti44 <ths8_reaclib__derived_from_inverse>)
Cr48 ⟶ p + V47 [Q = -8.11 MeV] (Cr48 --> p + V47 <nfis_reaclib__derived_from_inverse>)
Mn51 + p ⟶ He4 + Cr48 [Q = -0.56 MeV] (Mn51 + p --> He4 + Cr48 <ths8_reaclib__derived_from_inverse>)
Mn51 ⟶ He4 + V47 [Q = -8.66 MeV] (Mn51 --> He4 + V47 <ths8_reaclib__derived_from_inverse>)
Fe52 ⟶ He4 + Cr48 [Q = -7.94 MeV] (Fe52 --> He4 + Cr48 <ths8_reaclib__derived_from_inverse>)
Fe52 ⟶ p + Mn51 [Q = -7.38 MeV] (Fe52 --> p + Mn51 <ths8_reaclib__derived_from_inverse>)
Co55 + p ⟶ He4 + Fe52 [Q = -0.83 MeV] (Co55 + p --> He4 + Fe52 <ths8_reaclib__derived_from_inverse>)
Co55 ⟶ He4 + Mn51 [Q = -8.21 MeV] (Co55 --> He4 + Mn51 <ths8_reaclib__derived_from_inverse>)
Ni56 ⟶ He4 + Fe52 [Q = -8.00 MeV] (Ni56 --> He4 + Fe52 <ths8_reaclib__derived_from_inverse>)
Ni56 ⟶ p + Co55 [Q = -7.17 MeV] (Ni56 --> p + Co55 <ths8_reaclib__derived_from_inverse>)
Iron group#
Now we’ll add some more nuclei in the iron group.
iron_peak = ["n", "p", "he4",
"mn51",
"fe52", "fe53", "fe54", "fe55", "fe56",
"co55", "co56", "co57",
"ni56", "ni57", "ni58"]
iron_reaclib = reaclib_lib.linking_nuclei(iron_peak)
weak_lib = pyna.TabularLibrary()
iron_weak_lib = weak_lib.linking_nuclei(iron_peak)
warning: He4 was not able to be linked in TabularLibrary
warning: Fe53 was not able to be linked in TabularLibrary
warning: Mn51 was not able to be linked in TabularLibrary
warning: Ni58 was not able to be linked in TabularLibrary
warning: Fe52 was not able to be linked in TabularLibrary
warning: Fe54 was not able to be linked in TabularLibrary
all_lib = core_lib + iron_reaclib + iron_weak_lib
Detailed balance#
Finally, we replace the reverse rates from ReacLib by rederiving them via detailed balance, and including the partition functions.
rates_to_derive = []
for r in all_lib.get_rates():
if r.derived_from_inverse:
# this rate was computed using detailed balance, regardless
# of whether Q < 0 or not. We want to remove it and then
# recompute it
rates_to_derive.append(r)
# now for each of those derived rates, look to see if the pair exists
for r in rates_to_derive:
fr = all_lib.get_rate_by_nuclei(r.products, r.reactants)
if fr:
print(f"modifying {r} from {fr}")
all_lib.remove_rate(r)
d = pyna.DerivedRate(rate=fr, compute_Q=False, use_pf=True)
all_lib.add_rate(d)
modifying N13 ⟶ p + C12 from C12 + p ⟶ N13 + 𝛾
modifying O16 ⟶ He4 + C12 from C12 + He4 ⟶ O16 + 𝛾
modifying F18 ⟶ He4 + N14 from N14 + He4 ⟶ F18 + 𝛾
modifying Ne20 ⟶ He4 + O16 from O16 + He4 ⟶ Ne20 + 𝛾
modifying Na22 ⟶ p + Ne21 from Ne21 + p ⟶ Na22 + 𝛾
modifying Na22 ⟶ He4 + F18 from F18 + He4 ⟶ Na22 + 𝛾
modifying Mg24 ⟶ p + Na23 from Na23 + p ⟶ Mg24 + 𝛾
modifying Mg24 ⟶ He4 + Ne20 from Ne20 + He4 ⟶ Mg24 + 𝛾
modifying Si28 ⟶ p + Al27 from Al27 + p ⟶ Si28 + 𝛾
modifying Si28 ⟶ He4 + Mg24 from Mg24 + He4 ⟶ Si28 + 𝛾
modifying S32 ⟶ p + P31 from P31 + p ⟶ S32 + 𝛾
modifying S32 ⟶ He4 + Si28 from Si28 + He4 ⟶ S32 + 𝛾
modifying Cl35 ⟶ He4 + P31 from P31 + He4 ⟶ Cl35 + 𝛾
modifying Ar36 ⟶ p + Cl35 from Cl35 + p ⟶ Ar36 + 𝛾
modifying Ar36 ⟶ He4 + S32 from S32 + He4 ⟶ Ar36 + 𝛾
modifying K39 ⟶ He4 + Cl35 from Cl35 + He4 ⟶ K39 + 𝛾
modifying Ca40 ⟶ p + K39 from K39 + p ⟶ Ca40 + 𝛾
modifying Ca40 ⟶ He4 + Ar36 from Ar36 + He4 ⟶ Ca40 + 𝛾
modifying Sc43 ⟶ He4 + K39 from K39 + He4 ⟶ Sc43 + 𝛾
modifying Ti44 ⟶ p + Sc43 from Sc43 + p ⟶ Ti44 + 𝛾
modifying Ti44 ⟶ He4 + Ca40 from Ca40 + He4 ⟶ Ti44 + 𝛾
modifying V47 ⟶ He4 + Sc43 from Sc43 + He4 ⟶ V47 + 𝛾
modifying Cr48 ⟶ p + V47 from V47 + p ⟶ Cr48 + 𝛾
modifying Cr48 ⟶ He4 + Ti44 from Ti44 + He4 ⟶ Cr48 + 𝛾
modifying Mn51 ⟶ He4 + V47 from V47 + He4 ⟶ Mn51 + 𝛾
modifying Fe52 ⟶ p + Mn51 from Mn51 + p ⟶ Fe52 + 𝛾
modifying Fe52 ⟶ He4 + Cr48 from Cr48 + He4 ⟶ Fe52 + 𝛾
modifying Co55 ⟶ He4 + Mn51 from Mn51 + He4 ⟶ Co55 + 𝛾
modifying Ni56 ⟶ p + Co55 from Co55 + p ⟶ Ni56 + 𝛾
modifying Ni56 ⟶ He4 + Fe52 from Fe52 + He4 ⟶ Ni56 + 𝛾
modifying C12 ⟶ 3 He4 from 3 He4 ⟶ C12 + 𝛾
modifying O16 + p ⟶ He4 + N13 from N13 + He4 ⟶ p + O16
modifying Ne20 + He4 ⟶ p + Na23 from Na23 + p ⟶ He4 + Ne20
modifying Ne21 + p ⟶ He4 + F18 from F18 + He4 ⟶ p + Ne21
modifying Mg24 + He4 ⟶ p + Al27 from Al27 + p ⟶ He4 + Mg24
modifying Si28 + He4 ⟶ p + P31 from P31 + p ⟶ He4 + Si28
modifying S32 + He4 ⟶ p + Cl35 from Cl35 + p ⟶ He4 + S32
modifying Ar36 + He4 ⟶ p + K39 from K39 + p ⟶ He4 + Ar36
modifying Ca40 + He4 ⟶ p + Sc43 from Sc43 + p ⟶ He4 + Ca40
modifying V47 + p ⟶ He4 + Ti44 from Ti44 + He4 ⟶ p + V47
modifying Mn51 + p ⟶ He4 + Cr48 from Cr48 + He4 ⟶ p + Mn51
modifying Co55 + p ⟶ He4 + Fe52 from Fe52 + He4 ⟶ p + Co55
modifying Fe53 ⟶ n + Fe52 from Fe52 + n ⟶ Fe53 + 𝛾
modifying Fe54 ⟶ n + Fe53 from Fe53 + n ⟶ Fe54 + 𝛾
modifying Fe55 ⟶ n + Fe54 from Fe54 + n ⟶ Fe55 + 𝛾
modifying Fe56 ⟶ n + Fe55 from Fe55 + n ⟶ Fe56 + 𝛾
modifying Co55 ⟶ p + Fe54 from Fe54 + p ⟶ Co55 + 𝛾
modifying Co56 ⟶ n + Co55 from Co55 + n ⟶ Co56 + 𝛾
modifying Co56 ⟶ p + Fe55 from Fe55 + p ⟶ Co56 + 𝛾
modifying Co57 ⟶ n + Co56 from Co56 + n ⟶ Co57 + 𝛾
modifying Co57 ⟶ p + Fe56 from Fe56 + p ⟶ Co57 + 𝛾
modifying Ni57 ⟶ n + Ni56 from Ni56 + n ⟶ Ni57 + 𝛾
modifying Ni57 ⟶ p + Co56 from Co56 + p ⟶ Ni57 + 𝛾
modifying Ni57 ⟶ He4 + Fe53 from Fe53 + He4 ⟶ Ni57 + 𝛾
modifying Ni58 ⟶ n + Ni57 from Ni57 + n ⟶ Ni58 + 𝛾
modifying Ni58 ⟶ p + Co57 from Co57 + p ⟶ Ni58 + 𝛾
modifying Ni58 ⟶ He4 + Fe54 from Fe54 + He4 ⟶ Ni58 + 𝛾
modifying Fe53 + He4 ⟶ n + Ni56 from Ni56 + n ⟶ He4 + Fe53
modifying Fe54 + p ⟶ He4 + Mn51 from Mn51 + He4 ⟶ p + Fe54
modifying Fe54 + He4 ⟶ n + Ni57 from Ni57 + n ⟶ He4 + Fe54
modifying Fe54 + He4 ⟶ p + Co57 from Co57 + p ⟶ He4 + Fe54
modifying Fe55 + p ⟶ n + Co55 from Co55 + n ⟶ p + Fe55
modifying Fe55 + He4 ⟶ n + Ni58 from Ni58 + n ⟶ He4 + Fe55
modifying Fe56 + p ⟶ n + Co56 from Co56 + n ⟶ p + Fe56
modifying Co56 + p ⟶ n + Ni56 from Ni56 + n ⟶ p + Co56
modifying Co56 + p ⟶ He4 + Fe53 from Fe53 + He4 ⟶ p + Co56
modifying Co57 + p ⟶ n + Ni57 from Ni57 + n ⟶ p + Co57
modifying Ni58 + p ⟶ He4 + Co55 from Co55 + He4 ⟶ p + Ni58
Removing duplicates#
There will be some duplicate rates now because we pulled rates both from ReacLib and from tabulated sources. Here we keep the tabulated version of any duplicate rates.
all_lib.eliminate_duplicates()
Creating the network and rate approximations#
Now that we have our library, we can make a network. This can be a RateCollection, PythonNetwork, or AmrexAstroCxxNetwork.
Note
For a SimpleCxxNetwork or FortranNetwork, we do not currently support partition functions, so we would need to comment out the above cell that rederives the reverse rates via detailed balance.
net = pyna.PythonNetwork(libraries=[all_lib], verbose=False)
Next, we will do the \((\alpha,p)(p,\gamma)\) approximation and eliminate the intermediate nuclei
net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"])
net.remove_nuclei(["cl35", "k39", "sc43", "v47"])
We will also approximate some of the neutron captures:
net.make_nn_g_approx(intermediate_nuclei=["fe53", "fe55", "ni57"])
net.remove_nuclei(["fe53", "fe55", "ni57"])
and finally, make some of the protons into NSE protons
# make all rates with A >= 48 use NSE protons
net.make_nse_protons(48)
Visualizing the network#
Let’s visualize the network
comp = pyna.Composition(net.get_nuclei())
comp.set_all(0.1)
comp.set_nuc("he4", 0.95)
comp.normalize()
rho = 1.e6
T = 1.e9
fig = net.plot(rho, T, comp, outfile="subch_simple.png",
rotated=True, hide_xalpha=True, curved_edges=True,
size=(1500, 450),
node_size=550, node_font_size=11)
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N14 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: p_nse partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
Test integration#
We’ll write out the network to a python module and do a test integration with it.
net.write_network("he_burn.py")
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: C12 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N13 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: N14 partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
/opt/hostedtoolcache/Python/3.11.13/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:126: UserWarning: p_nse partition function is not supported by tables: set pf = 1.0 by default
warnings.warn(UserWarning(f'{nuc} partition function is not supported by tables: set pf = 1.0 by default'))
import he_burn
from scipy.integrate import solve_ivp
import numpy as np
rho = 1.e6
T = 3.e9
X0 = np.zeros(he_burn.nnuc)
X0[he_burn.jhe4] = 0.95
X0[he_burn.jc12] = 0.02
Y0 = X0/he_burn.A
from pynucastro.screening import chugunov_2007
tmax = 10.0
sol = solve_ivp(he_burn.rhs, [0, tmax], Y0, method="BDF", jac=he_burn.jacobian,
dense_output=True, args=(rho, T, chugunov_2007),
rtol=1.e-6, atol=1.e-6)
import matplotlib.pyplot as plt
Here we only plot the most abundant species
fig, ax = plt.subplots()
for i in range(he_burn.nnuc):
if sol.y[i, :].max() * he_burn.A[i] > 5.e-3:
ax.loglog(sol.t, sol.y[i,:] * he_burn.A[i],
label=f"X({he_burn.names[i].capitalize()})")
ax.set_ylim(1.e-6, 1.1)
ax.legend(fontsize="small")
ax.set_xlabel("t (s)")
ax.set_ylabel("X")
ax.grid(ls=":")
