diff --git a/docs/settings_examples/plot_mxns.py b/docs/settings_examples/plot_mxns.py new file mode 100644 index 000000000..3ecd7caa2 --- /dev/null +++ b/docs/settings_examples/plot_mxns.py @@ -0,0 +1,109 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from cosmic.sample.initialbinarytable import InitialBinaryTable +from cosmic.evolve import Evolve + +"""This it the mxns parameter. It sets """ + +n_grid = 30 +masses = np.linspace(5.0, 50.0, n_grid) + +SSEDict = {'stellar_engine': 'sse'} + +BSEDict_base = { + "pts1": 0.001, "pts2": 0.01, "pts3": 0.02, "zsun": 0.014, + "windflag": -1, + "eddlimflag": 0, "neta": 0.5, "bwind": 0.0, "hewind": 0.5, "beta": 0.125, + "xi": 0.5, "acc2": 1.5, "LBV_flag": 1, "alpha1": [1.0, 1.0], + "lambdaf": 0.0, "ceflag": 1, "cekickflag": 2, "cemergeflag": 1, + "cehestarflag": 0, "qcflag": 5, + "qcrit_array": [0.0]*16, + "kickflag": 5, "sigma": 265.0, "bhflag": 1, "bhsigmafrac": 1.0, + "sigmadiv": -20.0, "ecsn": 2.25, "ecsn_mlow": 1.6, "aic": 1, "ussn": 1, + "polar_kick_angle": 90.0, + "natal_kick_array": [[-100.0,-100.0,-100.0,-100.0,0.0],[-100.0,-100.0,-100.0,-100.0,0.0]], + "mm_mu_ns": 400.0, "mm_mu_bh": 200.0, + "remnantflag": 4, + "fryer_mass_limit": 0, + "mxns": 3.0, + "fryer_fmix": 1.0, + "fryer_mcrit_nsbh": 5.75, "rembar_massloss": 0.5, "wd_mass_lim": 1, + "maltsev_mode": 0, "maltsev_fallback": 0.5, "maltsev_pf_prob": 0.1, + "pisn": -2, "ppi_co_shift": 0.0, "ppi_extra_ml": 0.0, "bhspinflag": 0, + "bhspinmag": 0.0, "grflag": 1, "eddfac": 10, "gamma": -2, "don_lim": -1, + "acc_lim": [-1, -1], "smt_periastron_check": 0, "tflag": 1, "ST_tide": 1, + "fprimc_array": [2.0/21.0]*16, + "ifflag": 1, "wdflag": 1, "epsnov": 0.001, "bdecayfac": 1, + "bconst": 3000, "ck": 1000, "rejuv_fac": 1.0, "rejuvflag": 0, + "bhms_coll_flag": 0, "htpmb": 1, "ST_cr": 1, "rtmsflag": 0 +} + +binary_grid = InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.ones(n_grid) * 0.1, + porb=np.ones(n_grid) * 100000.0, + ecc=np.zeros(n_grid), + tphysf=np.ones(n_grid) * 13700.0, + kstar1=np.ones(n_grid), + kstar2=np.ones(n_grid), + metallicity=np.ones(n_grid) * 0.002 +) + +mxns_values = [1.5, 2.0, 3.0, 5.0, 7.0] + +ns_counts = [] +bh_counts = [] + +for mxns in mxns_values: + BSEDict = BSEDict_base.copy() + BSEDict['mxns'] = mxns + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, BSEDict=BSEDict, SSEDict=SSEDict + ) + + ns_count = 0 + bh_count = 0 + for i in range(n_grid): + rows = bpp.loc[i] + final_kstar = rows.iloc[-1]['kstar_1'] + if final_kstar == 13: + ns_count += 1 + elif final_kstar == 14: + bh_count += 1 + + ns_counts.append(ns_count) + bh_counts.append(bh_count) + print(f"mxns={mxns}: {ns_count} NS, {bh_count} BH") + +# bar chart +x = np.arange(len(mxns_values)) +width = 0.35 + +fig, ax = plt.subplots(figsize=(10, 6)) + +bars_ns = ax.bar(x - width/2, ns_counts, width, + label='Neutron Stars', color='steelblue') +bars_bh = ax.bar(x + width/2, bh_counts, width, + label='Black Holes', color='firebrick') + +for bar in bars_ns: + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.2, + str(int(bar.get_height())), ha='center', va='bottom', fontsize=11) +for bar in bars_bh: + ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.2, + str(int(bar.get_height())), ha='center', va='bottom', fontsize=11) + +ax.set_xlabel('mxns (Msun)', fontsize=12) +ax.set_ylabel('Number of Remnants', fontsize=12) +ax.set_title('Effect of mxns on NS vs BH Count\n(higher mxns = more NS, fewer BH)', fontsize=13) +ax.set_xticks(x) +ax.set_xticklabels([f'mxns={m}\n(default)' if m == 3.0 else f'mxns={m}' + for m in mxns_values]) +ax.legend(fontsize=11) +ax.grid(True, alpha=0.3, axis='y') + +plt.tight_layout() +plt.savefig('mxns.png', dpi=150) +plt.show() diff --git a/docs/settings_examples/rembar_massloss.py b/docs/settings_examples/rembar_massloss.py new file mode 100644 index 000000000..0ce84e2fa --- /dev/null +++ b/docs/settings_examples/rembar_massloss.py @@ -0,0 +1,145 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from cosmic.sample.initialbinarytable import InitialBinaryTable +from cosmic.evolve import Evolve + +n_grid = 30 +masses = np.linspace(8.0, 100.0, n_grid) + +SSEDict = {'stellar_engine': 'sse'} + +BSEDict_base = { + "pts1": 0.001, "pts2": 0.01, "pts3": 0.02, "zsun": 0.014, + "windflag": -1, + "eddlimflag": 0, "neta": 0.5, "bwind": 0.0, "hewind": 0.5, "beta": 0.125, + "xi": 0.5, "acc2": 1.5, "LBV_flag": 1, "alpha1": [1.0, 1.0], + "lambdaf": 0.0, "ceflag": 1, "cekickflag": 2, "cemergeflag": 1, + "cehestarflag": 0, "qcflag": 5, + "qcrit_array": [0.0]*16, + "kickflag": 5, "sigma": 265.0, "bhflag": 1, "bhsigmafrac": 1.0, + "sigmadiv": -20.0, "ecsn": 2.25, "ecsn_mlow": 1.6, "aic": 1, "ussn": 1, + "polar_kick_angle": 90.0, + "natal_kick_array": [[-100.0,-100.0,-100.0,-100.0,0.0],[-100.0,-100.0,-100.0,-100.0,0.0]], + "mm_mu_ns": 400.0, "mm_mu_bh": 200.0, + "remnantflag": 4, + "fryer_mass_limit": 0, + "mxns": 3.0, + "fryer_fmix": 1.0, + "fryer_mcrit_nsbh": 5.75, "rembar_massloss": 0.5, + "wd_mass_lim": 1, + "maltsev_mode": 0, "maltsev_fallback": 0.5, "maltsev_pf_prob": 0.1, + "pisn": -2, "ppi_co_shift": 0.0, "ppi_extra_ml": 0.0, "bhspinflag": 0, + "bhspinmag": 0.0, "grflag": 1, "eddfac": 10, "gamma": -2, "don_lim": -1, + "acc_lim": [-1, -1], "smt_periastron_check": 0, "tflag": 1, "ST_tide": 1, + "fprimc_array": [2.0/21.0]*16, + "ifflag": 1, "wdflag": 1, "epsnov": 0.001, "bdecayfac": 1, + "bconst": 3000, "ck": 1000, "rejuv_fac": 1.0, "rejuvflag": 0, + "bhms_coll_flag": 0, "htpmb": 1, "ST_cr": 1, "rtmsflag": 0 +} + +binary_grid = InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.ones(n_grid) * 0.1, + porb=np.ones(n_grid) * 100000.0, + ecc=np.zeros(n_grid), + tphysf=np.ones(n_grid) * 13700.0, + kstar1=np.ones(n_grid), + kstar2=np.ones(n_grid), + metallicity=np.ones(n_grid) * 0.002 +) + +def get_remnant_masses(BSEDict, n_grid): + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, BSEDict=BSEDict, SSEDict=SSEDict + ) + ns_init, ns_final = [], [] + bh_init, bh_final = [], [] + for i in range(n_grid): + rows = bpp.loc[i] + init_mass = rows.iloc[0]['mass_1'] + remnant_rows = rows[rows['kstar_1'].isin([13, 14])] + if len(remnant_rows) > 0: + final_mass = remnant_rows.iloc[-1]['mass_1'] + kstar = remnant_rows.iloc[-1]['kstar_1'] + if kstar == 13: + ns_init.append(init_mass) + ns_final.append(final_mass) + elif kstar == 14: + bh_init.append(init_mass) + bh_final.append(final_mass) + return ns_init, ns_final, bh_init, bh_final + +# two modes — positive and negative +positive_values = { + 0.1: 'red', + 0.5: 'orange', # default + 1.0: 'green', + 2.0: 'blue', +} + +negative_values = { + -0.1: 'red', + -0.3: 'orange', + -0.5: 'green', + -0.7: 'blue', +} + +fig, axes = plt.subplots(2, 2, figsize=(16, 12)) +fig.suptitle('Effect of rembar_massloss on Remnant Masses', fontsize=14) + +# ── top row: positive values ────────────────────────────────────────────────── +for rml, color in positive_values.items(): + BSEDict = BSEDict_base.copy() + BSEDict['rembar_massloss'] = rml + ns_init, ns_final, bh_init, bh_final = get_remnant_masses(BSEDict, n_grid) + label = f'rembar={rml} (default)' if rml == 0.5 else f'rembar={rml}' + + if len(ns_init) > 0: + axes[0, 0].plot(ns_init, ns_final, 'o-', color=color, label=label) + if len(bh_init) > 0: + axes[0, 1].plot(bh_init, bh_final, 'o-', color=color, label=label) + + print(f"rembar_massloss={rml}: {len(ns_init)} NS, {len(bh_init)} BH") + +axes[0, 0].set_xlabel('Initial Stellar Mass (Msun)', fontsize=11) +axes[0, 0].set_ylabel('Final NS Mass (Msun)', fontsize=11) +axes[0, 0].set_title('NS Masses — Positive values\n(higher = more neutrino mass loss)', fontsize=11) +axes[0, 0].legend(fontsize=9) +axes[0, 0].grid(True, alpha=0.3) + +axes[0, 1].set_xlabel('Initial Stellar Mass (Msun)', fontsize=11) +axes[0, 1].set_ylabel('Final BH Mass (Msun)', fontsize=11) +axes[0, 1].set_title('BH Masses — Positive values\n(higher = more neutrino mass loss)', fontsize=11) +axes[0, 1].legend(fontsize=9) +axes[0, 1].grid(True, alpha=0.3) + +# ── bottom row: negative values ─────────────────────────────────────────────── +for rml, color in negative_values.items(): + BSEDict = BSEDict_base.copy() + BSEDict['rembar_massloss'] = rml + ns_init, ns_final, bh_init, bh_final = get_remnant_masses(BSEDict, n_grid) + label = f'rembar={rml}' + + if len(ns_init) > 0: + axes[1, 0].plot(ns_init, ns_final, 'o-', color=color, label=label) + if len(bh_init) > 0: + axes[1, 1].plot(bh_init, bh_final, 'o-', color=color, label=label) + + print(f"rembar_massloss={rml}: {len(ns_init)} NS, {len(bh_init)} BH") + +axes[1, 0].set_xlabel('Initial Stellar Mass (Msun)', fontsize=11) +axes[1, 0].set_ylabel('Final NS Mass (Msun)', fontsize=11) +axes[1, 0].set_title('NS Masses — Negative values\n(more negative = less mass kept)', fontsize=11) +axes[1, 0].legend(fontsize=9) +axes[1, 0].grid(True, alpha=0.3) + +axes[1, 1].set_xlabel('Initial Stellar Mass (Msun)', fontsize=11) +axes[1, 1].set_ylabel('Final BH Mass (Msun)', fontsize=11) +axes[1, 1].set_title('BH Masses — Negative values\n(more negative = less mass kept)', fontsize=11) +axes[1, 1].legend(fontsize=9) +axes[1, 1].grid(True, alpha=0.3) + +plt.tight_layout() +plt.savefig('rembar_massloss.png', dpi=150) +plt.show() diff --git a/docs/settings_examples/remnant_flag.py b/docs/settings_examples/remnant_flag.py new file mode 100644 index 000000000..4dbc01a30 --- /dev/null +++ b/docs/settings_examples/remnant_flag.py @@ -0,0 +1,115 @@ +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from cosmic.sample.initialbinarytable import InitialBinaryTable +from cosmic.evolve import Evolve + +n_grid = 30 +masses = np.linspace(5.0, 100.0, n_grid) # extended down to 5 Msun + +SSEDict = {'stellar_engine': 'sse'} + +BSEDict_base = { + "pts1": 0.001, "pts2": 0.01, "pts3": 0.02, "zsun": 0.014, + "windflag": -1, + "eddlimflag": 0, "neta": 0.5, "bwind": 0.0, "hewind": 0.5, "beta": 0.125, + "xi": 0.5, "acc2": 1.5, "LBV_flag": 1, "alpha1": [1.0, 1.0], + "lambdaf": 0.0, "ceflag": 1, "cekickflag": 2, "cemergeflag": 1, + "cehestarflag": 0, "qcflag": 5, + "qcrit_array": [0.0]*16, + "kickflag": 5, "sigma": 265.0, "bhflag": 1, "bhsigmafrac": 1.0, + "sigmadiv": -20.0, "ecsn": 2.25, "ecsn_mlow": 1.6, "aic": 1, "ussn": 1, + "polar_kick_angle": 90.0, + "natal_kick_array": [[-100.0,-100.0,-100.0,-100.0,0.0],[-100.0,-100.0,-100.0,-100.0,0.0]], + "mm_mu_ns": 400.0, "mm_mu_bh": 200.0, + "remnantflag": 4, + "fryer_mass_limit": 0, + "mxns": 3.0, "fryer_fmix": 1.0, + "fryer_mcrit_nsbh": 5.75, "rembar_massloss": 0.5, "wd_mass_lim": 1, + "maltsev_mode": 0, "maltsev_fallback": 0.5, "maltsev_pf_prob": 0.1, + "pisn": -2, "ppi_co_shift": 0.0, "ppi_extra_ml": 0.0, "bhspinflag": 0, + "bhspinmag": 0.0, "grflag": 1, "eddfac": 10, "gamma": -2, "don_lim": -1, + "acc_lim": [-1, -1], "smt_periastron_check": 0, "tflag": 1, "ST_tide": 1, + "fprimc_array": [2.0/21.0]*16, + "ifflag": 1, "wdflag": 1, "epsnov": 0.001, "bdecayfac": 1, + "bconst": 3000, "ck": 1000, "rejuv_fac": 1.0, "rejuvflag": 0, + "bhms_coll_flag": 0, "htpmb": 1, "ST_cr": 1, "rtmsflag": 0 +} + +remnant_flags = { + 0: ('Hurley+2000', 'blue', 0.002), # use low metallicity + 1: ('Belczynski+2002', 'orange', 0.002), + 2: ('Belczynski+2008', 'green', 0.002), + 3: ('Fryer+2012 rapid', 'red', 0.002), + 4: ('Fryer+2012 delayed', 'purple', 0.002), + 5: ('Mandel & Mueller 2020', 'brown', 0.002), + 6: ('Maltsev+2025', 'pink', 0.02), # use solar metallicity + 7: ('Fryer+2022', 'gray', 0.002), +} + +fig, axes = plt.subplots(4, 2, figsize=(14, 20)) +fig.suptitle('Effect of remnantflag on NS and BH Masses', fontsize=16) +axes = axes.flatten() + +for idx, (rflag, (label, color, metallicity)) in enumerate(remnant_flags.items()): + ax = axes[idx] + + binary_grid = InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.ones(n_grid) * 0.1, + porb=np.ones(n_grid) * 100000.0, + ecc=np.zeros(n_grid), + tphysf=np.ones(n_grid) * 13700.0, + kstar1=np.ones(n_grid), + kstar2=np.ones(n_grid), + metallicity=np.ones(n_grid) * metallicity + ) + + BSEDict = BSEDict_base.copy() + BSEDict['remnantflag'] = rflag + + # remnantflag=5 needs specific settings per docs + if rflag == 5: + BSEDict['mxns'] = 2.0 + BSEDict['rembar_massloss'] = 0.1 + BSEDict['kickflag'] = 6 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, BSEDict=BSEDict, SSEDict=SSEDict + ) + + init_masses = [] + final_masses = [] + kstar_types = [] + for i in range(n_grid): + rows = bpp.loc[i] + init_mass = rows.iloc[0]['mass_1'] + remnant_rows = rows[rows['kstar_1'].isin([13, 14])] + if len(remnant_rows) > 0: + final_mass = remnant_rows.iloc[-1]['mass_1'] + kstar = remnant_rows.iloc[-1]['kstar_1'] + init_masses.append(init_mass) + final_masses.append(final_mass) + kstar_types.append(kstar) + + ns_init = [init_masses[i] for i in range(len(kstar_types)) if kstar_types[i] == 13] + ns_final = [final_masses[i] for i in range(len(kstar_types)) if kstar_types[i] == 13] + bh_init = [init_masses[i] for i in range(len(kstar_types)) if kstar_types[i] == 14] + bh_final = [final_masses[i] for i in range(len(kstar_types)) if kstar_types[i] == 14] + + if len(ns_init) > 0: + ax.plot(ns_init, ns_final, 'o-', color='steelblue', label='Neutron Star') + if len(bh_init) > 0: + ax.plot(bh_init, bh_final, 'o-', color=color, label='Black Hole') + + ax.axhspan(2, 5, alpha=0.1, color='gray', label='Mass gap (2-5 Msun)') + ax.set_xlabel('Initial Stellar Mass (Msun)', fontsize=10) + ax.set_ylabel('Final Remnant Mass (Msun)', fontsize=10) + ax.set_title(f'remnantflag={rflag}: {label} (Z={metallicity})', fontsize=11) + ax.legend(fontsize=8) + ax.grid(True, alpha=0.3) + print(f"remnantflag={rflag} ({label}): {len(ns_init)} NS, {len(bh_init)} BH found") + +plt.tight_layout() +plt.savefig('remnantflag.png', dpi=150) +plot.show()