5 new config. figures (pts, neta, hewind, zsun, and windflag)

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

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