mxns figure plot

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()

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()

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()