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
He4 + He4 + He4 ⟶ C12 + 𝛾 [Q = 7.28 MeV] (He4 + He4 + 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__reverse>)
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__reverse>)
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__reverse>)
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 ⟶ He4 + He4 + He4 [Q = -7.28 MeV] (C12 --> He4 + He4 + He4 <fy05_reaclib__reverse>)
N13 ⟶ p + C12 [Q = -1.94 MeV] (N13 --> p + C12 <ls09_reaclib__reverse>)
O16 + p ⟶ He4 + N13 [Q = -5.22 MeV] (O16 + p --> He4 + N13 <cf88_reaclib__reverse>)
O16 ⟶ He4 + C12 [Q = -7.16 MeV] (O16 --> He4 + C12 <nac2_reaclib__reverse>)
F18 ⟶ He4 + N14 [Q = -4.41 MeV] (F18 --> He4 + N14 <il10_reaclib__reverse>)
Ne20 + He4 ⟶ p + Na23 [Q = -2.38 MeV] (Ne20 + He4 --> p + Na23 <il10_reaclib__reverse>)
Ne20 ⟶ He4 + O16 [Q = -4.73 MeV] (Ne20 --> He4 + O16 <co10_reaclib__reverse>)
Ne21 + p ⟶ He4 + F18 [Q = -1.74 MeV] (Ne21 + p --> He4 + F18 <rpsm_reaclib__reverse>)
Na22 ⟶ He4 + F18 [Q = -8.48 MeV] (Na22 --> He4 + F18 <rpsm_reaclib__reverse>)
Na22 ⟶ p + Ne21 [Q = -6.74 MeV] (Na22 --> p + Ne21 <il10_reaclib__reverse>)
Mg24 + He4 ⟶ p + Al27 [Q = -1.60 MeV] (Mg24 + He4 --> p + Al27 <il10_reaclib__reverse>)
Mg24 ⟶ He4 + Ne20 [Q = -9.32 MeV] (Mg24 --> He4 + Ne20 <il10_reaclib__reverse>)
Mg24 ⟶ p + Na23 [Q = -11.69 MeV] (Mg24 --> p + Na23 <il10_reaclib__reverse>)
Si28 + He4 ⟶ p + P31 [Q = -1.92 MeV] (Si28 + He4 --> p + P31 <il10_reaclib__reverse>)
Si28 ⟶ He4 + Mg24 [Q = -9.98 MeV] (Si28 --> He4 + Mg24 <st08_reaclib__reverse>)
Si28 ⟶ p + Al27 [Q = -11.59 MeV] (Si28 --> p + Al27 <il10_reaclib__reverse>)
S32 + He4 ⟶ p + Cl35 [Q = -1.87 MeV] (S32 + He4 --> p + Cl35 <il10_reaclib__reverse>)
S32 ⟶ He4 + Si28 [Q = -6.95 MeV] (S32 --> He4 + Si28 <ths8_reaclib__reverse>)
S32 ⟶ p + P31 [Q = -8.86 MeV] (S32 --> p + P31 <il10_reaclib__reverse>)
Cl35 ⟶ He4 + P31 [Q = -7.00 MeV] (Cl35 --> He4 + P31 <ths8_reaclib__reverse>)
Ar36 + He4 ⟶ p + K39 [Q = -1.29 MeV] (Ar36 + He4 --> p + K39 <ths8_reaclib__reverse>)
Ar36 ⟶ He4 + S32 [Q = -6.64 MeV] (Ar36 --> He4 + S32 <ths8_reaclib__reverse>)
Ar36 ⟶ p + Cl35 [Q = -8.51 MeV] (Ar36 --> p + Cl35 <il10_reaclib__reverse>)
K39 ⟶ He4 + Cl35 [Q = -7.22 MeV] (K39 --> He4 + Cl35 <ths8_reaclib__reverse>)
Ca40 + He4 ⟶ p + Sc43 [Q = -3.52 MeV] (Ca40 + He4 --> p + Sc43 <ths8_reaclib__reverse>)
Ca40 ⟶ He4 + Ar36 [Q = -7.04 MeV] (Ca40 --> He4 + Ar36 <ths8_reaclib__reverse>)
Ca40 ⟶ p + K39 [Q = -8.33 MeV] (Ca40 --> p + K39 <lo18_reaclib__reverse>)
Sc43 ⟶ He4 + K39 [Q = -4.81 MeV] (Sc43 --> He4 + K39 <ths8_reaclib__reverse>)
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__reverse>)
Ti44 ⟶ p + Sc43 [Q = -8.65 MeV] (Ti44 --> p + Sc43 <ths8_reaclib__reverse>)
V47 ⟶ He4 + Sc43 [Q = -8.24 MeV] (V47 --> He4 + Sc43 <ths8_reaclib__reverse>)
Cr48 ⟶ He4 + Ti44 [Q = -7.70 MeV] (Cr48 --> He4 + Ti44 <ths8_reaclib__reverse>)
Cr48 ⟶ p + V47 [Q = -8.11 MeV] (Cr48 --> p + V47 <nfis_reaclib__reverse>)
Mn51 + p ⟶ He4 + Cr48 [Q = -0.56 MeV] (Mn51 + p --> He4 + Cr48 <ths8_reaclib__reverse>)
Mn51 ⟶ He4 + V47 [Q = -8.66 MeV] (Mn51 --> He4 + V47 <ths8_reaclib__reverse>)
Fe52 ⟶ He4 + Cr48 [Q = -7.94 MeV] (Fe52 --> He4 + Cr48 <ths8_reaclib__reverse>)
Fe52 ⟶ p + Mn51 [Q = -7.38 MeV] (Fe52 --> p + Mn51 <ths8_reaclib__reverse>)
Co55 + p ⟶ He4 + Fe52 [Q = -0.83 MeV] (Co55 + p --> He4 + Fe52 <ths8_reaclib__reverse>)
Co55 ⟶ He4 + Mn51 [Q = -8.21 MeV] (Co55 --> He4 + Mn51 <ths8_reaclib__reverse>)
Ni56 ⟶ He4 + Fe52 [Q = -8.00 MeV] (Ni56 --> He4 + Fe52 <ths8_reaclib__reverse>)
Ni56 ⟶ p + Co55 [Q = -7.17 MeV] (Ni56 --> p + Co55 <ths8_reaclib__reverse>)
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
warning: Fe53 was not able to be linked
warning: Mn51 was not able to be linked
warning: Ni58 was not able to be linked
warning: Fe52 was not able to be linked
warning: Fe54 was not able to be linked
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.reverse:
# 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 ⟶ He4 + He4 + He4 from He4 + He4 + 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.
dupes = all_lib.find_duplicate_links()
rates_to_remove = []
for d in dupes:
for r in d:
if isinstance(r, pyna.rates.ReacLibRate):
rates_to_remove.append(r)
for r in rates_to_remove:
all_lib.remove_rate(r)
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])
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"])
using approximate rate S32 + He4 ⟶ Ar36 + 𝛾
using approximate rate Ar36 ⟶ S32 + He4
using approximate rate Ar36 + He4 ⟶ Ca40 + 𝛾
using approximate rate Ca40 ⟶ Ar36 + He4
using approximate rate Ca40 + He4 ⟶ Ti44 + 𝛾
using approximate rate Ti44 ⟶ Ca40 + He4
using approximate rate Ti44 + He4 ⟶ Cr48 + 𝛾
using approximate rate Cr48 ⟶ Ti44 + He4
removing rate S32 + He4 ⟶ Ar36 + 𝛾
removing rate S32 + He4 ⟶ p + Cl35
removing rate Cl35 + p ⟶ Ar36 + 𝛾
removing rate Ar36 ⟶ He4 + S32
removing rate Ar36 ⟶ p + Cl35
removing rate Cl35 + p ⟶ He4 + S32
removing rate Ar36 + He4 ⟶ Ca40 + 𝛾
removing rate Ar36 + He4 ⟶ p + K39
removing rate K39 + p ⟶ Ca40 + 𝛾
removing rate Ca40 ⟶ He4 + Ar36
removing rate Ca40 ⟶ p + K39
removing rate K39 + p ⟶ He4 + Ar36
removing rate Ca40 + He4 ⟶ Ti44 + 𝛾
removing rate Ca40 + He4 ⟶ p + Sc43
removing rate Sc43 + p ⟶ Ti44 + 𝛾
removing rate Ti44 ⟶ He4 + Ca40
removing rate Ti44 ⟶ p + Sc43
removing rate Sc43 + p ⟶ He4 + Ca40
removing rate Ti44 + He4 ⟶ Cr48 + 𝛾
removing rate Ti44 + He4 ⟶ p + V47
removing rate V47 + p ⟶ Cr48 + 𝛾
removing rate Cr48 ⟶ He4 + Ti44
removing rate Cr48 ⟶ p + V47
removing rate V47 + p ⟶ He4 + Ti44
looking to remove P31 + He4 ⟶ Cl35 + 𝛾
looking to remove Cl35 + He4 ⟶ K39 + 𝛾
looking to remove Cl35 ⟶ He4 + P31
looking to remove K39 ⟶ He4 + Cl35
looking to remove Cl35 + He4 ⟶ K39 + 𝛾
looking to remove K39 + He4 ⟶ Sc43 + 𝛾
looking to remove K39 ⟶ He4 + Cl35
looking to remove Sc43 ⟶ He4 + K39
looking to remove K39 + He4 ⟶ Sc43 + 𝛾
looking to remove Sc43 + He4 ⟶ V47 + 𝛾
looking to remove Sc43 ⟶ He4 + K39
looking to remove V47 ⟶ He4 + Sc43
looking to remove Sc43 + He4 ⟶ V47 + 𝛾
looking to remove V47 + He4 ⟶ Mn51 + 𝛾
looking to remove V47 ⟶ He4 + Sc43
looking to remove Mn51 ⟶ He4 + 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"])
approximating out Fe53
using approximate rate Fe52 + n + n ⟶ Fe54 + 𝛾
using approximate rate Fe54 ⟶ Fe52 + n + n
approximating out Fe55
using approximate rate Fe54 + n + n ⟶ Fe56 + 𝛾
using approximate rate Fe56 ⟶ Fe54 + n + n
approximating out Ni57
using approximate rate Ni56 + n + n ⟶ Ni58 + 𝛾
using approximate rate Ni58 ⟶ Ni56 + n + n
removing rate Fe52 + n ⟶ Fe53 + 𝛾
removing rate Fe53 + n ⟶ Fe54 + 𝛾
removing rate Fe54 ⟶ n + Fe53
removing rate Fe53 ⟶ n + Fe52
removing rate Fe54 + n ⟶ Fe55 + 𝛾
removing rate Fe55 + n ⟶ Fe56 + 𝛾
removing rate Fe56 ⟶ n + Fe55
removing rate Fe55 ⟶ n + Fe54
removing rate Ni56 + n ⟶ Ni57 + 𝛾
removing rate Ni57 + n ⟶ Ni58 + 𝛾
removing rate Ni58 ⟶ n + Ni57
removing rate Ni57 ⟶ n + Ni56
looking to remove Fe53 + He4 ⟶ Ni57 + 𝛾
looking to remove Fe53 + He4 ⟶ p + Co56
looking to remove Ni56 + n ⟶ He4 + Fe53
looking to remove Ni57 ⟶ He4 + Fe53
looking to remove Fe53 + He4 ⟶ n + Ni56
looking to remove Co56 + p ⟶ He4 + Fe53
looking to remove Fe55 + p ⟶ Co56 + 𝛾
looking to remove Co55 + n ⟶ p + Fe55
looking to remove Ni58 + n ⟶ He4 + Fe55
looking to remove Co56 ⟶ p + Fe55
looking to remove Fe55 + p ⟶ n + Co55
looking to remove Fe55 + He4 ⟶ n + Ni58
looking to remove Co55 + e⁻ ⟶ Fe55 + 𝜈
looking to remove Fe55 ⟶ Co55 + e⁻ + 𝜈
looking to remove Fe53 + He4 ⟶ Ni57 + 𝛾
looking to remove Co56 + p ⟶ Ni57 + 𝛾
looking to remove Ni57 + n ⟶ p + Co57
looking to remove Ni57 + n ⟶ He4 + Fe54
looking to remove Ni57 ⟶ p + Co56
looking to remove Ni57 ⟶ He4 + Fe53
looking to remove Fe54 + He4 ⟶ n + Ni57
looking to remove Co57 + p ⟶ n + Ni57
looking to remove Co57 ⟶ Ni57 + e⁻ + 𝜈
looking to remove Ni57 + e⁻ ⟶ Co57 + 𝜈
and finally, make some of the protons into NSE protons
# make all rates with A >= 48 use NSE protons
net.make_nse_protons(48)
modifying p_Mn51__Fe52 to use NSE protons
modifying Fe52__p_Mn51__derived to use NSE protons
modifying p_Co55__Ni56 to use NSE protons
modifying Ni56__p_Co55__derived to use NSE protons
modifying He4_Cr48__p_Mn51 to use NSE protons
modifying p_Mn51__He4_Cr48__derived to use NSE protons
modifying He4_Fe52__p_Co55 to use NSE protons
modifying p_Co55__He4_Fe52__derived to use NSE protons
modifying p_Fe54__Co55 to use NSE protons
modifying Co55__p_Fe54__derived to use NSE protons
modifying p_Fe56__Co57 to use NSE protons
modifying Co57__p_Fe56__derived to use NSE protons
modifying p_Co57__Ni58 to use NSE protons
modifying Ni58__p_Co57__derived to use NSE protons
modifying He4_Mn51__p_Fe54 to use NSE protons
modifying p_Fe54__He4_Mn51__derived to use NSE protons
modifying He4_Co55__p_Ni58 to use NSE protons
modifying p_Ni58__He4_Co55__derived to use NSE protons
modifying n_Co56__p_Fe56 to use NSE protons
modifying p_Fe56__n_Co56__derived to use NSE protons
modifying p_Co57__He4_Fe54 to use NSE protons
modifying He4_Fe54__p_Co57__derived to use NSE protons
modifying n_Ni56__p_Co56 to use NSE protons
modifying p_Co56__n_Ni56__derived to use NSE protons
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=500, node_font_size=11, node_color="#337dff", node_shape="s")
/opt/hostedtoolcache/Python/3.11.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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.11/x64/lib/python3.11/site-packages/pynucastro/rates/derived_rate.py:86: 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
fig, ax = plt.subplots()
for i in range(he_burn.nnuc):
if sol.y[i, :].max() > 5.e-4:
ax.loglog(sol.t, sol.y[i,:] * he_burn.A[i], label=f"X({he_burn.names[i].capitalize()})")
ax.set_ylim(1.e-8, 1.1)
ax.legend(fontsize="small")
ax.set_xlabel("t (s)")
ax.set_ylabel("X")
Text(0, 0.5, 'X')
