From 28eab6232ee3a773b61be0c7a00368ec0b665e0f Mon Sep 17 00:00:00 2001
From: lukewilner <lwilner@andrew.cmu.edu>
Date: Wed, 24 Jun 2026 14:27:14 -0700
Subject: [PATCH] 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()