diff --git a/networks/he-burn/he-burn-31anp/he-burn-31anp.png b/networks/he-burn/he-burn-31anp/he-burn-31anp.png new file mode 100644 index 000000000..c869945c0 Binary files /dev/null and b/networks/he-burn/he-burn-31anp/he-burn-31anp.png differ diff --git a/networks/he-burn/he-burn-31anp/he-burn-31anp.py b/networks/he-burn/he-burn-31anp/he-burn-31anp.py deleted file mode 100644 index d8587c314..000000000 --- a/networks/he-burn/he-burn-31anp/he-burn-31anp.py +++ /dev/null @@ -1,130 +0,0 @@ -import pynucastro as pyna -from pynucastro.rates import ReacLibRate, TabularRate - -DO_DERIVED_RATES = True - -reaclib_lib = pyna.ReacLibLibrary() -weak_lib = pyna.TabularLibrary() - -# these are the nuclei we have in subch_simple -all_reactants = ["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"] - -# create a library of ReacLib rates -core_lib = reaclib_lib.linking_nuclei(all_reactants) - -# in this list, we have the reactants, the actual reactants, -# and modified products that we will use instead - -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.add_rate(_r) - -# finally, the aprox nets don't include the reverse rates for -# C12+C12, C12+O16, and O16+O16, so remove those - -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) - -# C12+Ne20 and reverse -# (a,g) links between Na23 and Al27 -# (a,g) links between Al27 and P31 - -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: - print("removing: ", r) - _r = core_lib.get_rate_by_name(r) - core_lib.remove_rate(_r) - -# now create a list of iron group nuclei and find both the -# ReacLib and weak / tabular rates linking these. - -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) -iron_weak_lib = weak_lib.linking_nuclei(iron_peak) - -# add the libraries - -all_lib = core_lib + iron_reaclib + iron_weak_lib - -if DO_DERIVED_RATES: - 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) - -# we will have duplicate rates -- we want to remove any ReacLib rates -# that we have tabular rates for - -dupes = all_lib.find_duplicate_links() - -rates_to_remove = [] -for d in dupes: - for r in d: - if isinstance(r, ReacLibRate): - rates_to_remove.append(r) - -for r in rates_to_remove: - all_lib.remove_rate(r) - -# combine all three libraries into a single network - -net = pyna.AmrexAstroCxxNetwork(libraries=[all_lib], - symmetric_screening=False) - -# now we approximate some (alpha, p)(p, gamma) links - -net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"]) -net.remove_nuclei(["cl35", "k39", "sc43", "v47"]) - -net.make_nn_g_approx(intermediate_nuclei=["fe53", "fe55", "ni57"]) -net.remove_nuclei(["fe53", "fe55", "ni57"]) - -# make all rates with A >= 48 use NSE protons -net.make_nse_protons(48) - -print(f"number of nuclei = {len(net.unique_nuclei)}") -print(f"number of ReacLib rates = {len(net.reaclib_rates)}") -print(f"number of tabular rates = {len(net.tabular_rates)}") - -fig = net.plot(rotated=True, curved_edges=True, size=(1500, 800), hide_xalpha=True, node_size=400, node_font_size=9) -fig.savefig("newnet.png") - -net.write_network() diff --git a/networks/he-burn/he-burn-31anp/he_burn_31anp.py b/networks/he-burn/he-burn-31anp/he_burn_31anp.py new file mode 100644 index 000000000..6774ab48a --- /dev/null +++ b/networks/he-burn/he-burn-31anp/he_burn_31anp.py @@ -0,0 +1,55 @@ +import pynucastro as pyna +from pynucastro.networks import AmrexAstroCxxNetwork + +import he_burn_core + + +DO_DERIVED_RATES = True + + +def doit(): + + lib = he_burn_core.get_core_library(include_n14_sequence=True, + include_zn=False, + include_iron_peak=True, + include_low_ye=False, + do_detailed_balance=DO_DERIVED_RATES) + + net = pyna.AmrexAstroCxxNetwork(libraries=[lib], + symmetric_screening=False) + + # now we approximate some (alpha, p)(p, gamma) links + + net.make_ap_pg_approx(intermediate_nuclei=["cl35", "k39", "sc43", "v47"]) + net.remove_nuclei(["cl35", "k39", "sc43", "v47"]) + + net.make_nn_g_approx(intermediate_nuclei=["fe53", "fe55", "ni57"]) + net.remove_nuclei(["fe53", "fe55", "ni57"]) + + # make all rates with A >= 48 use NSE protons + net.make_nse_protons(48) + + print(f"number of nuclei = {len(net.unique_nuclei)}") + print(f"number of ReacLib rates = {len(net.reaclib_rates)}") + print(f"number of tabular rates = {len(net.tabular_rates)}") + + # let's make a figure + + comp = pyna.Composition(net.unique_nuclei) + comp.set_equal() + + rho = 9.e7 + T = 6.e9 + + fig = net.plot(rho, T, comp, + rotated=True, curved_edges=True, hide_xalpha=True, + size=(1800, 900), + node_size=500, node_shape="s", node_color="#337dff", node_font_size=10) + + fig.savefig("he-burn-31anp.png") + + net.write_network() + + +if __name__ == "__main__": + doit() diff --git a/networks/he-burn/he-burn-31anp/he_burn_core.py b/networks/he-burn/he-burn-31anp/he_burn_core.py new file mode 120000 index 000000000..bc9cd3928 --- /dev/null +++ b/networks/he-burn/he-burn-31anp/he_burn_core.py @@ -0,0 +1 @@ +../he_burn_core.py \ No newline at end of file diff --git a/networks/he-burn/he-burn-36a/he_burn_36a.py b/networks/he-burn/he-burn-36a/he_burn_36a.py index 4c04aede7..04493ecea 100644 --- a/networks/he-burn/he-burn-36a/he_burn_36a.py +++ b/networks/he-burn/he-burn-36a/he_burn_36a.py @@ -12,6 +12,7 @@ def doit(): lib = he_burn_core.get_core_library(include_n14_sequence=True, include_zn=True, include_iron_peak=True, + include_low_ye=True, do_detailed_balance=DO_DERIVED_RATES) net = pyna.AmrexAstroCxxNetwork(libraries=[lib], diff --git a/networks/he-burn/he_burn_core.py b/networks/he-burn/he_burn_core.py index 7a47e24ca..722d631bf 100644 --- a/networks/he-burn/he_burn_core.py +++ b/networks/he-burn/he_burn_core.py @@ -8,6 +8,7 @@ def get_core_library(*, include_n14_sequence=False, include_zn=False, include_iron_peak=False, + include_low_ye=False, do_detailed_balance=False): reaclib_lib = pyna.ReacLibLibrary() @@ -71,11 +72,17 @@ def get_core_library(*, # ReacLib and weak / tabular rates linking these. iron_peak = ["n", "p", "he4", - "mn51", "mn55", + "mn51", "fe52", "fe53", "fe54", "fe55", "fe56", "co55", "co56", "co57", - "ni56", "ni57", "ni58", - "cu59", "zn60"] + "ni56", "ni57", "ni58"] + + if include_zn: + iron_peak += ["cu59", "zn60"] + + if include_low_ye: + iron_peak += ["mn55"] + iron_reaclib = reaclib_lib.linking_nuclei(iron_peak)