From 28eab6232ee3a773b61be0c7a00368ec0b665e0f Mon Sep 17 00:00:00 2001 From: lukewilner Date: Wed, 24 Jun 2026 14:27:14 -0700 Subject: [PATCH 1/3] 5 new config. figures --- docs/settings_examples/plot_hewind.py | 133 +++++++++++++++++++++ docs/settings_examples/plot_neta.py | 134 ++++++++++++++++++++++ docs/settings_examples/plot_pts.py | 140 +++++++++++++++++++++++ docs/settings_examples/plot_windflag.py | 146 ++++++++++++++++++++++++ docs/settings_examples/plot_zsun.py | 141 +++++++++++++++++++++++ 5 files changed, 694 insertions(+) create mode 100644 docs/settings_examples/plot_hewind.py create mode 100644 docs/settings_examples/plot_neta.py create mode 100644 docs/settings_examples/plot_pts.py create mode 100644 docs/settings_examples/plot_windflag.py create mode 100644 docs/settings_examples/plot_zsun.py diff --git a/docs/settings_examples/plot_hewind.py b/docs/settings_examples/plot_hewind.py new file mode 100644 index 000000000..a29aef579 --- /dev/null +++ b/docs/settings_examples/plot_hewind.py @@ -0,0 +1,133 @@ +""" +``hewind`` +========== + +This example shows the effect of the ``hewind`` setting on naked helium-star +evolution. ``hewind`` scales the wind mass loss for naked helium stars, so this +test initializes the primary as a helium main-sequence star with ``kstar_1 = 7``. + +The binaries are very wide, which keeps the example focused on helium-star wind +mass loss rather than Roche-lobe overflow or common-envelope evolution. The +notebook version of this test found the clearest direct ``hewind`` dependence +with ``windflag = 0``, so this script uses that setup for the final diagnostic. +""" + +import sys + +sys.path.append("..") +import generate_default_bsedict + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) +plt.style.use("../_static/gallery.mplstyle") + + +hewind_values = [0.0, 0.25, 0.5, 0.75, 1.0] +he_masses = np.array([3.0, 5.0, 8.0, 12.0, 20.0]) + + +def hewind_label(hewind): + """Format labels so the no-hewind and default cases are explicit.""" + if hewind == 0.0: + return "0 (none)" + if hewind == 0.5: + return "0.5 (default)" + return f"{hewind:g}" + + +def make_helium_star_grid(): + """Build wide binaries with naked helium-star primaries.""" + n_systems = len(he_masses) + return InitialBinaryTable.InitialBinaries( + m1=he_masses, + m2=np.ones(n_systems), + porb=np.ones(n_systems) * 1.0e6, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * 2000.0, + kstar1=np.ones(n_systems) * 7, + kstar2=np.ones(n_systems), + metallicity=np.ones(n_systems) * 0.014, + ) + + +def summarize_final_state(bcm, hewind): + """Collect final masses for one hewind value.""" + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["hewind"] = hewind + summary["initial_he_mass_1"] = he_masses[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + summary["mass_lost_1"] = summary["initial_he_mass_1"] - summary["final_mass_1"] + summary["final_kstar_1"] = final_bcm["kstar_1"] + return summary.reset_index(drop=True) + + +def evolve_hewind_grid(helium_grid, windflag=0, dtp=2.0): + """Run the helium-star grid for each hewind value.""" + summaries = [] + + for hewind in hewind_values: + settings = BSEDict.copy() + settings["windflag"] = windflag + settings["hewind"] = hewind + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=helium_grid, + BSEDict=settings, + dtp=dtp, + ) + summaries.append(summarize_final_state(bcm, hewind)) + + return pd.concat(summaries, ignore_index=True) + + +helium_grid = make_helium_star_grid() +summaries = evolve_hewind_grid(helium_grid, windflag=0, dtp=2.0) +baseline = summaries[summaries["hewind"].eq(0.5)].set_index("initial_he_mass_1") + +fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), sharex=True, constrained_layout=True) + +for hewind in hewind_values: + subset = summaries[summaries["hewind"].eq(hewind)].sort_values( + "initial_he_mass_1" + ) + label = hewind_label(hewind) + + axes[0].plot( + subset["initial_he_mass_1"], + subset["final_mass_1"], + marker="o", + label=label, + ) + + indexed = subset.set_index("initial_he_mass_1") + common = indexed.index.intersection(baseline.index) + delta = indexed.loc[common, "final_mass_1"] - baseline.loc[ + common, "final_mass_1" + ] + axes[1].plot(common, delta, marker="o", label=label) + +axes[1].axhline(0, linestyle="--", linewidth=1, color="black") +axes[0].set_xlabel("Initial He-star mass [$M_\\odot$]") +axes[0].set_ylabel("Final primary mass [$M_\\odot$]") +axes[0].set_title("Final mass") + +axes[1].set_xlabel("Initial He-star mass [$M_\\odot$]") +axes[1].set_ylabel("Delta final mass [$M_\\odot$]") +axes[1].set_title("Difference from default") + +for ax in axes: + ax.legend(title="hewind", fontsize=8, title_fontsize=9, markerscale=1.2) + ax.grid(True, alpha=0.3) + +fig.suptitle("Effect of hewind on naked helium stars, windflag=0") +plt.savefig("hewind_flagtest.png", dpi=150) +plt.show() diff --git a/docs/settings_examples/plot_neta.py b/docs/settings_examples/plot_neta.py new file mode 100644 index 000000000..007e442ce --- /dev/null +++ b/docs/settings_examples/plot_neta.py @@ -0,0 +1,134 @@ +""" +``neta`` +======== + +This example shows the effect of the ``neta`` setting on low- and +intermediate-mass single-star evolution. ``neta`` is the Reimers wind +coefficient, so it mainly affects cool evolved stars on the red giant branch +and asymptotic giant branch. + +COSMIC requires positive ``neta`` values, so ``neta = 0.01`` is used as a +nearly-off control. The left panel focuses on the narrow mass range where the +final mass is most sensitive to ``neta``. The right panel shows the same final +mass comparison over the full initial-mass grid. +""" + +import sys + +sys.path.append("..") +import generate_default_bsedict + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) +plt.style.use("../_static/gallery.mplstyle") + + +neta_values = [0.01, 0.25, 0.5, 1.0, 2.0] +masses = np.round(np.arange(0.8, 8.05, 0.1), 2) + + +def make_single_star_grid(): + """Build the low/intermediate-mass single-star grid.""" + n_systems = len(masses) + return InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.zeros(n_systems), + porb=np.ones(n_systems) * -1.0, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * 13700.0, + kstar1=np.ones(n_systems), + kstar2=np.zeros(n_systems), + metallicity=np.ones(n_systems) * 0.014, + ) + + +def neta_label(neta): + """Format labels so the default and nearly-off cases are explicit.""" + if neta == 0.01: + return "0.01 (near off)" + if neta == 0.5: + return "0.5 (default)" + return f"{neta:g}" + + +def summarize_final_state(bcm, neta): + """Collect final masses for one neta value.""" + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["neta"] = neta + summary["initial_mass_1"] = masses[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + summary["mass_lost_1"] = summary["initial_mass_1"] - summary["final_mass_1"] + summary["final_kstar_1"] = final_bcm["kstar_1"] + return summary.reset_index(drop=True) + + +def evolve_neta_grid(grid): + """Run the single-star grid for each neta value.""" + summaries = [] + + for neta in neta_values: + settings = BSEDict.copy() + settings["neta"] = neta + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=grid, + BSEDict=settings, + dtp=25, + ) + summaries.append(summarize_final_state(bcm, neta)) + + return pd.concat(summaries, ignore_index=True) + + +def plot_final_mass(ax, data, title): + """Plot final mass as a function of initial mass for all neta values.""" + for neta in neta_values: + subset = data[data["neta"].eq(neta)].sort_values("initial_mass_1") + ax.plot( + subset["initial_mass_1"], + subset["final_mass_1"], + marker="o", + markersize=4, + label=neta_label(neta), + ) + + ax.set_xlabel("Initial mass [$M_\\odot$]") + ax.set_ylabel("Final mass [$M_\\odot$]") + ax.set_title(title) + ax.legend(title="neta", fontsize=8, title_fontsize=9, markerscale=1.2) + ax.grid(True, alpha=0.3) + + +grid = make_single_star_grid() +summaries = evolve_neta_grid(grid) +zoom = summaries[ + (summaries["initial_mass_1"] >= 0.9) + & (summaries["initial_mass_1"] <= 1.5) +] + +fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) + +plot_final_mass( + axes[0], + zoom, + "Sensitive transition region", +) +plot_final_mass( + axes[1], + summaries, + "Full mass grid", +) + +fig.suptitle("Effect of Reimers wind coefficient neta") +plt.tight_layout() +plt.savefig("neta_flagtest.png", dpi=150) +plt.show() diff --git a/docs/settings_examples/plot_pts.py b/docs/settings_examples/plot_pts.py new file mode 100644 index 000000000..a682ebfbc --- /dev/null +++ b/docs/settings_examples/plot_pts.py @@ -0,0 +1,140 @@ +""" +``pts1``, ``pts2``, and ``pts3`` +================================ + +This example shows how the ``pts1``, ``pts2``, and ``pts3`` timestep modifiers +affect single-star evolution in COSMIC. These are numerical timestep-resolution +parameters, not physical prescriptions. Smaller values use finer timesteps. + +The broad grid compares final stellar masses from 1 to 50 solar masses for +three timestep choices. The zoomed grid repeats the comparison from 22 to 25 +solar masses and includes one extra-fine timestep choice to highlight a more +timestep-sensitive mass range. +""" + +import sys + +sys.path.append("..") +import generate_default_bsedict + +import matplotlib.pyplot as plt +import numpy as np + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) +plt.style.use("../_static/gallery.mplstyle") + + +N_GRID = 50 +TPHYSF = 13700.0 +METALLICITY = 0.014 + +masses = np.linspace(1.0, 50.0, N_GRID) +zoom_masses = np.linspace(22.0, 25.0, 31) + +pts_values = [ + {"pts1": 0.0005, "pts2": 0.005, "pts3": 0.01}, + {"pts1": 0.001, "pts2": 0.01, "pts3": 0.02}, + {"pts1": 0.002, "pts2": 0.02, "pts3": 0.04}, +] + +zoom_pts_values = [ + {"pts1": 0.00025, "pts2": 0.0025, "pts3": 0.005}, + {"pts1": 0.0005, "pts2": 0.005, "pts3": 0.01}, + {"pts1": 0.001, "pts2": 0.01, "pts3": 0.02}, + {"pts1": 0.002, "pts2": 0.02, "pts3": 0.04}, +] + + +def make_single_star_grid(mass_grid): + """Build a COSMIC single-star grid using the standard single-star setup.""" + n_systems = len(mass_grid) + return InitialBinaryTable.InitialBinaries( + m1=mass_grid, + m2=np.zeros(n_systems), + porb=np.ones(n_systems) * -1.0, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * TPHYSF, + kstar1=np.ones(n_systems), + kstar2=np.zeros(n_systems), + metallicity=np.ones(n_systems) * METALLICITY, + ) + + +def pts_label(pts): + """Format a legend label from the three timestep modifiers.""" + return f"pts=({pts['pts1']}, {pts['pts2']}, {pts['pts3']})" + + +def extract_final_mass_1(bcm, n_systems): + """Return final primary masses, ordered by COSMIC binary number.""" + final = ( + bcm.sort_values(["bin_num", "tphys"]) + .groupby("bin_num")["mass_1"] + .last() + .reindex(np.arange(n_systems)) + ) + + if final.isna().any(): + missing = final[final.isna()].index.to_list() + raise ValueError(f"Missing final bcm rows for bin_num values: {missing}") + + return final.to_numpy() + + +def final_masses_for_pts(mass_grid, pts_grid): + """Evolve one mass grid for each pts choice and collect final masses.""" + binary_grid = make_single_star_grid(mass_grid) + final_masses = {} + + for pts in pts_grid: + settings = BSEDict.copy() + settings.update(pts) + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, + BSEDict=settings, + dtp=0, + ) + + final_masses[pts_label(pts)] = extract_final_mass_1(bcm, len(mass_grid)) + + return final_masses + + +def plot_final_mass_panel(ax, mass_grid, final_masses, title): + """Plot final mass as a function of initial mass for a set of pts choices.""" + for label, final_mass in final_masses.items(): + ax.plot(mass_grid, final_mass, marker="o", markersize=4, label=label) + + ax.set_xlabel("Initial Mass [$M_\\odot$]", fontsize=12) + ax.set_ylabel("Final Mass [$M_\\odot$]", fontsize=12) + ax.set_title(title, fontsize=13) + ax.legend(fontsize=10) + ax.grid(True, alpha=0.3) + + +final_masses = final_masses_for_pts(masses, pts_values) +zoom_final_masses = final_masses_for_pts(zoom_masses, zoom_pts_values) + +fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + +plot_final_mass_panel( + axes[0], + masses, + final_masses, + "Effect of pts Timesteps on Final Mass", +) +plot_final_mass_panel( + axes[1], + zoom_masses, + zoom_final_masses, + "Zoomed Final Mass Near Sensitive Region", +) + +fig.suptitle("pts1, pts2, pts3 timestep convergence test", fontsize=14) +plt.tight_layout() +plt.savefig("pts_flagtest.png", dpi=150) +plt.show() diff --git a/docs/settings_examples/plot_windflag.py b/docs/settings_examples/plot_windflag.py new file mode 100644 index 000000000..40e5f9929 --- /dev/null +++ b/docs/settings_examples/plot_windflag.py @@ -0,0 +1,146 @@ +""" +``windflag`` +============ + +This example shows the effect of the ``windflag`` setting on massive-star +evolution. ``windflag`` selects the stellar wind mass-loss prescription used by +COSMIC. The initial metallicity and binary setup are held fixed, while only +``windflag`` is changed. + +The systems are very wide binaries with low-mass companions, so the example is +focused on stellar winds rather than Roche-lobe overflow, common-envelope +evolution, or tides. +""" + +import sys + +sys.path.append("..") +import generate_default_bsedict + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) +plt.style.use("../_static/gallery.mplstyle") + + +windflag_labels = { + -1: "-1 no winds", + 0: "0 SSE/BSE", + 1: "1 StarTrack", + 2: "2 Vink OB/WR", + 3: "3 Vink+LBV", + 5: "5 0.33x flag 3", + 6: "6 Bjorklund+23", + 7: "7 Krticka+24", +} + +windflag_values = list(windflag_labels) +masses = np.arange(5.0, 81.0, 1.0) + +kstar_labels = { + 1: "1: MS", + 2: "2: HG", + 3: "3: GB", + 4: "4: CHeB", + 5: "5: EAGB", + 6: "6: TPAGB", + 7: "7: HeMS", + 8: "8: HeHG", + 9: "9: HeGB", + 10: "10: HeWD", + 11: "11: COWD", + 12: "12: ONeWD", + 13: "13: NS", + 14: "14: BH", + 15: "15: MR", +} + + +def make_wide_binary_grid(): + """Build the wide-binary grid used to isolate stellar wind effects.""" + return InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.ones_like(masses) * 0.1, + porb=np.ones_like(masses) * 1e6, + ecc=np.zeros_like(masses), + tphysf=np.ones_like(masses) * 13700.0, + kstar1=np.ones_like(masses), + kstar2=np.ones_like(masses), + metallicity=np.ones_like(masses) * 0.014, + ) + + +def summarize_final_state(bpp, bcm, windflag): + """Collect final stellar quantities for one wind prescription.""" + final_bpp = bpp.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["windflag"] = windflag + summary["windflag_label"] = windflag_labels[windflag] + summary["initial_mass_1"] = masses[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + summary["final_co_core_1"] = final_bcm["massc_co_layer_1"] + summary["final_kstar_1"] = final_bcm["kstar_1"] + summary["final_evol_type"] = final_bpp["evol_type"] + return summary.reset_index(drop=True) + + +def evolve_windflag_grid(binary_grid): + """Run the same initial grid for each windflag value.""" + summaries = [] + + for windflag in windflag_values: + settings = BSEDict.copy() + settings["windflag"] = windflag + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, + BSEDict=settings, + ) + summaries.append(summarize_final_state(bpp, bcm, windflag)) + + return pd.concat(summaries, ignore_index=True) + + +binary_grid = make_wide_binary_grid() +summaries = evolve_windflag_grid(binary_grid) + +fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharex=True, constrained_layout=True) + +for offset, windflag in zip(np.linspace(-0.14, 0.14, len(windflag_values)), windflag_values): + subset = summaries[summaries["windflag"].eq(windflag)] + label = windflag_labels[windflag] + axes[0].plot(subset["initial_mass_1"], subset["final_mass_1"], label=label) + axes[1].scatter( + subset["initial_mass_1"], + subset["final_kstar_1"] + offset, + s=16, + label=label, + ) + +axes[0].set_xlabel("Initial mass [$M_\\odot$]") +axes[0].set_ylabel("Final mass [$M_\\odot$]") +axes[0].set_title("Final mass") + +axes[1].set_xlabel("Initial mass [$M_\\odot$]") +axes[1].set_ylabel("Final kstar") +axes[1].set_title("Final kstar type") + +used_kstars = sorted(summaries["final_kstar_1"].astype(int).unique()) +axes[1].set_yticks(used_kstars) +axes[1].set_yticklabels([kstar_labels.get(k, str(k)) for k in used_kstars]) + +for ax in axes: + ax.legend(title="windflag", fontsize=7, title_fontsize=8, markerscale=1.2) + ax.grid(True, alpha=0.3) + +fig.suptitle("Effect of windflag at fixed Z=0.014 and zsun=0.014") +plt.savefig("windflag_flagtest.png", dpi=150) +plt.show() diff --git a/docs/settings_examples/plot_zsun.py b/docs/settings_examples/plot_zsun.py new file mode 100644 index 000000000..3c58c8e50 --- /dev/null +++ b/docs/settings_examples/plot_zsun.py @@ -0,0 +1,141 @@ +""" +``zsun`` +======== + +This example shows the effect of the ``zsun`` setting at fixed absolute +metallicity. ``zsun`` sets the reference value COSMIC uses for solar +metallicity. It does not change the stellar metallicity input directly; +instead, it changes the ratio ``Z / zsun`` used by metallicity-dependent wind +prescriptions. + +The stars all have ``Z = 0.014``. Changing ``zsun`` therefore changes how +metal-rich the same stars appear relative to solar, which can change wind mass +loss and the final stellar or remnant type. +""" + +import sys + +sys.path.append("..") +import generate_default_bsedict + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) +plt.style.use("../_static/gallery.mplstyle") + + +zsun_values = [0.010, 0.014, 0.019, 0.020, 0.030] +masses = np.arange(5.0, 81.0, 1.0) + +kstar_labels = { + 1: "1: MS", + 2: "2: HG", + 3: "3: GB", + 4: "4: CHeB", + 5: "5: EAGB", + 6: "6: TPAGB", + 7: "7: HeMS", + 8: "8: HeHG", + 9: "9: HeGB", + 10: "10: HeWD", + 11: "11: COWD", + 12: "12: ONeWD", + 13: "13: NS", + 14: "14: BH", + 15: "15: MR", +} + + +def make_wide_binary_grid(): + """Build a very wide binary grid to isolate wind/remnant effects.""" + n_systems = len(masses) + return InitialBinaryTable.InitialBinaries( + m1=masses, + m2=np.ones(n_systems) * 0.1, + porb=np.ones(n_systems) * 1e6, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * 13700.0, + kstar1=np.ones(n_systems), + kstar2=np.ones(n_systems), + metallicity=np.ones(n_systems) * 0.014, + ) + + +def summarize_final_state(bpp, bcm, zsun): + """Collect final primary quantities for one zsun value.""" + final_bpp = bpp.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["zsun"] = zsun + summary["initial_mass_1"] = masses[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + summary["mass_lost_1"] = summary["initial_mass_1"] - summary["final_mass_1"] + summary["final_kstar_1"] = final_bcm["kstar_1"] + summary["final_evol_type"] = final_bpp["evol_type"] + return summary.reset_index(drop=True) + + +def evolve_zsun_grid(binary_grid): + """Run the same initial grid for each zsun value.""" + summaries = [] + + for zsun in zsun_values: + settings = BSEDict.copy() + settings["zsun"] = zsun + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=binary_grid, + BSEDict=settings, + ) + summaries.append(summarize_final_state(bpp, bcm, zsun)) + + return pd.concat(summaries, ignore_index=True) + + +binary_grid = make_wide_binary_grid() +summaries = evolve_zsun_grid(binary_grid) + +fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharex=True, constrained_layout=True) + +for offset, zsun in zip(np.linspace(-0.10, 0.10, len(zsun_values)), zsun_values): + subset = summaries[summaries["zsun"].eq(zsun)].sort_values("initial_mass_1") + label = f"zsun={zsun:g}" + + axes[0].plot( + subset["initial_mass_1"], + subset["final_mass_1"], + label=label, + ) + axes[1].scatter( + subset["initial_mass_1"], + subset["final_kstar_1"] + offset, + s=14, + label=label, + ) + +axes[0].set_xlabel("Initial mass [$M_\\odot$]") +axes[0].set_ylabel("Final mass [$M_\\odot$]") +axes[0].set_title("Final mass") + +axes[1].set_xlabel("Initial mass [$M_\\odot$]") +axes[1].set_ylabel("Final kstar") +axes[1].set_title("Final kstar type") + +used_kstars = sorted(summaries["final_kstar_1"].astype(int).unique()) +axes[1].set_yticks(used_kstars) +axes[1].set_yticklabels([kstar_labels.get(k, str(k)) for k in used_kstars]) + +for ax in axes: + ax.legend(title="zsun", fontsize=8, title_fontsize=9, markerscale=1.2) + ax.grid(True, alpha=0.3) + +fig.suptitle("Effect of zsun at fixed absolute metallicity Z=0.014") +plt.savefig("zsun_flagtest.png", dpi=150) +plt.show() From 5fb90e2b13328f4ae220b11fc7dc173cb3e7c0a5 Mon Sep 17 00:00:00 2001 From: lukewilner Date: Sat, 11 Jul 2026 16:20:11 -0700 Subject: [PATCH 2/3] Address documentation review comments --- docs/settings_examples/plot_hewind.py | 12 +- docs/settings_examples/plot_neta.py | 12 +- docs/settings_examples/plot_pts.py | 1 - docs/settings_examples/plot_windflag.py | 92 +++++++++++++--- docs/settings_examples/plot_zsun.py | 141 ------------------------ 5 files changed, 96 insertions(+), 162 deletions(-) delete mode 100644 docs/settings_examples/plot_zsun.py diff --git a/docs/settings_examples/plot_hewind.py b/docs/settings_examples/plot_hewind.py index a29aef579..2264fb968 100644 --- a/docs/settings_examples/plot_hewind.py +++ b/docs/settings_examples/plot_hewind.py @@ -26,6 +26,15 @@ BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) plt.style.use("../_static/gallery.mplstyle") +plt.rcParams.update({ + "font.size": 9, + "axes.titlesize": 11, + "axes.labelsize": 10, + "xtick.labelsize": 8, + "ytick.labelsize": 8, + "legend.fontsize": 8, + "legend.title_fontsize": 8, +}) hewind_values = [0.0, 0.25, 0.5, 0.75, 1.0] @@ -93,7 +102,7 @@ def evolve_hewind_grid(helium_grid, windflag=0, dtp=2.0): summaries = evolve_hewind_grid(helium_grid, windflag=0, dtp=2.0) baseline = summaries[summaries["hewind"].eq(0.5)].set_index("initial_he_mass_1") -fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), sharex=True, constrained_layout=True) +fig, axes = plt.subplots(1, 2, figsize=(10, 4.0), sharex=True, constrained_layout=True) for hewind in hewind_values: subset = summaries[summaries["hewind"].eq(hewind)].sort_values( @@ -128,6 +137,5 @@ def evolve_hewind_grid(helium_grid, windflag=0, dtp=2.0): ax.legend(title="hewind", fontsize=8, title_fontsize=9, markerscale=1.2) ax.grid(True, alpha=0.3) -fig.suptitle("Effect of hewind on naked helium stars, windflag=0") plt.savefig("hewind_flagtest.png", dpi=150) plt.show() diff --git a/docs/settings_examples/plot_neta.py b/docs/settings_examples/plot_neta.py index 007e442ce..5a9ab96a3 100644 --- a/docs/settings_examples/plot_neta.py +++ b/docs/settings_examples/plot_neta.py @@ -27,6 +27,15 @@ BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) plt.style.use("../_static/gallery.mplstyle") +plt.rcParams.update({ + "font.size": 9, + "axes.titlesize": 11, + "axes.labelsize": 10, + "xtick.labelsize": 8, + "ytick.labelsize": 8, + "legend.fontsize": 8, + "legend.title_fontsize": 8, +}) neta_values = [0.01, 0.25, 0.5, 1.0, 2.0] @@ -115,7 +124,7 @@ def plot_final_mass(ax, data, title): & (summaries["initial_mass_1"] <= 1.5) ] -fig, axes = plt.subplots(1, 2, figsize=(13, 4.8)) +fig, axes = plt.subplots(1, 2, figsize=(10, 4.0)) plot_final_mass( axes[0], @@ -128,7 +137,6 @@ def plot_final_mass(ax, data, title): "Full mass grid", ) -fig.suptitle("Effect of Reimers wind coefficient neta") plt.tight_layout() plt.savefig("neta_flagtest.png", dpi=150) plt.show() diff --git a/docs/settings_examples/plot_pts.py b/docs/settings_examples/plot_pts.py index a682ebfbc..b5d67b0fe 100644 --- a/docs/settings_examples/plot_pts.py +++ b/docs/settings_examples/plot_pts.py @@ -134,7 +134,6 @@ def plot_final_mass_panel(ax, mass_grid, final_masses, title): "Zoomed Final Mass Near Sensitive Region", ) -fig.suptitle("pts1, pts2, pts3 timestep convergence test", fontsize=14) plt.tight_layout() plt.savefig("pts_flagtest.png", dpi=150) plt.show() diff --git a/docs/settings_examples/plot_windflag.py b/docs/settings_examples/plot_windflag.py index 40e5f9929..fb4f95071 100644 --- a/docs/settings_examples/plot_windflag.py +++ b/docs/settings_examples/plot_windflag.py @@ -9,7 +9,8 @@ The systems are very wide binaries with low-mass companions, so the example is focused on stellar winds rather than Roche-lobe overflow, common-envelope -evolution, or tides. +evolution, or tides. The inset in the final-mass panel shows the same +comparison on log-log axes. """ import sys @@ -18,6 +19,7 @@ import generate_default_bsedict import matplotlib.pyplot as plt +from matplotlib.ticker import LogFormatterMathtext, LogLocator, NullFormatter import numpy as np import pandas as pd @@ -26,17 +28,27 @@ BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) plt.style.use("../_static/gallery.mplstyle") +plt.rcParams.update({ + "font.size": 9, + "axes.titlesize": 12, + "axes.labelsize": 10, + "xtick.labelsize": 8, + "ytick.labelsize": 8, + "legend.fontsize": 7, + "legend.title_fontsize": 8, + "figure.titlesize": 12, +}) windflag_labels = { - -1: "-1 no winds", + -1: "-1 none", 0: "0 SSE/BSE", 1: "1 StarTrack", - 2: "2 Vink OB/WR", + 2: "2 Vink", 3: "3 Vink+LBV", - 5: "5 0.33x flag 3", - 6: "6 Bjorklund+23", - 7: "7 Krticka+24", + 5: "5 0.33x", + 6: "6 Bjorklund", + 7: "7 Krticka", } windflag_values = list(windflag_labels) @@ -112,24 +124,64 @@ def evolve_windflag_grid(binary_grid): binary_grid = make_wide_binary_grid() summaries = evolve_windflag_grid(binary_grid) -fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharex=True, constrained_layout=True) +fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.2), constrained_layout=False) +fig.subplots_adjust(left=0.08, right=0.98, bottom=0.28, top=0.82, wspace=0.38) +log_inset = axes[0].inset_axes([0.09, 0.55, 0.42, 0.36]) for offset, windflag in zip(np.linspace(-0.14, 0.14, len(windflag_values)), windflag_values): subset = summaries[summaries["windflag"].eq(windflag)] label = windflag_labels[windflag] - axes[0].plot(subset["initial_mass_1"], subset["final_mass_1"], label=label) + log_subset = subset[ + subset["initial_mass_1"].gt(0) & subset["final_mass_1"].gt(0) + ] + axes[0].plot( + log_subset["initial_mass_1"], + log_subset["final_mass_1"], + label=label, + ) + log_inset.plot( + log_subset["initial_mass_1"], + log_subset["final_mass_1"], + lw=1.0, + ) axes[1].scatter( subset["initial_mass_1"], subset["final_kstar_1"] + offset, - s=16, + s=12, label=label, ) -axes[0].set_xlabel("Initial mass [$M_\\odot$]") -axes[0].set_ylabel("Final mass [$M_\\odot$]") +axes[0].set_xlabel("Initial Mass [$M_\\odot$]") +axes[0].set_ylabel("Final Mass [$M_\\odot$]") axes[0].set_title("Final mass") - -axes[1].set_xlabel("Initial mass [$M_\\odot$]") +axes[0].set_xscale("linear") +axes[0].set_yscale("linear") +axes[0].set_xticks([10, 20, 40, 60, 80]) +axes[0].set_yticks([0, 15, 30, 45, 60]) + +log_positive = summaries[ + summaries["initial_mass_1"].gt(0) & summaries["final_mass_1"].gt(0) +] +log_inset.set_xscale("log", base=10) +log_inset.set_yscale("log", base=10) +log_inset.set_xlim(masses.min(), masses.max()) +log_inset.set_ylim( + log_positive["final_mass_1"].min() * 0.8, + log_positive["final_mass_1"].max() * 1.2, +) +log_inset.xaxis.set_major_locator(LogLocator(base=10, numticks=3)) +log_inset.yaxis.set_major_locator(LogLocator(base=10, numticks=3)) +log_inset.xaxis.set_major_formatter(LogFormatterMathtext(base=10)) +log_inset.yaxis.set_major_formatter(LogFormatterMathtext(base=10)) +log_inset.xaxis.set_minor_locator(LogLocator(base=10, subs=np.arange(2, 10) * 0.1)) +log_inset.yaxis.set_minor_locator(LogLocator(base=10, subs=np.arange(2, 10) * 0.1)) +log_inset.xaxis.set_minor_formatter(NullFormatter()) +log_inset.yaxis.set_minor_formatter(NullFormatter()) +log_inset.tick_params(axis="both", which="major", labelsize=8, pad=1) +log_inset.tick_params(axis="both", which="minor", labelsize=0) +log_inset.grid(True, which="both", alpha=0.25) + +axes[1].set_xlabel("Initial Mass [$M_\\odot$]") axes[1].set_ylabel("Final kstar") axes[1].set_title("Final kstar type") @@ -138,9 +190,17 @@ def evolve_windflag_grid(binary_grid): axes[1].set_yticklabels([kstar_labels.get(k, str(k)) for k in used_kstars]) for ax in axes: - ax.legend(title="windflag", fontsize=7, title_fontsize=8, markerscale=1.2) ax.grid(True, alpha=0.3) -fig.suptitle("Effect of windflag at fixed Z=0.014 and zsun=0.014") -plt.savefig("windflag_flagtest.png", dpi=150) +handles, labels = axes[0].get_legend_handles_labels() +fig.legend( + handles, + labels, + title="windflag", + loc="lower center", + ncol=4, + fontsize=7, + title_fontsize=8, + bbox_to_anchor=(0.5, 0.03), +) plt.show() diff --git a/docs/settings_examples/plot_zsun.py b/docs/settings_examples/plot_zsun.py deleted file mode 100644 index 3c58c8e50..000000000 --- a/docs/settings_examples/plot_zsun.py +++ /dev/null @@ -1,141 +0,0 @@ -""" -``zsun`` -======== - -This example shows the effect of the ``zsun`` setting at fixed absolute -metallicity. ``zsun`` sets the reference value COSMIC uses for solar -metallicity. It does not change the stellar metallicity input directly; -instead, it changes the ratio ``Z / zsun`` used by metallicity-dependent wind -prescriptions. - -The stars all have ``Z = 0.014``. Changing ``zsun`` therefore changes how -metal-rich the same stars appear relative to solar, which can change wind mass -loss and the final stellar or remnant type. -""" - -import sys - -sys.path.append("..") -import generate_default_bsedict - -import matplotlib.pyplot as plt -import numpy as np -import pandas as pd - -from cosmic.evolve import Evolve -from cosmic.sample import InitialBinaryTable - -BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) -plt.style.use("../_static/gallery.mplstyle") - - -zsun_values = [0.010, 0.014, 0.019, 0.020, 0.030] -masses = np.arange(5.0, 81.0, 1.0) - -kstar_labels = { - 1: "1: MS", - 2: "2: HG", - 3: "3: GB", - 4: "4: CHeB", - 5: "5: EAGB", - 6: "6: TPAGB", - 7: "7: HeMS", - 8: "8: HeHG", - 9: "9: HeGB", - 10: "10: HeWD", - 11: "11: COWD", - 12: "12: ONeWD", - 13: "13: NS", - 14: "14: BH", - 15: "15: MR", -} - - -def make_wide_binary_grid(): - """Build a very wide binary grid to isolate wind/remnant effects.""" - n_systems = len(masses) - return InitialBinaryTable.InitialBinaries( - m1=masses, - m2=np.ones(n_systems) * 0.1, - porb=np.ones(n_systems) * 1e6, - ecc=np.zeros(n_systems), - tphysf=np.ones(n_systems) * 13700.0, - kstar1=np.ones(n_systems), - kstar2=np.ones(n_systems), - metallicity=np.ones(n_systems) * 0.014, - ) - - -def summarize_final_state(bpp, bcm, zsun): - """Collect final primary quantities for one zsun value.""" - final_bpp = bpp.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() - final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() - - summary = pd.DataFrame(index=final_bcm.index) - summary["zsun"] = zsun - summary["initial_mass_1"] = masses[summary.index.astype(int)] - summary["final_mass_1"] = final_bcm["mass_1"] - summary["mass_lost_1"] = summary["initial_mass_1"] - summary["final_mass_1"] - summary["final_kstar_1"] = final_bcm["kstar_1"] - summary["final_evol_type"] = final_bpp["evol_type"] - return summary.reset_index(drop=True) - - -def evolve_zsun_grid(binary_grid): - """Run the same initial grid for each zsun value.""" - summaries = [] - - for zsun in zsun_values: - settings = BSEDict.copy() - settings["zsun"] = zsun - settings["random_seed"] = 1 - - bpp, bcm, initC, kick_info = Evolve.evolve( - initialbinarytable=binary_grid, - BSEDict=settings, - ) - summaries.append(summarize_final_state(bpp, bcm, zsun)) - - return pd.concat(summaries, ignore_index=True) - - -binary_grid = make_wide_binary_grid() -summaries = evolve_zsun_grid(binary_grid) - -fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharex=True, constrained_layout=True) - -for offset, zsun in zip(np.linspace(-0.10, 0.10, len(zsun_values)), zsun_values): - subset = summaries[summaries["zsun"].eq(zsun)].sort_values("initial_mass_1") - label = f"zsun={zsun:g}" - - axes[0].plot( - subset["initial_mass_1"], - subset["final_mass_1"], - label=label, - ) - axes[1].scatter( - subset["initial_mass_1"], - subset["final_kstar_1"] + offset, - s=14, - label=label, - ) - -axes[0].set_xlabel("Initial mass [$M_\\odot$]") -axes[0].set_ylabel("Final mass [$M_\\odot$]") -axes[0].set_title("Final mass") - -axes[1].set_xlabel("Initial mass [$M_\\odot$]") -axes[1].set_ylabel("Final kstar") -axes[1].set_title("Final kstar type") - -used_kstars = sorted(summaries["final_kstar_1"].astype(int).unique()) -axes[1].set_yticks(used_kstars) -axes[1].set_yticklabels([kstar_labels.get(k, str(k)) for k in used_kstars]) - -for ax in axes: - ax.legend(title="zsun", fontsize=8, title_fontsize=9, markerscale=1.2) - ax.grid(True, alpha=0.3) - -fig.suptitle("Effect of zsun at fixed absolute metallicity Z=0.014") -plt.savefig("zsun_flagtest.png", dpi=150) -plt.show() From 8ae256dd8f22a2977d3dab5dac4cd0fd2ef3012d Mon Sep 17 00:00:00 2001 From: lukewilner Date: Tue, 21 Jul 2026 15:59:21 -0700 Subject: [PATCH 3/3] Add CE, LBV, alpha1, and lambdaf setting examples --- docs/settings_examples/plot_LBV_flag.py | 113 +++++- docs/settings_examples/plot_alpha1_lambdaf.py | 340 +++++++++++++++++ docs/settings_examples/plot_ce_flags.py | 348 ++++++++++++++++++ 3 files changed, 787 insertions(+), 14 deletions(-) create mode 100644 docs/settings_examples/plot_alpha1_lambdaf.py create mode 100644 docs/settings_examples/plot_ce_flags.py diff --git a/docs/settings_examples/plot_LBV_flag.py b/docs/settings_examples/plot_LBV_flag.py index 8b37b8411..016c58f9d 100644 --- a/docs/settings_examples/plot_LBV_flag.py +++ b/docs/settings_examples/plot_LBV_flag.py @@ -4,7 +4,12 @@ This example shows the effect of the ``LBV_flag`` on the evolution of massive stars. -The ``LBV_flag`` controls which prescription to use for LBV-like mass loss for stars that are above the Humphreys-Davidson limit. The plots below show the evolution of a few single stars and highlight the LBV regime in which we'd expect to see differences with this flag. +The ``LBV_flag`` controls which prescription to use for LBV-like mass loss for +stars that are above the Humphreys-Davidson limit. The first figures show the +evolution of a few single stars on the HR diagram and highlight the LBV regime +where we'd expect to see differences with this flag. The final figure compares +the final primary mass over a wider initial-mass grid, showing how the LBV wind +prescription changes the amount of mass retained by massive single stars. """ #---------------------------------------------------------------------------------- @@ -22,6 +27,7 @@ #---------------------------------------------------------------------------------- import numpy as np +import pandas as pd import astropy.units as u import astropy.constants as consts @@ -35,25 +41,26 @@ def LBV_limit(T_eff): # create a small grid of single stars -masses = [30, 33, 36, 40, 45, 50] +hr_masses = [30, 33, 36, 40, 45, 50] grid = InitialBinaryTable.InitialBinaries( - m1=masses, - m2=np.zeros(len(masses)), - porb=np.ones(len(masses))*-1.0, - ecc=np.ones(len(masses))*0.0, - tphysf=np.ones(len(masses))*13700.0, - kstar1=np.ones(len(masses)), - kstar2=np.ones(len(masses)), - metallicity=np.ones(len(masses))*0.014*0.01 + m1=hr_masses, + m2=np.zeros(len(hr_masses)), + porb=np.ones(len(hr_masses))*-1.0, + ecc=np.ones(len(hr_masses))*0.0, + tphysf=np.ones(len(hr_masses))*13700.0, + kstar1=np.ones(len(hr_masses)), + kstar2=np.ones(len(hr_masses)), + metallicity=np.ones(len(hr_masses))*0.014*0.01 ) # loop over different flag choices for LBV_flag, label in zip([0, 1, 2], ["No LBV mass loss", "LBV mass loss following Hurley+2000", "LBV mass loss following Belczynski+2008"]): # evolve with updated BSEDict - BSEDict["LBV_flag"] = LBV_flag + settings = BSEDict.copy() + settings["LBV_flag"] = LBV_flag bpp, bcm, initC, kick_info = Evolve.evolve( - initialbinarytable=grid, BSEDict=BSEDict, + initialbinarytable=grid, BSEDict=settings, dtp=0 ) @@ -72,7 +79,7 @@ def LBV_limit(T_eff): ax.plot(log_teff[mask], log_lum[mask], color='grey', alpha=0.5, lw=0.5, zorder=-1) # annotate the first point with the initial mass - ax.annotate(f'{bcm[bcm["bin_num"] == bin_num]["mass_1"].values[0]:.0f} M$_\odot$', (log_teff[0], log_lum[0]), + ax.annotate(f'{bcm[bcm["bin_num"] == bin_num]["mass_1"].values[0]:.0f} M$_\\odot$', (log_teff[0], log_lum[0]), fontsize=14, ha="right", color='black', va='top') # annotate the various kstar types @@ -101,4 +108,82 @@ def LBV_limit(T_eff): ax.invert_xaxis() - plt.show() \ No newline at end of file + plt.show() + + +# compare the final masses over a wider initial-mass grid +LBV_flag_labels = { + 0: "0: no LBV winds", + 1: "1: Hurley+2000", + 2: "2: Belczynski+2008", +} +LBV_flag_values = list(LBV_flag_labels) +mass_grid = np.arange(25.0, 81.0, 1.0) + + +def make_final_mass_grid(): + """Build the massive single-star grid used for the final-mass comparison.""" + n_systems = len(mass_grid) + return InitialBinaryTable.InitialBinaries( + m1=mass_grid, + m2=np.zeros(n_systems), + porb=np.ones(n_systems) * -1.0, + ecc=np.zeros(n_systems), + tphysf=np.ones(n_systems) * 13700.0, + kstar1=np.ones(n_systems), + kstar2=np.zeros(n_systems), + metallicity=np.ones(n_systems) * 0.014, + ) + + +def summarize_final_masses(bcm, LBV_flag): + """Collect final primary masses for one LBV_flag value.""" + final_bcm = bcm.sort_values(["bin_num", "tphys"]).groupby("bin_num").last() + + summary = pd.DataFrame(index=final_bcm.index) + summary["LBV_flag"] = LBV_flag + summary["initial_mass_1"] = mass_grid[summary.index.astype(int)] + summary["final_mass_1"] = final_bcm["mass_1"] + return summary.reset_index(drop=True) + + +def evolve_final_mass_grid(final_mass_grid): + """Run the same single-star grid for each LBV_flag value.""" + summaries = [] + + for LBV_flag in LBV_flag_values: + settings = BSEDict.copy() + settings["LBV_flag"] = LBV_flag + settings["random_seed"] = 1 + + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=final_mass_grid, + BSEDict=settings, + dtp=0, + ) + summaries.append(summarize_final_masses(bcm, LBV_flag)) + + return pd.concat(summaries, ignore_index=True) + + +final_mass_grid = make_final_mass_grid() +final_mass_summaries = evolve_final_mass_grid(final_mass_grid) + +fig, ax = plt.subplots(figsize=(8.2, 6.0)) + +for LBV_flag in LBV_flag_values: + subset = final_mass_summaries[final_mass_summaries["LBV_flag"].eq(LBV_flag)] + ax.plot( + subset["initial_mass_1"], + subset["final_mass_1"], + label=LBV_flag_labels[LBV_flag], + ) + +ax.set_xlabel("Initial primary mass [$M_\\odot$]") +ax.set_ylabel("Final primary mass [$M_\\odot$]", fontsize=14) +ax.set_title("Effect of LBV_flag on final primary mass", fontsize=15) +ax.legend(fontsize=8) +ax.grid(True, alpha=0.3) + +plt.tight_layout() +plt.show() diff --git a/docs/settings_examples/plot_alpha1_lambdaf.py b/docs/settings_examples/plot_alpha1_lambdaf.py new file mode 100644 index 000000000..8a5e976ed --- /dev/null +++ b/docs/settings_examples/plot_alpha1_lambdaf.py @@ -0,0 +1,340 @@ +""" +alpha1 and lambdaf +================== + +This example compares two common-envelope parameters using real COSMIC +separation histories. ``alpha1`` changes the efficiency with which orbital +energy ejects the common envelope. ``lambdaf`` changes the envelope +binding-energy lambda prescription; the negative values used here select fixed +lambda values. + +Each panel holds the initial binary fixed and changes only the parameter shown +in that panel. Open circles mark surviving binaries, x markers mark mergers, +and triangles mark disruptions at the final real COSMIC output point shown. +Tracks that overlap exactly are given tiny horizontal plotting offsets for +visibility only; the COSMIC output times and separations are unchanged. +""" + +import copy +import sys +import warnings + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +sys.path.append("..") +import generate_default_bsedict + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) + + +EXAMPLES = [ + { + "title": "alpha1: CE efficiency", + "flag_key": "alpha1", + "values": [0.2, 1.0, 5.0], + "labels": { + 0.2: r"$\alpha$ = 0.2", + 1.0: r"$\alpha$ = 1.0", + 5.0: r"$\alpha$ = 5.0", + }, + "initial": { + "m1": 8.0, + "m2": 2.8, + "porb": 1584.893192461114, + "tphysf": 70.0, + }, + "plot_end_pad": 0.35, + }, + { + "title": "lambdaf: fixed envelope lambda", + "flag_key": "lambdaf", + "values": [-0.1, -0.5, -1.0], + "labels": { + -0.1: r"fixed $\lambda$ = 0.1", + -0.5: r"fixed $\lambda$ = 0.5", + -1.0: r"fixed $\lambda$ = 1.0", + }, + "initial": { + "m1": 13.0, + "m2": 3.25, + "porb": 2238.72113856834, + "tphysf": 50.0, + }, + "plot_end_pad": 0.5, + }, +] + +LINE_WIDTHS = { + "alpha1": { + 0.2: 3.0, + 1.0: 2.2, + 5.0: 1.4, + }, + "lambdaf": { + -0.1: 3.0, + -0.5: 2.2, + -1.0: 1.4, + }, +} + +LINE_ZORDERS = { + "alpha1": { + 0.2: 3, + 1.0: 4, + 5.0: 5, + }, + "lambdaf": { + -0.1: 3, + -0.5: 4, + -1.0: 5, + }, +} + +DISPLAY_X_OFFSETS = { + "alpha1": { + 0.2: 0.00, + 1.0: 0.15, + 5.0: 0.30, + }, + "lambdaf": { + -0.1: 0.00, + -0.5: 0.15, + -1.0: 0.30, + }, +} + + +def make_binary(initial): + """Build a one-row InitialBinaryTable.""" + return InitialBinaryTable.InitialBinaries( + m1=np.array([initial["m1"]]), + m2=np.array([initial["m2"]]), + porb=np.array([initial["porb"]]), + ecc=np.zeros(1), + tphysf=np.ones(1) * initial["tphysf"], + kstar1=np.ones(1), + kstar2=np.ones(1), + metallicity=np.ones(1) * 0.014, + ) + + +def set_parameter(settings, key, value): + """Set one CE parameter while preserving COSMIC's expected type.""" + if key == "alpha1": + settings[key] = [value, value] + else: + settings[key] = value + + +def evolve_one(example, value): + """Evolve one binary while changing only the selected CE parameter.""" + settings = copy.deepcopy(BSEDict) + set_parameter(settings, example["flag_key"], value) + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="At least one of your initial binaries is starting in Roche Lobe Overflow:*", + category=UserWarning, + ) + bpp, bcm, initC, kick_info = Evolve.evolve( + initialbinarytable=make_binary(example["initial"]), + BSEDict=settings, + dtp=0.5, + randomseed=1, + ) + + return bpp, bcm + + +def first_ce_time(bpp): + """Return the first confirmed common-envelope event time.""" + ce_rows = bpp[bpp["evol_type"].eq(7)] + if ce_rows.empty: + raise RuntimeError("Selected binary does not enter common-envelope evolution.") + return ce_rows["tphys"].iloc[0] + + +def terminal_time(bpp): + """Return the first terminal-event time, or NaN if none is present.""" + terminal_rows = bpp[bpp["evol_type"].eq(6)] + if terminal_rows.empty: + return np.nan + return terminal_rows["tphys"].iloc[0] + + +def final_status(bpp, bcm): + """Read the final outcome from COSMIC outputs.""" + if np.isfinite(terminal_time(bpp)): + return "merges" + + final_row = bcm.sort_values("tphys").iloc[-1] + if final_row["bin_state"] == 0 and final_row["sep"] > 0: + return "survives" + if final_row["bin_state"] == 1: + return "merges" + return "disrupts" + + +def separation_track(bpp, bcm): + """Combine sampled bcm rows with bpp event rows into a real separation track.""" + sampled_rows = bcm[["tphys", "sep"]].copy() + sampled_rows["source_order"] = 0 + + event_rows = bpp[["tphys", "sep"]].copy() + event_rows["source_order"] = np.arange(len(event_rows)) + 1 + + track = pd.concat([sampled_rows, event_rows], ignore_index=True) + track = track[np.isfinite(track["tphys"]) & np.isfinite(track["sep"])] + track = track[track["sep"].gt(0)] + track = track.drop_duplicates(subset=["tphys", "sep"], keep="first") + track = track.sort_values(["tphys", "source_order"]).reset_index(drop=True) + + stop_time = terminal_time(bpp) + if np.isfinite(stop_time): + track = track[track["tphys"].le(stop_time)].copy() + + assert len(track) > 0 + assert np.all(np.isfinite(track["tphys"])) + assert np.all(np.isfinite(track["sep"])) + assert np.all(track["sep"] > 0) + assert np.all(np.diff(track["tphys"]) >= 0) + + return track + + +def mark_outcome(ax, track, status, color, x_offset=0.0, marker_scale=1.0, zorder=12): + """Place the endpoint marker at the final real plotted point.""" + final = track.iloc[-1] + marker_time = final["tphys"] + x_offset + if status == "survives": + ax.scatter( + marker_time, + final["sep"], + marker="o", + s=52 * marker_scale, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + elif status == "merges": + ax.scatter( + marker_time, + final["sep"], + marker="x", + s=62 * marker_scale, + color=color, + linewidths=2.1, + zorder=zorder, + ) + else: + ax.scatter( + marker_time, + final["sep"], + marker="^", + s=62 * marker_scale, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + + +fig, axes = plt.subplots( + 1, + 2, + figsize=(12, 5), + constrained_layout=True, +) + +for ax, example in zip(axes, EXAMPLES): + colors = plt.get_cmap("tab10").colors[: len(example["values"])] + records = [] + + for value, color in zip(example["values"], colors): + bpp, bcm = evolve_one(example, value) + ce_time = first_ce_time(bpp) + track = separation_track(bpp, bcm) + status = final_status(bpp, bcm) + records.append( + { + "value": value, + "label": example["labels"][value], + "color": color, + "track": track, + "status": status, + "ce_time": ce_time, + } + ) + + ce_times = np.array([record["ce_time"] for record in records]) + if not np.allclose(ce_times, ce_times[0]): + raise RuntimeError("Selected runs do not share the same CE onset time.") + ce_time = ce_times[0] + + for record in records: + linewidth = LINE_WIDTHS[example["flag_key"]][record["value"]] + zorder = LINE_ZORDERS[example["flag_key"]][record["value"]] + display_offset = DISPLAY_X_OFFSETS[example["flag_key"]][record["value"]] + marker_offset = display_offset + marker_scale = 1.0 + marker_zorder = 12 + zorder + if example["flag_key"] == "lambdaf" and record["status"] == "merges": + marker_offset += {-0.1: -0.06, -0.5: 0.06}.get(record["value"], 0.0) + marker_scale = {-0.1: 1.15, -0.5: 0.95}.get(record["value"], 1.0) + ax.plot( + record["track"]["tphys"] + display_offset, + record["track"]["sep"], + color=record["color"], + lw=linewidth, + linestyle="-", + marker=None, + label=record["label"], + zorder=zorder, + ) + mark_outcome( + ax, + record["track"], + record["status"], + record["color"], + x_offset=marker_offset, + marker_scale=marker_scale, + zorder=marker_zorder, + ) + + combined = pd.concat([record["track"] for record in records], ignore_index=True) + max_time = max(record["track"]["tphys"].max() for record in records) + ax.set_xlim(0, max_time + example["plot_end_pad"]) + ax.set_ylim(combined["sep"].min() * 0.85, combined["sep"].max() * 1.35) + ax.set_yscale("log") + ax.axvline(ce_time, color="0.45", lw=0.9, zorder=0) + x0, x1 = ax.get_xlim() + ax.text( + ce_time + 0.028 * (x1 - x0), + 0.93, + "CE onset", + transform=ax.get_xaxis_transform(), + ha="left", + va="top", + fontsize=7, + color="0.25", + ) + ax.set_title(example["title"], fontsize=11) + ax.set_xlabel("Time [Myr]") + ax.set_ylabel(r"Orbital separation [R$_\odot$]") + ax.grid(True, which="both", alpha=0.25) + ax.legend( + loc="best", + fontsize=8, + frameon=True, + title=None, + handlelength=1.8, + ) + +plt.show() diff --git a/docs/settings_examples/plot_ce_flags.py b/docs/settings_examples/plot_ce_flags.py new file mode 100644 index 000000000..d729d1852 --- /dev/null +++ b/docs/settings_examples/plot_ce_flags.py @@ -0,0 +1,348 @@ +""" +ceflag and cemergeflag +====================== + +This example compares two common-envelope-related COSMIC flags using one +representative binary per flag. Each panel evolves the same initial binary +multiple times while changing only the flag shown in that panel, then plots +orbital period as a function of time. + +Single-binary examples are useful for these flags because the important +behavior often happens at a boundary between common-envelope survival and +merger. + +Each panel varies one COSMIC flag while holding the initial binary fixed. +Tracks show orbital-period evolution; endpoint markers indicate whether the +binary ultimately survives or terminates through merger/disruption. The time +axis is limited to the interaction window where the flag differences are most +visible. +""" + +import copy +import sys +import warnings + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from matplotlib.lines import Line2D + +from cosmic.evolve import Evolve +from cosmic.sample import InitialBinaryTable + +sys.path.append("..") +import generate_default_bsedict + +BSEDict = generate_default_bsedict.get_default_BSE_settings(to_python=True) + + +EXAMPLES = [ + { + "title": "ceflag: CE orbital energy", + "flag_key": "ceflag", + "flags": { + 0: "0: core", + 1: "1: total (default)", + }, + "initial": { + "m1": 32.0, + "m2": 8.0, + "porb": 200.0, + "kstar1": 1, + "kstar2": 1, + "tphysf": 300.0, + }, + "annotation": "same binary, different CE energy prescription", + "plot_end": 7.0, + }, + { + "title": "cemergeflag: CE merger treatment", + "flag_key": "cemergeflag", + "flags": { + 0: "0: alternate", + 1: "1: default", + }, + "initial": { + "m1": 8.0, + "m2": 8.0, + "porb": 1500.0, + "kstar1": 1, + "kstar2": 1, + "tphysf": 500.0, + }, + "annotation": "0 permits CE evolution; 1 forces merger for this system", + "plot_end": 45.0, + }, +] + + +def make_binary(initial): + """Build one InitialBinaryTable row from an example dictionary.""" + return InitialBinaryTable.InitialBinaries( + m1=np.array([initial["m1"]]), + m2=np.array([initial["m2"]]), + porb=np.array([initial["porb"]]), + ecc=np.zeros(1), + tphysf=np.ones(1) * initial["tphysf"], + kstar1=np.ones(1) * initial["kstar1"], + kstar2=np.ones(1) * initial["kstar2"], + metallicity=np.ones(1) * 0.014, + ) + + +def evolve_example(example, flag_value): + """Evolve one example binary for one flag value.""" + settings = copy.deepcopy(BSEDict) + settings[example["flag_key"]] = flag_value + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", + message="At least one of your initial binaries is starting in Roche Lobe Overflow:*", + category=UserWarning, + ) + bpp, bcm, initC, _ = Evolve.evolve( + initialbinarytable=make_binary(example["initial"]), + BSEDict=settings, + dtp=0.2, + randomseed=1, + ) + + return bpp, bcm + + +def orbital_period_track(bpp, bcm): + """Combine sampled evolution rows with event rows into one real track.""" + sampled_rows = bcm[["tphys", "porb", "bin_state"]].copy() + sampled_rows["source_order"] = 0 + + event_rows = bpp[["tphys", "porb"]].copy() + event_rows["bin_state"] = np.nan + event_rows["source_order"] = np.arange(len(event_rows)) + 1 + + track = pd.concat([sampled_rows, event_rows], ignore_index=True) + track = track[track["porb"].gt(0)].drop_duplicates( + subset=["tphys", "porb", "bin_state"], keep="first" + ) + track = track.sort_values(["tphys", "source_order"]).reset_index(drop=True) + + time = track["tphys"].to_numpy() + period = track["porb"].to_numpy() + if not np.all(np.diff(time) >= 0): + raise RuntimeError("Orbital-period track is not time sorted.") + if not np.all(np.isfinite(period)) or not np.all(period > 0): + raise RuntimeError("Orbital-period track contains invalid period values.") + + return track + + +def final_status(bcm): + """Return a short final-status string.""" + final_row = bcm.sort_values("tphys").iloc[-1] + if final_row["bin_state"] == 0 and final_row["porb"] > 0: + return "survives" + if final_row["bin_state"] == 1: + return "merges" + return "disrupts" + + +def plotted_track(track, plot_end): + """Return only real COSMIC rows inside the panel time window.""" + if track.empty: + return track + + plot_track = track[track["tphys"].le(plot_end)].copy() + if plot_track.empty: + plot_track = track.iloc[[0]].copy() + + return plot_track + + +def mark_final_outcome( + ax, + track, + status, + color, + marker_size=70, + zorder=5, +): + """Mark the plotted endpoint with a consistent simple symbol.""" + if track.empty: + return + + final_point = track.iloc[-1] + if status == "survives": + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="o", + s=marker_size * 0.75, + facecolor="white", + edgecolor=color, + linewidths=1.8, + zorder=zorder, + ) + elif status == "merges": + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="x", + s=marker_size, + linewidths=2.2, + color=color, + zorder=zorder, + ) + else: + ax.scatter( + final_point["tphys"], + final_point["porb"], + marker="^", + s=marker_size, + facecolor="white", + edgecolor=color, + linewidths=2.0, + zorder=zorder, + ) + + +def first_event_time(bpp, event_type): + """Return the first time for an event type, or NaN if absent.""" + rows = bpp[bpp["evol_type"].eq(event_type)] + if rows.empty: + return np.nan + return rows["tphys"].iloc[0] + + +def first_rlo_time(bpp): + """Return the first recorded RLO time, or NaN if absent.""" + rrlo = bpp[["RRLO_1", "RRLO_2"]].max(axis=1) + rows = bpp[rrlo.ge(1.0)] + if rows.empty: + return np.nan + return rows["tphys"].iloc[0] + + +def primary_event(bpp): + """Return the single event label/time most relevant to the panel.""" + ce_time = first_event_time(bpp, 7) + if np.isfinite(ce_time): + return "CE/RLO", ce_time + + rlo_time = first_rlo_time(bpp) + if np.isfinite(rlo_time): + return "CE/RLO", rlo_time + + return "", np.nan + + +def annotate_event(ax, label, event_time): + """Draw one readable event marker.""" + if not label or not np.isfinite(event_time): + return + + ax.axvline(event_time, color="0.35", ls=":", lw=1.4, zorder=4) + + +def line_width(flag_value): + """Use a consistent width for the two CE-flag tracks.""" + return 2.0 + + +def endpoint_marker_size(flag_value): + """Use a consistent endpoint marker size.""" + return 70 + + +fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), constrained_layout=True) +axes = axes.ravel() + +for ax, example in zip(axes, EXAMPLES): + records = [] + reference_record = None + + for flag_value, label in example["flags"].items(): + bpp, bcm = evolve_example(example, flag_value) + track = orbital_period_track(bpp, bcm) + status = final_status(bcm) + record = { + "flag_value": flag_value, + "label": label, + "bpp": bpp, + "bcm": bcm, + "track": track, + "status": status, + } + records.append(record) + + if flag_value == max(example["flags"]): + reference_record = record + + event_label, event_time = primary_event(reference_record["bpp"]) + all_plotted_tracks = [] + plot_end = example["plot_end"] + for plot_index, record in enumerate(records): + track = plotted_track(record["track"], plot_end) + all_plotted_tracks.append(track) + linewidth = line_width(record["flag_value"]) + zorder = 3 + record["flag_value"] + line, = ax.plot( + track["tphys"], + track["porb"], + lw=linewidth, + ls="-", + label=record["label"], + zorder=zorder, + ) + record["color"] = line.get_color() + mark_final_outcome( + ax, + track, + record["status"], + line.get_color(), + marker_size=endpoint_marker_size(record["flag_value"]), + zorder=5, + ) + + annotate_event(ax, event_label, event_time) + combined = pd.concat(all_plotted_tracks, ignore_index=True) + ax.set_yscale("log") + if not combined.empty: + max_time = combined["tphys"].max() + ax.set_xlim(0, max_time * 1.08) + min_porb = combined["porb"].min() + max_porb = combined["porb"].max() + ax.set_ylim(min_porb * 0.65, max_porb * 1.6) + ax.set_title(example["title"], fontsize=11) + ax.set_xlabel("Time [Myr]") + ax.set_ylabel("Orbital period [days]") + ax.grid(True, which="both", alpha=0.25) + flag_legend = ax.legend( + fontsize=7, + loc="best", + title=example["flag_key"], + title_fontsize=8, + ) + ax.text( + 0.02, + 0.03, + example["annotation"], + transform=ax.transAxes, + fontsize=7, + color="0.25", + ) +fig.legend( + handles=[ + Line2D([0], [0], marker="o", color="0.2", markerfacecolor="white", + lw=0, markersize=6, label="survives"), + Line2D([0], [0], marker="x", color="0.2", + lw=0, markersize=7, label="merges"), + Line2D([0], [0], marker="^", color="0.2", markerfacecolor="white", + lw=0, markersize=7, label="disrupts"), + Line2D([0], [0], color="0.35", ls=":", lw=1.4, label="CE/RLO"), + ], + loc="lower center", + ncol=5, + bbox_to_anchor=(0.5, -0.02), + fontsize=7, +) +plt.show()