MovingUniverseLab · jluastro · Aug 13, 2026 · Aug 13, 2026 · Aug 13, 2026 · Aug 13, 2026
diff --git a/docs/figures/_offner_vs_spisea25.py b/docs/figures/_offner_vs_spisea25.py
@@ -0,0 +1,136 @@
+"""Shared layout and Table 1/2 data for Offner vs SPISEA v2.5 plots."""
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Patch
+from matplotlib.lines import Line2D
+
+
+_OFFNER_COLOR = '#8b3a2a'
+_TABLE1_COLOR = '#2f6db3'
+_BD_SHADE = '#e8d5b5'
+_FIGSIZE = (11.2, 4.6)
+_DPI = 160
+_BD_XLIM = (0.012, 0.20)
+_FULL_XLIM = (0.015, 40.0)
+
+# Table 1 γ_trunc (1–100 au / 1–102 au) with 1σ. Masses are geometric
+# means of the tabulated M1 intervals. The L/early-T γ=2.5 text value
+# is an interpolation knot only and is not plotted as a Table 1 point.
+_TABLE1_GAMMA = [
+    # name, M_lo, M_hi, gamma, gamma_err
+    ('Fontanive+2018', 0.019, 0.058, 4.8, 2.2),
+    ('Close+2003', 0.080, 0.095, 3.3, 1.2),
+    ('Allen+2007', 0.06, 0.15, 1.7, 0.5),
+    ('Winters mid-M', 0.15, 0.30, 0.7, 0.5),
+    ('Winters early-M', 0.3, 0.6, 0.1, 0.4),
+    ('Raghavan+2010', 0.75, 1.25, 0.2, 0.4),
+    ('De Rosa A', 1.6, 2.4, -1.3, 0.4),
+    ('MDS 3-5', 3.0, 5.0, -1.0, 0.5),
+    ('MDS 5-8', 5.0, 8.0, -1.7, 0.5),
+    ('MDS 8-17', 8.0, 17.0, -1.6, 0.5),
+    ('Sana O', 17.0, 50.0, -1.4, 0.4),
+]
+
+# Table 1 ã_all (au) with 1σ.
+_TABLE1_A_ALL = [
+    ('Fontanive+2018', 0.019, 0.058, 2.9, 1.1),
+    ('Close+2003', 0.080, 0.095, 3.7, 1.3),
+    ('Allen+2007', 0.06, 0.15, 6.9, 1.4),
+    ('Winters late-M', 0.075, 0.15, 3.9, 1.2),
+    ('Winters mid-M', 0.15, 0.30, 10.0, 3.0),
+    ('Winters early-M', 0.3, 0.6, 26.0, 4.0),
+    ('Raghavan+2010', 0.75, 1.25, 49.0, 6.0),
+    ('Tokovinin 2014b', 0.85, 1.5, 31.0, 5.0),
+    ('Moe & Kratter', 1.6, 2.4, 32.0, 8.0),
+    ('MDS 3-5', 3.0, 5.0, 28.0, 7.0),
+    ('MDS 5-8', 5.0, 8.0, 25.0, 7.0),
+    ('MDS 8-17', 8.0, 17.0, 23.0, 7.0),
+    ('Sana O', 17.0, 50.0, 19.0, 6.0),
+]
+
+# Table 2 lognormal μ (au) at the three published bins.
+_TABLE2_MU = [
+    ('late-M', 0.075, 0.15, 4.0),
+    ('early-M', 0.3, 0.6, 25.0),
+    ('FGK', 0.75, 1.25, 40.0),
+]
+
+
+def geom(lo, hi):
+    return float(np.sqrt(lo * hi))
+
+
+def table_xy(rows, y_idx=3, e_idx=4):
+    m = np.array([geom(r[1], r[2]) for r in rows])
+    y = np.array([r[y_idx] for r in rows], dtype=float)
+    err = np.array([r[e_idx] for r in rows], dtype=float)
+    return m, y, err
+
+
+def mean_q_from_gamma(gamma, q_min=0.01):
+    """⟨q⟩ for P(q) ∝ q^γ on [q_min, 1]."""
+    g = np.asarray(gamma, dtype=float)
+    qmin = float(q_min)
+    g_flat = np.atleast_1d(g).astype(float)
+    out_flat = np.empty(g_flat.shape, dtype=float)
+    near_m1 = np.abs(g_flat + 1.0) < 1e-12
+    near_m2 = np.abs(g_flat + 2.0) < 1e-12
+    ok = ~near_m1 & ~near_m2
+    if np.any(near_m1):
+        out_flat[near_m1] = (1.0 - qmin) / (-np.log(qmin))
+    if np.any(near_m2):
+        out_flat[near_m2] = -np.log(qmin) / (1.0 / qmin - 1.0)
+    if np.any(ok):
+        gp = g_flat[ok]
+        num = (1.0 - np.power(qmin, gp + 2.0)) / (gp + 2.0)
+        den = (1.0 - np.power(qmin, gp + 1.0)) / (gp + 1.0)
+        out_flat[ok] = num / den
+    out = out_flat.reshape(np.shape(g))
+    return float(out) if np.isscalar(gamma) else out
+
+
+def gamma_step_masses():
+    """Dense sampling so the 0.08 Msun γ step renders as a vertical jump."""
+    return np.concatenate([
+        np.logspace(np.log10(0.012), np.log10(0.07999), 300),
+        np.array([0.08, 0.08001]),
+        np.logspace(np.log10(0.081), np.log10(40.0), 300),
+    ])
+
+
+def bd_shade(ax, xlim):
+    ax.axvspan(xlim[0], 0.08, color=_BD_SHADE, alpha=0.55, zorder=0)
+    ax.axvline(0.08, color='#c4a574', ls='--', lw=1.2, zorder=1)
+
+
+def finish_panel(ax, xlim, ylim, title, ylabel, ylog=False):
+    ax.set_xscale('log')
+    if ylog:
+        ax.set_yscale('log')
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
+    ax.set_title(title, fontsize=11)
+    ax.set_xlabel(r'Primary mass $M_1$ ($M_\odot$)')
+    ax.set_ylabel(ylabel)
+    ax.tick_params(which='both', direction='in', top=True, right=True)
+
+
+def two_axes(suptitle):
+    fig, axes = plt.subplots(1, 2, figsize=_FIGSIZE, gridspec_kw={'wspace': 0.28})
+    fig.suptitle(suptitle, fontsize=13, y=1.02)
+    return fig, axes
+
+
+def save(fig, filename):
+    out = os.path.join(os.path.dirname(os.path.abspath(__file__)), filename)
+    fig.savefig(out, dpi=_DPI, bbox_inches='tight', facecolor='white')
+    plt.close(fig)
+    print('Wrote', out)
+    return out
+
+
+# Re-export legend artists so plot scripts need one import.
+Patch = Patch
+Line2D = Line2D
diff --git a/docs/figures/csf_offner_vs_spisea2.5.png b/docs/figures/csf_offner_vs_spisea2.5.png
diff --git a/docs/figures/meanq_offner_vs_spisea2.5.png b/docs/figures/meanq_offner_vs_spisea2.5.png
diff --git a/docs/figures/mf_offner_vs_spisea2.5.png b/docs/figures/mf_offner_vs_spisea2.5.png
diff --git a/docs/figures/plot_csf_offner_vs_spisea2.5.py b/docs/figures/plot_csf_offner_vs_spisea2.5.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+"""
+Generate docs/figures/csf_offner_vs_spisea2.5.png
+
+Two-panel comparison of companion star fraction vs primary mass:
+SPISEA v2.5 ``MultiplicityUnresolved.companion_star_fraction``,
+Offner et al. 2023 logistic in log-mass, and Offner Table 1 CF
+points.
+
+Run from the repository root::
+
+    python docs/figures/plot_csf_offner_vs_spisea2.5.py
+"""
+import os
+import sys
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.patches import Patch
+from matplotlib.lines import Line2D
+
+# Allow running without installing the package.
+_REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__),
+                                          '..', '..'))
+if _REPO_ROOT not in sys.path:
+    sys.path.insert(0, _REPO_ROOT)
+
+from spisea.imf import multiplicity
+
+
+# Offner et al. 2023 Table 1: (M_lo, M_hi, CF)
+_TABLE1_CF = [
+    (0.019, 0.058, 0.08),
+    (0.05, 0.08, 0.16),
+    (0.080, 0.095, 0.19),
+    (0.06, 0.15, 0.20),
+    (0.075, 0.15, 0.21),
+    (0.15, 0.30, 0.27),
+    (0.3, 0.6, 0.38),
+    (0.75, 1.25, 0.60),
+    (0.85, 1.5, 0.62),
+    (1.6, 2.4, 0.99),
+    (3.0, 5.0, 1.28),
+    (5.0, 8.0, 1.55),
+    (8.0, 17.0, 1.80),
+    (17.0, 50.0, 2.10),
+]
+
+
+def _table1_xy():
+    m = np.array([np.sqrt(lo * hi) for lo, hi, _ in _TABLE1_CF])
+    cf = np.array([row[2] for row in _TABLE1_CF])
+    return m, cf
+
+
+def _style_panel(ax, m_off, csf_off, m_lu, csf_lu, m_tab, cf_tab,
+                 xlim, ylim, title):
+    ax.axvspan(xlim[0], 0.08, color='#e8d5b5', alpha=0.55, zorder=0)
+    ax.axvline(0.08, color='#c4a574', ls='--', lw=1.2, zorder=1)
+    ax.plot(m_lu, csf_lu, color='0.25', ls='--', lw=1.6, zorder=3,
+            label=r'SPISEA v2.5 $0.50\,M^{0.45}$')
+    ax.plot(m_off, csf_off, color='#8b3a2a', ls='-', lw=2.4, zorder=4,
+            label='Offner logistic in log M')
+    ax.plot(m_tab, cf_tab, 'o', color='#2f6db3', ms=5.5, mfc='white',
+            mew=1.3, zorder=5, label='Offner Table 1 CF')
+    ax.set_xscale('log')
+    ax.set_xlim(*xlim)
+    ax.set_ylim(*ylim)
+    ax.set_title(title, fontsize=11)
+    ax.set_xlabel(r'Primary mass $M_1$ ($M_\odot$)')
+    ax.set_ylabel('Companion star fraction')
+    ax.tick_params(which='both', direction='in', top=True, right=True)
+
+
+def main():
+    offner = multiplicity.MultiplicityUnresolvedOffner2023()
+    lu = multiplicity.MultiplicityUnresolved()
+
+    m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800)
+    csf_off = offner.companion_star_fraction(m_wide)
+    csf_lu = lu.companion_star_fraction(m_wide)
+    m_tab, cf_tab = _table1_xy()
+
+    fig, axes = plt.subplots(1, 2, figsize=(11.2, 4.6),
+                             gridspec_kw={'wspace': 0.28})
+    fig.suptitle('Offner 2023 vs SPISEA v2.5: companion star fraction',
+                 fontsize=13, y=1.02)
+
+    _style_panel(
+        axes[0], m_wide, csf_off, m_wide, csf_lu, m_tab, cf_tab,
+        xlim=(0.012, 0.20), ylim=(0.0, 0.50),
+        title='Brown-dwarf regime')
+    _style_panel(
+        axes[1], m_wide, csf_off, m_wide, csf_lu, m_tab, cf_tab,
+        xlim=(0.015, 20.0), ylim=(0.0, 3.0),
+        title='BD through early B')
+
+    legend_handles = [
+        Line2D([0], [0], color='0.25', ls='--', lw=1.6,
+               label=r'SPISEA v2.5 $0.50\,M^{0.45}$'),
+        Line2D([0], [0], color='#8b3a2a', ls='-', lw=2.4,
+               label='Offner logistic in log M'),
+        Line2D([0], [0], marker='o', color='#2f6db3', ls='none',
+               mfc='white', mew=1.3, ms=6, label='Offner Table 1 CF'),
+        Patch(facecolor='#e8d5b5', edgecolor='none', alpha=0.8,
+              label=r'BD ($M\leq 0.08$)'),
+    ]
+    axes[0].legend(handles=legend_handles, loc='upper left', fontsize=8,
+                   frameon=True, fancybox=False, edgecolor='0.7')
+
+    out = os.path.join(os.path.dirname(__file__),
+                       'csf_offner_vs_spisea2.5.png')
+    fig.savefig(out, dpi=160, bbox_inches='tight', facecolor='white')
+    plt.close(fig)
+    print('Wrote', out)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/docs/figures/plot_meanq_offner_vs_spisea2.5.py b/docs/figures/plot_meanq_offner_vs_spisea2.5.py
@@ -0,0 +1,68 @@
+#!/usr/bin/env python
+"""
+Generate docs/figures/meanq_offner_vs_spisea2.5.png
+
+Mean q implied by ``q_power_at_mass`` on Offner vs SPISEA v2.5
+``MultiplicityUnresolved``.
+
+Run from the repository root::
+
+    python docs/figures/plot_meanq_offner_vs_spisea2.5.py
+"""
+import os
+import sys
+
+import numpy as np
+
+_REPO_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..'))
+if _REPO_ROOT not in sys.path:
+    sys.path.insert(0, _REPO_ROOT)
+_FIGDIR = os.path.dirname(os.path.abspath(__file__))
+if _FIGDIR not in sys.path:
+    sys.path.insert(0, _FIGDIR)
+
+from spisea.imf import multiplicity
+from _offner_vs_spisea25 import (
+    Line2D, Patch, _BD_SHADE, _BD_XLIM, _FULL_XLIM, _OFFNER_COLOR,
+    bd_shade, finish_panel, gamma_step_masses, mean_q_from_gamma, save,
+    two_axes,
+)
+
+
+def main():
+    offner = multiplicity.MultiplicityUnresolvedOffner2023()
+    lu = multiplicity.MultiplicityUnresolved()
+
+    m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800)
+    q_off = mean_q_from_gamma(offner.q_power_at_mass(m_wide),
+                             q_min=offner.q_min)
+    m_step = gamma_step_masses()
+    q_lu = mean_q_from_gamma(lu.q_power_at_mass(m_step), q_min=lu.q_min)
+
+    fig, axes = two_axes(
+        r'Offner 2023 vs SPISEA v2.5: mean mass ratio $\langle q\rangle$')
+    ylabel = r'$\langle q\rangle$  on $[0.01,\,1]$'
+    for ax, xlim, title in (
+        (axes[0], _BD_XLIM, 'Brown-dwarf regime'),
+        (axes[1], _FULL_XLIM, 'BD through O'),
+    ):
+        bd_shade(ax, xlim)
+        ax.plot(m_step, q_lu, color='0.25', ls='--', lw=1.6, zorder=3)
+        ax.plot(m_wide, q_off, color=_OFFNER_COLOR, ls='-', lw=2.4, zorder=4)
+        finish_panel(ax, xlim, (0.0, 1.0), title, ylabel)
+
+    legend_handles = [
+        Line2D([0], [0], color='0.25', ls='--', lw=1.6,
+               label=r'SPISEA v2.5 from $\gamma$ step'),
+        Line2D([0], [0], color=_OFFNER_COLOR, ls='-', lw=2.4,
+               label=r'Offner  from $\gamma(M)$ logistic'),
+        Patch(facecolor=_BD_SHADE, edgecolor='none', alpha=0.8,
+              label=r'BD ($M\leq 0.08$)'),
+    ]
+    axes[1].legend(handles=legend_handles, loc='upper right', fontsize=8,
+                   frameon=True, fancybox=False, edgecolor='0.7')
+    save(fig, 'meanq_offner_vs_spisea2.5.png')
+
+
+if __name__ == '__main__':
+    main()