diff --git a/docs/figures/_offner_vs_spisea25.py b/docs/figures/_offner_vs_spisea25.py new file mode 100644 index 00000000..eba01e87 --- /dev/null +++ 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 new file mode 100644 index 00000000..d9207772 Binary files /dev/null and b/docs/figures/csf_offner_vs_spisea2.5.png differ diff --git a/docs/figures/meanq_offner_vs_spisea2.5.png b/docs/figures/meanq_offner_vs_spisea2.5.png new file mode 100644 index 00000000..b00e27a2 Binary files /dev/null and b/docs/figures/meanq_offner_vs_spisea2.5.png differ diff --git a/docs/figures/mf_offner_vs_spisea2.5.png b/docs/figures/mf_offner_vs_spisea2.5.png new file mode 100644 index 00000000..20da5170 Binary files /dev/null and b/docs/figures/mf_offner_vs_spisea2.5.png differ diff --git a/docs/figures/plot_csf_offner_vs_spisea2.5.py b/docs/figures/plot_csf_offner_vs_spisea2.5.py new file mode 100644 index 00000000..685973b9 --- /dev/null +++ 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 new file mode 100644 index 00000000..9891dc4b --- /dev/null +++ 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() diff --git a/docs/figures/plot_mf_offner_vs_spisea2.5.py b/docs/figures/plot_mf_offner_vs_spisea2.5.py new file mode 100644 index 00000000..b4999d13 --- /dev/null +++ b/docs/figures/plot_mf_offner_vs_spisea2.5.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python +""" +Generate docs/figures/mf_offner_vs_spisea2.5.png + +Two-panel comparison of multiplicity fraction vs primary mass: +SPISEA v2.5 ``MultiplicityUnresolved`` (array power law and +scalar BD staircase), Offner et al. 2023 logistic in log-mass, +and Offner Table 1 points with error bars. + +Run from the repository root:: + + python docs/figures/plot_mf_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, MF, MF_err) +_TABLE1 = [ + (0.019, 0.058, 0.08, 0.06), + (0.05, 0.08, 0.15, 0.04), + (0.080, 0.095, 0.19, 0.07), + (0.06, 0.15, 0.20, 0.04), + (0.075, 0.15, 0.19, 0.03), + (0.15, 0.30, 0.23, 0.02), + (0.3, 0.6, 0.30, 0.02), + (0.75, 1.25, 0.46, 0.03), + (0.85, 1.5, 0.47, 0.03), + (1.6, 2.4, 0.68, 0.07), + (3.0, 5.0, 0.81, 0.06), + (5.0, 8.0, 0.89, 0.05), + (8.0, 17.0, 0.93, 0.04), + (17.0, 50.0, 0.96, 0.04), +] + + +def _table1_xy(): + m = np.array([np.sqrt(lo * hi) for lo, hi, _, _ in _TABLE1]) + mf = np.array([row[2] for row in _TABLE1]) + err = np.array([row[3] for row in _TABLE1]) + return m, mf, err + + +def _style_panel(ax, m_off, mf_off, m_lu, mf_lu, m_step, mf_step, + m_tab, mf_tab, err_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, mf_lu, color='0.25', ls='--', lw=1.6, zorder=3, + label=r'SPISEA v2.5 $0.44\,M^{0.51}$') + ax.plot(m_step, mf_step, color='0.45', ls=':', lw=1.8, zorder=3, + label=r'SPISEA v2.5 scalar BD bins (0 / 8% / 16%)') + ax.plot(m_off, mf_off, color='#8b3a2a', ls='-', lw=2.4, zorder=4, + label='Offner logistic in log M') + ax.errorbar(m_tab, mf_tab, yerr=err_tab, fmt='o', color='#2f6db3', + ms=5.5, mfc='white', mew=1.3, elinewidth=1.1, capsize=2.5, + zorder=5, label='Offner Table 1') + 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('Multiplicity 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) + mf_off = offner.multiplicity_fraction(m_wide) + # Array path: stellar power law only (no BD staircase). + mf_lu = lu.multiplicity_fraction(m_wide) + + # Staircase sampled densely so the steps render as vertical jumps. + # Scalar path applies the Aberasturi/Fontanive BD bins. + m_step = np.concatenate([ + np.array([0.012, 0.01999, 0.02001, 0.05999, 0.06001, 0.07999, 0.08001]), + np.logspace(np.log10(0.081), np.log10(40.0), 200), + ]) + mf_step = np.array([ + float(lu.multiplicity_fraction(float(m))) for m in m_step]) + + m_tab, mf_tab, err_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: multiplicity fraction', + fontsize=13, y=1.02) + + _style_panel( + axes[0], m_wide, mf_off, m_wide, mf_lu, m_step, mf_step, + m_tab, mf_tab, err_tab, + xlim=(0.012, 0.20), ylim=(0.0, 0.45), + title='Brown-dwarf regime') + _style_panel( + axes[1], m_wide, mf_off, m_wide, mf_lu, m_step, mf_step, + m_tab, mf_tab, err_tab, + xlim=(0.015, 20.0), ylim=(0.0, 1.0), + title='BD through early B') + + legend_handles = [ + Line2D([0], [0], color='0.25', ls='--', lw=1.6, + label=r'SPISEA v2.5 $0.44\,M^{0.51}$'), + Line2D([0], [0], color='0.45', ls=':', lw=1.8, + label=r'SPISEA v2.5 scalar BD bins (0 / 8% / 16%)'), + 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'), + 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__), 'mf_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_q_offner_vs_spisea2.5.py b/docs/figures/plot_q_offner_vs_spisea2.5.py new file mode 100644 index 00000000..a67fd204 --- /dev/null +++ b/docs/figures/plot_q_offner_vs_spisea2.5.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python +""" +Generate docs/figures/q_offner_vs_spisea2.5.png + +Offner ``q_power_at_mass`` vs SPISEA v2.5 ``MultiplicityUnresolved``. + +Run from the repository root:: + + python docs/figures/plot_q_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, + _TABLE1_COLOR, _TABLE1_GAMMA, bd_shade, finish_panel, gamma_step_masses, + save, table_xy, two_axes, +) + + +def main(): + offner = multiplicity.MultiplicityUnresolvedOffner2023() + lu = multiplicity.MultiplicityUnresolved() + + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + g_off = offner.q_power_at_mass(m_wide) + m_step = gamma_step_masses() + g_lu = lu.q_power_at_mass(m_step) + m_tab, g_tab, err_tab = table_xy(_TABLE1_GAMMA) + + fig, axes = two_axes(r'Offner 2023 vs SPISEA v2.5: mass-ratio index $\gamma$') + ylabel = r'$\gamma$ ($P(q)\propto q^{\gamma}$)' + for ax, xlim, title in ( + (axes[0], _BD_XLIM, 'Brown-dwarf regime'), + (axes[1], _FULL_XLIM, 'BD through O'), + ): + bd_shade(ax, xlim) + ax.axhline(0.0, color='0.55', ls=':', lw=1.1, zorder=2) + ax.plot(m_step, g_lu, color='0.25', ls='--', lw=1.6, zorder=3) + ax.plot(m_wide, g_off, color=_OFFNER_COLOR, ls='-', lw=2.4, zorder=4) + ax.errorbar(m_tab, g_tab, yerr=err_tab, fmt='o', color=_TABLE1_COLOR, + ms=5.5, mfc='white', mew=1.3, elinewidth=1.1, capsize=2.5, + zorder=5) + finish_panel(ax, xlim, (-2.2, 7.0), title, ylabel) + + legend_handles = [ + Line2D([0], [0], color='0.25', ls='--', lw=1.6, + label=r'SPISEA v2.5 ($\gamma=6.1$ / $-0.4$)'), + Line2D([0], [0], color=_OFFNER_COLOR, ls='-', lw=2.4, + label=r'Offner err-wt logistic in log $M$'), + Line2D([0], [0], marker='o', color=_TABLE1_COLOR, ls='none', + mfc='white', mew=1.3, ms=6, label=r'Table 1 $\gamma_\mathrm{trunc}$'), + 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, 'q_offner_vs_spisea2.5.png') + + +if __name__ == '__main__': + main() diff --git a/docs/figures/plot_q_sep_offner_vs_spisea2.5.py b/docs/figures/plot_q_sep_offner_vs_spisea2.5.py new file mode 100644 index 00000000..ae97b497 --- /dev/null +++ b/docs/figures/plot_q_sep_offner_vs_spisea2.5.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python +""" +Write the q / separation / σ / mean-q comparison figures vs SPISEA v2.5. + +This is a thin runner. Each PNG is generated by its matching script: + + docs/figures/plot_q_offner_vs_spisea2.5.py + docs/figures/plot_sep_offner_vs_spisea2.5.py + docs/figures/plot_sig_loga_offner_vs_spisea2.5.py + docs/figures/plot_meanq_offner_vs_spisea2.5.py + +Run from the repository root:: + + python docs/figures/plot_q_sep_offner_vs_spisea2.5.py +""" +import os +import runpy + + +_SCRIPTS = ( + 'plot_q_offner_vs_spisea2.5.py', + 'plot_sep_offner_vs_spisea2.5.py', + 'plot_sig_loga_offner_vs_spisea2.5.py', + 'plot_meanq_offner_vs_spisea2.5.py', +) + + +def main(): + here = os.path.dirname(os.path.abspath(__file__)) + for name in _SCRIPTS: + runpy.run_path(os.path.join(here, name), run_name='__main__') + + +if __name__ == '__main__': + main() diff --git a/docs/figures/plot_sep_offner_vs_spisea2.5.py b/docs/figures/plot_sep_offner_vs_spisea2.5.py new file mode 100644 index 00000000..2f7f100f --- /dev/null +++ b/docs/figures/plot_sep_offner_vs_spisea2.5.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python +""" +Generate docs/figures/sep_offner_vs_spisea2.5.png + +Offner ``a_mean`` vs SPISEA v2.5 ``MultiplicityResolvedDK.a_mean``. + +Run from the repository root:: + + python docs/figures/plot_sep_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, + _TABLE1_A_ALL, _TABLE1_COLOR, _TABLE2_MU, bd_shade, finish_panel, geom, + save, table_xy, two_axes, +) + + +def main(): + resolved = multiplicity.MultiplicityResolvedOffner2023() + dk = multiplicity.MultiplicityResolvedDK() + + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + a_off = resolved.a_mean(m_wide) + a_lu = dk.a_mean(m_wide) + m_tab, a_tab, err_tab = table_xy(_TABLE1_A_ALL) + m_t2 = np.array([geom(r[1], r[2]) for r in _TABLE2_MU]) + a_t2 = np.array([r[3] for r in _TABLE2_MU], dtype=float) + + fig, axes = two_axes( + r'Offner 2023 vs SPISEA v2.5: characteristic separation') + ylabel = r'$\mu(a)$ (AU)' + 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_wide, a_lu, color='0.25', ls='--', lw=1.6, zorder=3) + ax.plot(m_wide, a_off, color=_OFFNER_COLOR, ls='-', lw=2.4, zorder=4) + ax.errorbar(m_tab, a_tab, yerr=err_tab, fmt='o', color=_TABLE1_COLOR, + ms=5.5, mfc='white', mew=1.3, elinewidth=1.1, capsize=2.5, + zorder=5) + ax.plot(m_t2, a_t2, 's', color=_OFFNER_COLOR, mfc=_OFFNER_COLOR, + ms=7, zorder=6, mew=0.6) + finish_panel(ax, xlim, (1.0, 400.0), title, ylabel, ylog=True) + + legend_handles = [ + Line2D([0], [0], color='0.25', ls='--', lw=1.6, + label=r'SPISEA v2.5 mean $a$'), + Line2D([0], [0], color=_OFFNER_COLOR, ls='-', lw=2.4, + label=r'Offner smooth broken PL ($s=0.1$ dex)'), + Line2D([0], [0], marker='o', color=_TABLE1_COLOR, ls='none', + mfc='white', mew=1.3, ms=6, label=r'Table 1 $\tilde{a}_\mathrm{all}$'), + Line2D([0], [0], marker='s', color=_OFFNER_COLOR, ls='none', + mfc=_OFFNER_COLOR, ms=7, label=r'Table 2 $\mu$ knots'), + Patch(facecolor=_BD_SHADE, 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') + save(fig, 'sep_offner_vs_spisea2.5.png') + + +if __name__ == '__main__': + main() diff --git a/docs/figures/plot_sig_loga_offner_vs_spisea2.5.py b/docs/figures/plot_sig_loga_offner_vs_spisea2.5.py new file mode 100644 index 00000000..dc1da8d0 --- /dev/null +++ b/docs/figures/plot_sig_loga_offner_vs_spisea2.5.py @@ -0,0 +1,71 @@ +#!/usr/bin/env python +""" +Generate docs/figures/sig_loga_offner_vs_spisea2.5.png + +Offner ``sigma_log_a`` vs SPISEA v2.5 +``MultiplicityResolvedDK.sigma_log_a``. + +Run from the repository root:: + + python docs/figures/plot_sig_loga_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, save, two_axes, +) + + +def main(): + resolved = multiplicity.MultiplicityResolvedOffner2023() + dk = multiplicity.MultiplicityResolvedDK() + + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + sig_off = resolved.sigma_log_a(m_wide) + sig_lu = dk.sigma_log_a(m_wide) + m_t2 = np.array(resolved.sep_sig_mass, dtype=float) + sig_t2 = np.array(resolved.sep_sig, dtype=float) + + fig, axes = two_axes( + r'Offner 2023 vs SPISEA v2.5: $\sigma(\log_{10} a)$') + ylabel = r'$\sigma(\log_{10} a)$' + 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_wide, sig_lu, color='0.25', ls='--', lw=1.6, zorder=3) + ax.plot(m_wide, sig_off, color=_OFFNER_COLOR, ls='-', lw=2.4, zorder=4) + ax.plot(m_t2, sig_t2, 's', color=_OFFNER_COLOR, mfc=_OFFNER_COLOR, + ms=7, zorder=6, mew=0.6) + finish_panel(ax, xlim, (0.0, 2.05), title, ylabel) + + legend_handles = [ + Line2D([0], [0], color='0.25', ls='--', lw=1.6, + label=r'SPISEA v2.5 DK $\sigma_{\log a}$'), + Line2D([0], [0], color=_OFFNER_COLOR, ls='-', lw=2.4, + label=r'Offner 2-param logistic $\sigma$'), + Line2D([0], [0], marker='s', color=_OFFNER_COLOR, ls='none', + mfc=_OFFNER_COLOR, ms=7, label=r'Table 2 $\sigma$ knots'), + Patch(facecolor=_BD_SHADE, edgecolor='none', alpha=0.8, + label=r'BD ($M\leq 0.08$)'), + ] + axes[1].legend(handles=legend_handles, loc='upper left', fontsize=8, + frameon=True, fancybox=False, edgecolor='0.7') + save(fig, 'sig_loga_offner_vs_spisea2.5.png') + + +if __name__ == '__main__': + main() diff --git a/docs/figures/q_offner_vs_spisea2.5.png b/docs/figures/q_offner_vs_spisea2.5.png new file mode 100644 index 00000000..059743bd Binary files /dev/null and b/docs/figures/q_offner_vs_spisea2.5.png differ diff --git a/docs/figures/sep_offner_vs_spisea2.5.png b/docs/figures/sep_offner_vs_spisea2.5.png new file mode 100644 index 00000000..e34f2ba3 Binary files /dev/null and b/docs/figures/sep_offner_vs_spisea2.5.png differ diff --git a/docs/figures/sig_loga_offner_vs_spisea2.5.png b/docs/figures/sig_loga_offner_vs_spisea2.5.png new file mode 100644 index 00000000..792d2433 Binary files /dev/null and b/docs/figures/sig_loga_offner_vs_spisea2.5.png differ diff --git a/docs/imf.rst b/docs/imf.rst index ed601775..91c2946a 100644 --- a/docs/imf.rst +++ b/docs/imf.rst @@ -17,7 +17,10 @@ and exponents of the IMF. The IMF object is an input for the :ref:`cluster_objects`, and will be used to draw the inital stellar mass distribution for the cluster. A :ref:`multi_obj` may be passed to the IMF object to -form multiple star systems. +form multiple star systems. The default is the SPISEA v2.5 +:class:`~imf.multiplicity.MultiplicityUnresolved` / +:class:`~imf.multiplicity.MultiplicityResolvedDK` objects; +Offner et al. 2023 is opt-in (see :ref:`multi_obj`). Base IMF Class -------------- diff --git a/docs/index.rst b/docs/index.rst index 9a0f7408..fdfe0055 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -97,6 +97,8 @@ Change Log * New tutorial for creating SPISEA clusters using COSMIC: docs/Cluster_w_COSMIC.ipynb * *Minor Changes* + * Added opt-in Offner et al. 2023 (Protostars and Planets VII; arXiv:2203.10066) multiplicity model, including brown dwarfs. MF/CSF is a logistic in log-mass; :math:`\gamma(M)` is an error-weighted logistic; :math:`\mu(a)` is a smooth broken power law (:math:`s=0.1` dex); :math:`\sigma(\log_{10} a)` is a 2-parameter logistic. Comparison figures vs SPISEA v2.5 are in :ref:`multi_obj`. The SPISEA v2.5 :class:`~imf.multiplicity.MultiplicityUnresolved` / :class:`~imf.multiplicity.MultiplicityResolvedDK` objects remain the default. + * Companion mass and separation draws (including brown-dwarf q and binaries-only BD systems) now live on the multiplicity object rather than being special-cased in ``imf.py``. * The MISTv1.2-synthpop model extension was modified to include denser sampling in the gap between the base MISTv1.2 grids and 0.1Msun. * Modified default for MISTv1.2 isochrones: synthpop_extension will be True by default to keep a consistent lower mass limit of 0.1Msun across all ages and metallicities. * Added option to return synthetic photometry in terms of AB or ST mag units in IsochronePhot. Vega mag units remains the default. New meta keyword `MAGSYS` added to output tables to specify magnitude unit system. diff --git a/docs/multiplicity.rst b/docs/multiplicity.rst index 54177383..4bef07bc 100644 --- a/docs/multiplicity.rst +++ b/docs/multiplicity.rst @@ -7,13 +7,25 @@ The properties of multiple systems in the stellar population is defined by the stellar multiplicity object. The multiplicity classes are defined in ``spisea/imf/multiplicity.py``. -To call a multiplicity class:: +To call a multiplicity class and wire it into a cluster:: - from spisea.imf import multiplicity - multi_obj = multiplicity. + from spisea.imf import imf, multiplicity + from spisea import synthetic -The multiplicity object is an input for the :ref:`imf_objects`, as it -impacts how the stellar masses are drawn. + multi = multiplicity.(...) + imf_obj = imf.Kroupa_2001(multiplicity=multi) + cluster = synthetic.ResolvedCluster(iso, imf_obj, Mcl) + +The multiplicity object provides the following functions used by the IMF: + +* ``multiplicity_fraction(mass)`` +* ``companion_star_fraction(mass)`` +* ``random_q(x, mass=None)`` — pass ``mass`` for mass-dependent q + (brown-dwarf vs stellar). ``random_q(x)`` with no mass assumes a uniform distribution. +* ``random_companion_count(x, CSF, MF, mass=None, rng=None)`` — + companion-count draw. If ``companion_max`` is True, counts are + capped at ``CSF_max`` at all masses. +* attributes ``companion_max``, ``CSF_max``, ``q_min`` The user can choose either an unresolved or a resolved multiplicity object. If a resolved @@ -27,21 +39,102 @@ returned in the ``star_systems`` table off the cluster object is the same for both unresolved and resolved multiplicity classes: it represents the combined photometry of all stars within a given system. -For most selected evolution models, the multiples are evolved as single stars. -To evolve binaries (does not support higher order multiples), you should use one of the ``MultiplicityResolved`` classes +For most selected evolution models, the companions are evolved as single stars. +To evolve binaries with mass exchange (does not support triples and higher order multiples), +you should use one of the ``MultiplicityResolved`` classes and the ``COSMIC`` evolution model. See the example jupyter notebook `Cluster_w_COSMIC.ipynb `_ for an example. -Note that currently COSMIC due to being external evolution is significantly slower than the other evolution options. +Note that currently COSMIC, due to being external evolution is significantly slower (2-10x) than the other evolution options. + +The default is the SPISEA v2.5 +:class:`~imf.multiplicity.MultiplicityUnresolved` / +:class:`~imf.multiplicity.MultiplicityResolvedDK` pair. Offner et al. +2023 is opt-in; see those class docstrings. Unresolved Multiplicity Classes ------------------------------------------ .. autoclass:: imf.multiplicity.MultiplicityUnresolved + :show-inheritance: :members: companion_star_fraction, - multiplicity_fraction, random_q + multiplicity_fraction, random_q, + q_power_at_mass, random_companion_count + +.. autoclass:: imf.multiplicity.MultiplicityPiecewisePowerLaw + :show-inheritance: + :members: multiplicity_fraction, companion_star_fraction + +.. autoclass:: imf.multiplicity.MultiplicityLogistic + :show-inheritance: + :members: multiplicity_fraction, companion_star_fraction + +.. autoclass:: imf.multiplicity.MultiplicityUnresolvedOffner2023 + :show-inheritance: + :members: multiplicity_fraction, companion_star_fraction, + q_power_at_mass, random_q, log_a_mean, a_mean, + sigma_log_a Resolved Multiplicity Classes ------------------------------------------ .. autoclass:: imf.multiplicity.MultiplicityResolvedDK :show-inheritance: + :members: log_semimajoraxis, log_a_mean, a_mean, sigma_log_a + +.. autoclass:: imf.multiplicity.MultiplicityResolvedOffner2023 + :show-inheritance: + :members: log_semimajoraxis, log_a_mean, a_mean, sigma_log_a + + +Comparison figures +------------------------------------------ +Offner 2023 vs SPISEA v2.5. Model details are on the class +docstrings above. + +.. figure:: figures/mf_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: multiplicity fraction vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: multiplicity fraction vs primary mass. + +.. figure:: figures/csf_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: companion star fraction vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: companion star fraction vs primary mass. + +.. figure:: figures/q_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: mass-ratio index vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: mass-ratio index :math:`\gamma` vs primary mass. + +.. figure:: figures/meanq_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: mean mass ratio vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: mean mass ratio :math:`\langle q \rangle` vs primary mass. + +.. figure:: figures/sep_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: characteristic separation vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: characteristic :math:`\mu(a)` vs primary mass. + +.. figure:: figures/sig_loga_offner_vs_spisea2.5.png + :alt: Offner 2023 vs SPISEA v2.5: sigma of log10 a vs primary mass + :align: center + + Offner 2023 vs SPISEA v2.5: :math:`\sigma(\log_{10} a)` vs primary mass. + +From the repository root:: + + python docs/figures/plot_mf_offner_vs_spisea2.5.py + python docs/figures/plot_csf_offner_vs_spisea2.5.py + python docs/figures/plot_q_offner_vs_spisea2.5.py + python docs/figures/plot_meanq_offner_vs_spisea2.5.py + python docs/figures/plot_sep_offner_vs_spisea2.5.py + python docs/figures/plot_sig_loga_offner_vs_spisea2.5.py + +``python docs/figures/plot_q_sep_offner_vs_spisea2.5.py`` writes the +last four PNGs. diff --git a/docs/paper_examples/Begbie+26/Figure 12.ipynb b/docs/paper_examples/Begbie+26/Figure 12.ipynb index 929b6a92..9ab4fe9f 100644 --- a/docs/paper_examples/Begbie+26/Figure 12.ipynb +++ b/docs/paper_examples/Begbie+26/Figure 12.ipynb @@ -56,7 +56,8 @@ ], "source": [ "# semimajor axis distribution\n", - "log_a = np.array([mult.log_semimajoraxis(m) for m in masses])\n", + "rng = np.random.default_rng()\n", + "log_a = np.array([mult.log_semimajoraxis(m, rng=rng) for m in masses])\n", "\n", "plt.figure()\n", "plt.scatter(masses, 10**log_a, s=1, alpha=0.3)\n", diff --git a/spisea/imf/imf.py b/spisea/imf/imf.py index 241e04d9..fbc3a445 100755 --- a/spisea/imf/imf.py +++ b/spisea/imf/imf.py @@ -249,62 +249,13 @@ def generate_cluster(self, totalMass): def calc_multi(self, newMasses, newIsMultiple, CSF, MF): """ Helper function to calculate multiples more efficiently. - We will use array operations as much as possible. - Uses Fontanive+18 parameters for brown dwarf masses - (M <= 0.08 M_sun) while keeping default parameters for - all other stellar primaries. + Companion counts and companion-mass draws (including brown-dwarf + q distributions and binaries-only BD systems) are delegated to + the multiplicity object. """ - # Copy over the primary masses. Eventually add the companions. - newSystemMasses = newMasses.copy() - - # Identify multiple systems, calculate number of companions for each - multiple_idx = np.where(newIsMultiple)[0] - comp_nums = 1 + self.rng.poisson((CSF[multiple_idx] / MF[multiple_idx]) - 1) - if self._multi_props.companion_max: - too_many = np.where(comp_nums > self._multi_props.CSF_max)[0] - comp_nums[too_many] = self._multi_props.CSF_max - primary = newMasses[multiple_idx] - - # limit BD primaries to 1 companion (Fontanive+18) - bd_mask = primary <= 0.08 - comp_nums[bd_mask & (comp_nums > 1)] = 1 - - # We will deal with each number of multiple system independently. This is - # so we can put in uniform arrays in _multi_props.random_q. - comp_unique = np.unique(comp_nums) - comp_indices = [np.where(comp_nums == i)[0] for i in comp_unique] - if np.any(newIsMultiple): - compMasses = np.zeros((len(newMasses), max(comp_unique))) - else: - compMasses = np.zeros((len(newMasses), 1)) - - for comp_num, comp_index in zip(comp_unique, comp_indices): - prim_subset = primary[comp_index] - bd_sub_mask = prim_subset <= 0.08 - star_sub_mask = ~bd_sub_mask - - q_values = np.empty((len(comp_index), comp_num)) - - # Stellar primaries: use default Duchene & Kraus distribution - if np.any(star_sub_mask): - q_values[star_sub_mask] = self._multi_props.random_q(self.rng.random((star_sub_mask.sum(), comp_num))) - - # BD primaries: use Fontanive+18 power-law distribution - if np.any(bd_sub_mask): - b = 1.0 + 6.1 # gamma from Fontanive+18 - rand_vals = self.rng.random((bd_sub_mask.sum(), comp_num)) - q_values[bd_sub_mask] = (rand_vals * (1.0 - self._multi_props.q_min ** b) + - self._multi_props.q_min ** b) ** (1.0 / b) - - m_comp = np.multiply(q_values, np.transpose([prim_subset])) - compMasses[multiple_idx[comp_index], :comp_num] = m_comp - - # Mask out companions below the minimum mass - compMasses = np.ma.MaskedArray(compMasses, mask=compMasses < self._mass_limits[0]) - newSystemMasses[multiple_idx] += compMasses[multiple_idx].sum(axis=1) - newIsMultiple = np.any(~compMasses.mask, axis=1) - - return compMasses, newSystemMasses, newIsMultiple + return self._multi_props.draw_companion_masses( + newMasses, newIsMultiple, CSF, MF, + rng=self.rng, mass_min=self._mass_limits[0]) class IMF_broken_powerlaw(IMF): diff --git a/spisea/imf/multiplicity.py b/spisea/imf/multiplicity.py index 4258cb4c..8ea79c81 100755 --- a/spisea/imf/multiplicity.py +++ b/spisea/imf/multiplicity.py @@ -1,6 +1,5 @@ import numpy as np import astropy.modeling -from random import choice from scipy.stats import truncnorm defaultMF_amp = 0.44 @@ -13,104 +12,186 @@ default_aMean = 100.0 # log (AU) default_aSigma = 0.1 # log (AU) +# Hydrogen-burning limit used for BD-primary (binaries-only) logic. +# Offner et al. 2023 use M_comp > 0.075 Msun as the MS companion cut; +# SPISEA keeps 0.08 Msun for consistency with existing BD handling. +H_BURNING_MASS = 0.08 + # Eventually we should add in separation properties. (a_mean, a_sigma) + +class _ResolvedOrbitalMixin(object): + """Eccentricity and Keplerian angles shared by resolved multiplicity classes.""" + + def random_e(self, x): + """ + Generate random eccentricity from the inverse of the CDF where the PDF is f(e) = 2e from Duchene and Kraus 2013 + + Parameters + ---------- + x : float or array_like + Random number between 0 and 1. + + Returns + ------- + e : float or array_like + eccentricity + """ + e = np.sqrt(x) + + return e + + def random_keplarian_parameters(self, x, y, z, rng=None): + """ + Generate random inclination and angles of a binary system. + + Inclination uses the inverse CDF of isotropic orientations + (``i = arccos(s * x)`` with a random sign ``s``). ``Omega`` + and ``omega`` are uniform on [0, 360) deg from ``y`` and ``z``. + + Parameters + ---------- + x : float or array_like + Random number between 0 and 1, used for inclination. + y : float or array_like + Random number between 0 and 1, used for Omega. + z : float or array_like + Random number between 0 and 1, used for omega. + rng : numpy.random.Generator, optional + Random number generator used for the inclination sign. + Default is a new ``numpy.random.default_rng()``. + + Returns + ------- + inclination : float or array_like + Inclination angle in degrees. + Omega : float or array_like + Longitude of the ascending node in degrees. + omega : float or array_like + Argument of periastron in degrees. + """ + if rng is None: + rng = np.random.default_rng() + sign = rng.choice([-1, 1], size=len(x)) + x = sign * x + inclination = np.arccos(x) * 180 / np.pi + + Omega = 360 * y + omega = 360 * z + + return inclination, Omega, omega + class MultiplicityUnresolved(object): """ - The properties of stellar companions (see notes below). - The default parameters are as described in + SPISEA v2.5 default unresolved multiplicity (companions, no orbits). + + The default parameters are as described in `Lu et al. 2013 `_. These parameters are most appropriate for stellar populations - with ages <10 Myr. + with ages < 10 Myr. This is the unresolved default for backwards + compatibility. For a scientifically preferred model that includes + brown dwarfs, see :class:`MultiplicityUnresolvedOffner2023` + (opt-in; not the default). Notes ----- - The number of stellar companions, their masses, and separations - are be described by the following functions: + The number of stellar companions and their masses are described + by the following functions. - **Multiplicity Fraction** -- the number of stellar systems that host + **Multiplicity Fraction** -- the number of stellar systems that host multiple stars. In other words, the number of primary stars with companions. The multiplicity fraction (MF) is typically described as:: + B + T + Q + ... MF = --------------------- S + B + T + Q + ... where S = single, B = binary, T = triple, Q = quadruple, etc. - The MF also changes with mass and this dependency can be + The MF also changes with mass and this dependency can be described as a power-law:: - + MF(mass) = MF_amp * (mass ** MF_power) - However, in the brown dwarf mass regime, it is currently recognized - that only binaries are possible, and the MF decreases dissimilarly - to higher masses (> 0.08 solar masses). The values for this range - are given by Aberasturi et al. (2014) and Fontanive et al. (2023). + Defaults are MF_amp = 0.44, MF_power = 0.51 (Lu et al. 2013). + MF is clipped to [0, 1]. + + In the brown-dwarf regime, only binaries are expected and the MF + does not follow the stellar power law. For **scalar** mass, + :meth:`multiplicity_fraction` applies a staircase from + Aberasturi et al. (2014) and Fontanive et al. (2018):: + + M < 0.02 Msun MF = 0 + 0.02 < M <= 0.06 Msun MF = 0.08 + 0.06 < M <= 0.08 Msun MF = 0.16 + + Array masses use the stellar power law only (no staircase). + Cluster generation on this class therefore still uses + ``0.44 M**0.51`` for brown-dwarf primaries. **Companion Star Fraction** -- the expected number of companions in - a multiple system. The companion star fraction (CSF) also + a multiple system. The companion star fraction (CSF) also changes with mass and this dependency can be described as a power-law:: - + CSF(mass) = CSF_amp * (mass ** CSF_power) - The companion star fraction is clipped to some maximum - value, CSF_max. The actual number of companions is drawn + Defaults are CSF_amp = 0.50, CSF_power = 0.45. CSF is clipped + to CSF_max (default 3). The actual number of companions is drawn from a Poisson distribution with an expectation value of CSF. - - In the brown dwarf regime we impose an assumption that only - binary systems are possible due to current literature trends. + Higher-order multiples (triples+) are allowed at all masses, + including brown dwarfs. If ``companion_max`` is True, the + companion count is capped at CSF_max at all masses. **Mass Ratio (Q)** -- The ratio between the companion star - mass and primary star mass, Q = (m_comp / m_prim ) has - a probability density function described by a powerlaw:: + mass and primary star mass, Q = (m_comp / m_prim) has + a probability density function described by a power law:: P(Q) = Q ** q_power for q_min <= Q <= 1 - Current observations show no significant mass dependence. - + Stellar primaries use q_power = −0.4 (Lu et al. 2013). + Brown-dwarf primaries (M <= 0.08 Msun) use γ = 6.1 from + Fontanive et al. (2018) when ``mass`` is passed to + :meth:`q_power_at_mass` or :meth:`random_q`. That 0.08 Msun + switch is a mass-ratio index only, not a companion-count + policy. ``random_q(x)`` with no mass keeps the stellar-only + power law. + Parameters ---------- MF_amp : float, optional - The amplitude of the power-law describing the Multiplicity - Fraction as a function of stellar mass. - + Amplitude of the MF power law, dimensionless + (units of MF / Msun**MF_power). Default 0.44. MF_power : float, optional - The power of the power-law describing the Multiplicity - Fraction as a function of stellar mass. - + Power-law index of MF(M), dimensionless. Default 0.51. CSF_amp : float, optional - The amplitude of the power-law describing the companion star - fraction as a function of stellar mass. - + Amplitude of the CSF power law, dimensionless mean + companion count / Msun**CSF_power. Default 0.50. CSF_power : float, optional - The power of the power-law describing the companion star - fraction as a function of stellar mass. - + Power-law index of CSF(M), dimensionless. Default 0.45. CSF_max : float, optional - The maximum allowed companion star fraction, which is the - expectation value for the number of companion stars. Given - a CSF_max = 3, some systems will still have more than 3 - companions. - + Maximum companion star fraction, dimensionless mean + companion count (not bounded by 1). Default 3. q_power : float, optional - The power of the power-law describing the probability - density function for the mass ratio. - + Stellar mass-ratio power-law index γ, dimensionless. + Default −0.4 (Lu et al. 2013). q_min : float, optional - The minimum allowed Q value for the probability - density function of the mass ratio. - + Minimum mass ratio m_comp/m_prim, dimensionless. + Default 0.01. companion_max : bool, optional - Sets CSF_max is the max as the max number of companions. + If True, cap companion counts at CSF_max at all masses. Default False. - + bd_q_power : float, optional + Mass-ratio power-law index γ for brown-dwarf primaries + (M <= 0.08 Msun), dimensionless. Default 6.1 + (Fontanive et al. 2018). """ - def __init__(self, + def __init__(self, MF_amp=0.44, MF_power=0.51, CSF_amp=0.50, CSF_power=0.45, CSF_max=3, - q_power=-0.4, q_min=0.01, companion_max = False): - + q_power=-0.4, q_min=0.01, companion_max=False, + bd_q_power=6.1): + self.MF_amp = MF_amp self.MF_pow = MF_power self.CSF_amp = CSF_amp @@ -119,25 +200,29 @@ def __init__(self, self.q_pow = q_power self.q_min = q_min self.companion_max = companion_max + self.bd_q_power = bd_q_power def multiplicity_fraction(self, mass): """ Given a star's mass, determine the probability that the star is in a multiple system (multiplicity fraction = MF). - Modified to allow binary fraction to decrease in brown dwarf regime. - Supported by Aberasturi et al. (2014) and Fontanive et al. (2018). + Arrays use ``MF = MF_amp * M**MF_power`` clipped to 1. + Scalars also apply the Aberasturi/Fontanive brown-dwarf + staircase (0 / 8% / 16% below 0.02 / 0.06 / 0.08 Msun). Parameters ---------- - mass : float or numpy array - Mass of primary star. + mass : float or array_like + Primary mass n solar masses (Msun). Returns ------- - mf : float or numpy array - Multiplicity Fraction, the fraction of stars at this mass - that will have one or more companions. + mf : float or ndarray + Multiplicity fraction, dimensionless, in [0, 1]. + The fraction of stars at this mass that will have one or + more companions. Python float if ``mass`` is scalar, + ndarray otherwise. """ # Multiplicity Fraction mf = self.MF_amp * mass ** self.MF_pow @@ -160,57 +245,91 @@ def multiplicity_fraction(self, mass): def companion_star_fraction(self, mass): """ Given a star's mass, determine the average number of - companion stars (companion star fraction = CSF). For - brown dwarfs we impose a hard limit of one companion. + companion stars (companion star fraction = CSF). Parameters ---------- - mass : float or numpy array - Mass of primary star + mass : float or array_like + Primary mass in solar masses (Msun). Returns ------- - csf : float or numpy array - Companion Star Fraction, the expected number of companions - for a star at this mass. + csf : float or ndarray + Companion star fraction, the expected number of companions + for a star at this mass. Dimensionless mean companion count + (not bounded by 1). Python float if ``mass`` is scalar, + ndarray otherwise. """ # Companion Star Fraction csf = self.CSF_amp * mass ** self.CSF_pow - + if np.isscalar(csf): if csf > self.CSF_max: csf = self.CSF_max - if (mass <= 0.08): - csf = self.multiplicity_fraction(mass) else: csf[csf > self.CSF_max] = self.CSF_max - bd = mass <= 0.08 - csf[bd] = self.multiplicity_fraction(mass[bd]) return csf - def random_q(self, x): + def q_power_at_mass(self, mass): + """ + Mass-ratio power-law index, P(q) ∝ q ** q_power. + + Lu et al. (2013) use a single ``q_power`` for stellar primaries. + Brown-dwarf primaries (M <= 0.08 Msun) use ``bd_q_power`` + (default 6.1, Fontanive et al. 2018), matching the + companion-mass draw previously special-cased in + ``imf.calc_multi``. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + q_power : float or ndarray + Mass-ratio power-law index γ, dimensionless. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + q_pow = np.full(mass_arr.shape, self.q_pow, dtype=float) + q_pow[mass_arr <= H_BURNING_MASS] = self.bd_q_power + if np.isscalar(mass): + return float(q_pow[0]) + return q_pow + + def random_q(self, x, mass=None): """ Generative function for companion mass ratio, equivalent to the inverse of the CDF. - `q = m_compnaion / m_primary` + `q = m_companion / m_primary` `P(q) = q ** beta` for q_min <= q <= 1 Parameters ---------- x : float or array_like - Random number between 0 and 1. + Uniform random draw, dimensionless, in [0, 1]. Inverse CDF + sample for q. + + mass : float or array_like, optional + Primary mass in solar masses (Msun). If given, the + power-law index is ``q_power_at_mass(mass)`` (brown-dwarf + vs stellar for the SPISEA v2.5 default; mass-dependent for + Offner et al. 2023). If omitted, ``self.q_pow`` is used + for all companions. Returns ------- - q : float or array_like - companion mass ratio(s) + q : float or ndarray + Companion mass ratio m_comp/m_prim, dimensionless, in + [q_min, 1]. Python float if ``x`` is scalar, ndarray + otherwise. """ - b = 1.0 + self.q_pow - q = (x * (1.0 - self.q_min ** b) + self.q_min ** b) ** (1.0 / b) - - return q + if mass is None: + return _q_from_powerlaw(x, self.q_pow, self.q_min) + return _q_from_powerlaw(x, self.q_power_at_mass(mass), self.q_min) def random_is_multiple(self, x, MF): """ @@ -218,48 +337,253 @@ def random_is_multiple(self, x, MF): """ return x < MF - def random_companion_count(self, x, CSF, MF): + def random_companion_count(self, x, CSF, MF, mass=None, rng=None): """ - Helper function: calculate number of companions. - """ - # bd stipulation since mf=0 - if MF <= 0: - return 0 + Number of companions for primaries already identified as multiple. - n_comp = 1 + np.random.poisson((CSF / MF) - 1) - - if self.companion_max == True: - if n_comp > self.CSF_max: + The count is drawn from a Poisson with expectation CSF/MF - 1, + then 1 is added so every multiple has at least one companion. + ``x`` is unused and kept for API compatibility. + + Parameters + ---------- + x : float or array_like + Unused (historical signature). Dimensionless uniform + draw in [0, 1] if provided. + CSF : float or array_like + Companion star fraction, dimensionless mean companion + count (not bounded by 1). + MF : float or array_like + Multiplicity fraction, dimensionless, in [0, 1]. + mass : float or array_like, optional + Primary mass in solar masses (Msun). + rng : numpy.random.Generator, optional + Random generator. If omitted, uses ``numpy.random`` (the + historical scalar helper). + + Returns + ------- + n_comp : int or ndarray of int + Number of companions, integer count. Python int if ``CSF`` + and ``MF`` are scalar, ndarray otherwise. + """ + return_scalar = np.isscalar(CSF) and np.isscalar(MF) + if return_scalar and rng is None: + if MF <= 0: + return 0 + n_comp = 1 + np.random.poisson((CSF / MF) - 1) + if self.companion_max and n_comp > self.CSF_max: n_comp = self.CSF_max - + return int(n_comp) + + CSF = np.atleast_1d(np.asarray(CSF, dtype=float)) + MF = np.atleast_1d(np.asarray(MF, dtype=float)) + if rng is None: + n_comp = 1 + np.random.poisson((CSF / MF) - 1) + else: + n_comp = 1 + rng.poisson((CSF / MF) - 1) + + if self.companion_max: + n_comp = np.minimum(n_comp, self.CSF_max) + + if return_scalar: + return int(n_comp[0]) + return n_comp + + def draw_n_companions(self, mass, CSF, MF, rng): + """ + Vectorized companion counts for primaries that are already + identified as multiple. Delegates to + :meth:`random_companion_count`. If ``companion_max`` is True, + counts are capped at CSF_max at all masses. + + Parameters + ---------- + mass : array_like + Primary masses of systems already identified as multiple. + Must be positive, in solar masses (Msun). + CSF : array_like + Companion star fraction at each primary, dimensionless + mean companion count (not bounded by 1). + MF : array_like + Multiplicity fraction at each primary, dimensionless, + in [0, 1]. + rng : numpy.random.Generator + Random generator used for the Poisson companion-count draw. + + Returns + ------- + n_comp : ndarray of int + Number of companions per primary, integer count, shape + matching ``mass``. + """ + n_comp = self.random_companion_count(None, CSF, MF, mass=mass, rng=rng) + + return np.atleast_1d(n_comp) + + def _q_values_for_primaries(self, prim_subset, n_comp, rng): + """ + Draw mass ratios for ``n_comp`` companions of each primary. + + The stellar / brown-dwarf split (two separate RNG draws) preserves + the historical SPISEA v2.5 random sequence used by + ``imf.calc_multi``. + + Parameters + ---------- + prim_subset : array_like + Primary masses for this companion-count group. Must be + positive, in solar masses (Msun). + n_comp : int + Number of companions per primary, integer count. + rng : numpy.random.Generator + Random generator used for inverse-CDF q draws. + + Returns + ------- + q_values : ndarray + Companion mass ratios m_comp/m_prim, dimensionless, in + [q_min, 1]. Shape (len(prim_subset), n_comp). + """ + q_values = np.empty((len(prim_subset), n_comp)) + bd_mask = prim_subset <= H_BURNING_MASS + star_mask = ~bd_mask + + if np.any(star_mask): + q_values[star_mask] = self.random_q(rng.random((star_mask.sum(), n_comp)), mass=prim_subset[star_mask]) + + if np.any(bd_mask): + q_values[bd_mask] = self.random_q(rng.random((bd_mask.sum(), n_comp)), mass=prim_subset[bd_mask]) + + return q_values + + def draw_companion_masses(self, primary_masses, is_multiple, CSF, MF, + rng, mass_min): + """ + Assign companion masses for a set of primaries. + + This is the multiplicity-object entry point used by + ``IMF.calc_multi``. Companion-mass draws, including brown-dwarf + q distributions and the binaries-only BD companion count, live + here rather than in ``imf.py``. + + Parameters + ---------- + primary_masses : array_like + Primary masses must be positive, in solar masses (Msun). + is_multiple : array_like of bool + True for primaries drawn as multiple systems. + CSF : array_like + Companion star fraction at each primary, dimensionless + mean companion count (not bounded by 1). + MF : array_like + Multiplicity fraction at each primary, dimensionless, + in [0, 1]. + rng : numpy.random.Generator + Random generator. + mass_min : float + Minimum companion mass in solar masses (Msun); lighter + companions are masked. + + Returns + ------- + comp_masses : numpy.ma.MaskedArray + Companion masses in solar masses (Msun), shape + (n_primaries, max_n_comp). + system_masses : ndarray + Primary plus unmasked companion mass, in solar masses + (Msun). + is_multiple : ndarray of bool + Updated multiplicity flags after masking sub-minimum + companions. + """ + primary_masses = np.asarray(primary_masses, dtype=float) + is_multiple = np.asarray(is_multiple, dtype=bool) + CSF = np.asarray(CSF, dtype=float) + MF = np.asarray(MF, dtype=float) + + system_masses = primary_masses.copy() + multiple_idx = np.where(is_multiple)[0] + primary = primary_masses[multiple_idx] + n_comp = self.draw_n_companions(primary, CSF[multiple_idx], MF[multiple_idx], rng) + + if len(multiple_idx) == 0: + comp_masses = np.zeros((len(primary_masses), 1), dtype=float) + comp_masses = np.ma.MaskedArray(comp_masses, mask=comp_masses < mass_min) + + return comp_masses, system_masses, is_multiple + + n_unique = np.unique(n_comp) + n_indices = [np.where(n_comp == i)[0] for i in n_unique] + comp_masses = np.zeros((len(primary_masses), int(np.max(n_unique)))) + + for n_c, idx in zip(n_unique, n_indices): + prim_subset = primary[idx] + q_values = self._q_values_for_primaries(prim_subset, int(n_c), rng) + m_comp = q_values * prim_subset[:, np.newaxis] + comp_masses[multiple_idx[idx], :int(n_c)] = m_comp + + comp_masses = np.ma.MaskedArray(comp_masses, mask=comp_masses < mass_min) + system_masses[multiple_idx] += comp_masses[multiple_idx].sum(axis=1) + is_multiple = np.any(~comp_masses.mask, axis=1) + + return comp_masses, system_masses, is_multiple -class MultiplicityResolvedDK(MultiplicityUnresolved): + +class MultiplicityResolvedDK(MultiplicityUnresolved, _ResolvedOrbitalMixin): """ - Sub-class of MultiplicityUnresolved that adds semimajor axis and eccentricity information - for multiple objects from distributions described in Duchene and Kraus 2013 + SPISEA v2.5 default resolved multiplicity. + + Same MF/CSF/q as :class:`MultiplicityUnresolved` (Lu et al. 2013 + stellar power laws, scalar BD staircase, Fontanive γ = 6.1), plus + semimajor axis and eccentricity from Duchêne & Kraus (2013). + This is the resolved default for backwards compatibility. It is + **not** meant for the brown-dwarf regime the way + :class:`MultiplicityResolvedOffner2023` is. - For brown dwarf regime, mean separation and std are given by Fontanive et al. (2018). + Notes + ----- + Characteristic μ(a) is a broken power law in a vs M + (``astropy.modeling.powerlaws.BrokenPowerLaw1D`` with the fitted + ``a_amp``, ``a_break``, ``a_slope1``, ``a_slope2``). + σ(log10 a) is a linear fit in log M (``a_std_slope``, + ``a_std_intercept``) that saturates above 2.9 Msun and is + clipped to ≥ 0.1 dex. The dip near 0.08 Msun is the BD/stellar + blend: Fontanive et al. (2018) interpolation of the BD mean + (2.5–8 AU) and width (0.25–0.5 dex) over 0.01–0.08 Msun, + combined with the stellar law by a sigmoid at 0.08 Msun. + + :meth:`log_a_mean`, :meth:`a_mean`, and :meth:`sigma_log_a` + return that same characteristic mean and width. + :meth:`log_semimajoraxis` draws from a truncated lognormal + with ``loc = log_a_mean(mass)`` and + ``scale = sigma_log_a(mass)``, truncated to 0.01–2000 AU. + Eccentricity follows f(e) = 2e; inclination and angles are + random (Duchêne & Kraus 2013). Parameters - -------------- - a_amp: float, optional - Ampltiude of the broken power law describing the log_semimajoraxis - - a_break: float, optional - Break location on the x-axis of the broken power law describing the log_semimajoraxis - - a_slope1: float, optional - Slope of the left side of the broken power law describing the log_semimajoraxis - - a_slope2: float, optional - Slope of the right side of the broken power law describing the log_semimajoraxis - - a_std_slope: float, optional - Slope of the line that fit sigma_log_semimajoraxis vs log_mass - - a_std_intercept: float, optional - Intercept of the line that fit sigma_log_semimajoraxis vs log_mass + ---------- + a_amp : float, optional + Amplitude of the broken power law for μ(a), in AU. + Default 379.79953034. + a_break : float, optional + Break mass of the broken power law, in solar masses (Msun). + Default 4.90441533. + a_slope1 : float, optional + Power-law index below ``a_break``, dimensionless. + Default −1.80171539. + a_slope2 : float, optional + Power-law index above ``a_break``, dimensionless. + Default 4.23325571. + a_std_slope : float, optional + Slope of σ(log10 a) vs log10(M / Msun), in dex per dex. + Default 1.19713084. + a_std_intercept : float, optional + Intercept of σ(log10 a) vs log10(M / Msun), in dex. + Default 1.28974264. + **kwargs + Passed to :class:`MultiplicityUnresolved` (MF/CSF/q). """ def __init__(self, a_amp = 379.79953034, a_break = 4.90441533, a_slope1 = -1.80171539, a_slope2 = 4.23325571, a_std_slope = 1.19713084, a_std_intercept = 1.28974264, **kwargs): @@ -270,113 +594,1142 @@ def __init__(self, a_amp = 379.79953034, a_break = 4.90441533, a_slope1 = -1.801 self.a_slope2 = a_slope2 self.a_std_slope = a_std_slope self.a_std_intercept = a_std_intercept - - def log_semimajoraxis(self, mass): + + return + + def log_a_mean(self, mass): """ - Generate the semimajor axis for a given mass. The mean and standard deviation of a given mass are determined - by fitting the data from fitting the semimajor axis data as a function of mass in table 1 of Duchene and Kraus 2013. - Then a random semimajor axis is drawn from a log normal distribution with that mean and standard deviation. + Characteristic log10(a/AU) used as the loc of + :meth:`log_semimajoraxis`. - The brown dwarf range is covered by mass-dependent scaling of both the characteristic separation and dispersion - matching trends described in Fontanive et al. (2018). + Duchêne & Kraus (2013) broken power law for stellar primaries, + Fontanive et al. (2018) interpolation for brown dwarfs, and a + sigmoid blend at 0.08 Msun. Parameters ---------- - mass : array-like - Mass array of primary star + mass : float or array_like + Primary mass must be positive, in solar masses (Msun). Returns ------- - log_semimajoraxis : array-like - Log of the semimajor axis/separation between the stars in units of AU + log_a_mean : float or ndarray + Characteristic log10(a / 1 AU) in dex (not ln, not AU). + Python float if ``mass`` is scalar, ndarray otherwise. """ - mass = np.atleast_1d(mass) + return_scalar = np.isscalar(mass) + mass = np.atleast_1d(np.asarray(mass, dtype=float)) logm = np.log10(mass) - # Stellar mean and std (Duchene & Kraus 2013) - a_mean_func = astropy.modeling.powerlaws.BrokenPowerLaw1D(amplitude=self.a_amp, x_break=self.a_break, - alpha_1=self.a_slope1, alpha_2=self.a_slope2) - log_a_mean_star = np.log10(a_mean_func(mass)) # mean log(a) - log_a_std_func = astropy.modeling.models.Linear1D(slope=self.a_std_slope, intercept=self.a_std_intercept) - log_a_std_star = log_a_std_func(logm) # sigma_log(a) - log_a_std_star[mass >= 2.9] = log_a_std_func(np.log10(2.9)) # sigma_log(a) - log_a_std_star = np.clip(log_a_std_star, 0.1, None) - - # BD mean and std (Fontanive+18): interpolated over substellar range + a_mean_func = astropy.modeling.powerlaws.BrokenPowerLaw1D( + amplitude=self.a_amp, x_break=self.a_break, + alpha_1=self.a_slope1, alpha_2=self.a_slope2) + log_a_mean_star = np.log10(a_mean_func(mass)) log_a_mean_bd = np.interp( logm, [np.log10(0.01), np.log10(0.08)], [np.log10(2.5), np.log10(8.0)] ) + w = 1.0 / (1.0 + np.exp(-(logm - np.log10(0.08)) / 0.15)) + log_a_mean = (1 - w) * log_a_mean_bd + w * log_a_mean_star + if return_scalar: + return float(log_a_mean[0]) + return log_a_mean + + def a_mean(self, mass): + """ + Characteristic μ(a) in AU, ``10 ** log_a_mean(mass)``. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, in solar masses (Msun). + + Returns + ------- + a_mean : float or ndarray + Characteristic separation μ(a) in AU. Python float if + ``mass`` is scalar, ndarray otherwise. + """ + log_a = self.log_a_mean(mass) + if np.isscalar(log_a): + return 10.0 ** log_a + return 10.0 ** np.asarray(log_a, dtype=float) + + def sigma_log_a(self, mass): + """ + σ(log10 a) used as the scale of :meth:`log_semimajoraxis`. + + Linear fit in log-mass for stellar primaries (saturates above + 2.9 Msun, clipped to ≥ 0.1), Fontanive et al. (2018) + interpolation for brown dwarfs, and a sigmoid blend at + 0.08 Msun. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, in solar masses (Msun). + + Returns + ------- + sigma_log_a : float or ndarray + Standard deviation of log10(a / 1 AU), in dex. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + return_scalar = np.isscalar(mass) + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + logm = np.log10(mass) + + log_a_std_func = astropy.modeling.models.Linear1D( + slope=self.a_std_slope, intercept=self.a_std_intercept) + log_a_std_star = log_a_std_func(logm) + log_a_std_star[mass >= 2.9] = log_a_std_func(np.log10(2.9)) + log_a_std_star = np.clip(log_a_std_star, 0.1, None) log_a_std_bd = np.interp( logm, [np.log10(0.01), np.log10(0.08)], [0.25, 0.5] ) - - # Sigmoid blend: smoothly transitions from BD to stellar regime at 0.08 M_sun w = 1.0 / (1.0 + np.exp(-(logm - np.log10(0.08)) / 0.15)) - log_a_mean = (1 - w) * log_a_mean_bd + w * log_a_mean_star log_a_std = (1 - w) * log_a_std_bd + w * log_a_std_star + if return_scalar: + return float(log_a_std[0]) + return log_a_std + + def log_semimajoraxis(self, mass, rng=None): + """ + Draw log10(a/AU) from a mass-dependent truncated lognormal. + + The mean and standard deviation at a given mass come from + fitting Table 1 of Duchêne and Kraus 2013. The brown dwarf + range uses mass-dependent scaling of both the characteristic + separation and dispersion (Fontanive et al. 2018). + + Draws use ``loc = log_a_mean(mass)`` and + ``scale = sigma_log_a(mass)``. + + Parameters + ---------- + mass : array-like + Primary mass must be positive, in solar masses (Msun). + rng : numpy.random.Generator, optional + Random number generator passed to ``truncnorm.rvs``. + Default is a new ``numpy.random.default_rng()``. + + Returns + ------- + log_semimajoraxis : array-like + Drawn log10(a / 1 AU) in dex, truncated so a is between + 0.01 AU and 2000 AU. + """ + if rng is None: + rng = np.random.default_rng() + mass = np.atleast_1d(mass) + log_a_mean = np.atleast_1d(np.asarray(self.log_a_mean(mass), + dtype=float)) + log_a_std = np.atleast_1d(np.asarray(self.sigma_log_a(mass), + dtype=float)) - # Trunc normal distribution between log10(0.01) AU and log10(2000) AU + # Trunc normal between log10(0.01) AU and log10(2000) AU log_a_lower = np.log10(0.01) log_a_upper = np.log10(2000) a_lower_std = (log_a_lower - log_a_mean) / log_a_std a_upper_std = (log_a_upper - log_a_mean) / log_a_std - log_semimajoraxis = truncnorm.rvs(a_lower_std, a_upper_std, loc=log_a_mean, scale=log_a_std) + log_semimajoraxis = truncnorm.rvs( + a_lower_std, a_upper_std, loc=log_a_mean, scale=log_a_std, + random_state=rng) return log_semimajoraxis - def random_e(self, x): + +class MultiplicityPiecewisePowerLaw(MultiplicityUnresolved): + """ + Generic helper: multiplicity as a piecewise power law in primary mass. + + Use this when MF and CSF should change slope at user-specified + mass edges (survey-specific segments), rather than a single + power law or a logistic. Offner et al. 2023 does **not** use + this class; it uses :class:`MultiplicityLogistic`. + + Notes + ----- + On each mass segment i, with edges + ``mass_limits[i] <= M < mass_limits[i+1]``:: + + MF(M) = MF_amp[i] * M ** MF_power[i] + CSF(M) = CSF_amp[i] * M ** CSF_power[i] + + MF is clipped to [0, 1]. CSF is clipped to [0, CSF_max] and + raised to at least MF so the Poisson companion-count draw is + well defined. Higher-order multiples are allowed at all masses. + If ``companion_max`` is True, counts are capped at CSF_max at + all masses. Mass-ratio draws use a single ``q_power`` (same as + :class:`MultiplicityUnresolved`) unless a subclass overrides + :meth:`q_power_at_mass`. There is no scalar-only brown-dwarf + staircase; each segment is evaluated for both scalar and array + mass. + + Parameters + ---------- + mass_limits : array_like + Segment edges in solar masses (Msun), length N+1, strictly + increasing. + MF_amps : array_like + Length-N amplitudes for the multiplicity fraction, + dimensionless (units of MF / Msun**MF_power). + MF_powers : array_like + Length-N power-law indices for the multiplicity fraction, + dimensionless. + CSF_amps : array_like + Length-N amplitudes for the companion star fraction, + dimensionless (mean companion count / Msun**CSF_power). + CSF_powers : array_like + Length-N power-law indices for the companion star fraction, + dimensionless. + CSF_max : float, optional + Maximum companion star fraction, dimensionless mean companion + count (not bounded by 1). Default 3. + q_power : float, optional + Mass-ratio power-law index, dimensionless. Default -0.4. + q_min : float, optional + Minimum mass ratio m_comp/m_prim, dimensionless. Default 0.01. + companion_max : bool, optional + If True, cap companion counts at CSF_max at all masses. + Default False. + """ + def __init__(self, mass_limits, MF_amps, MF_powers, CSF_amps, CSF_powers, + CSF_max=3, q_power=-0.4, q_min=0.01, companion_max=False): + mass_limits = np.asarray(mass_limits, dtype=float) + MF_amps = np.asarray(MF_amps, dtype=float) + MF_powers = np.asarray(MF_powers, dtype=float) + CSF_amps = np.asarray(CSF_amps, dtype=float) + CSF_powers = np.asarray(CSF_powers, dtype=float) + + nseg = len(MF_amps) + + if len(mass_limits) != nseg + 1: + raise ValueError('len(mass_limits) must be len(MF_amps) + 1') + + if not (len(MF_powers) == len(CSF_amps) == len(CSF_powers) == nseg): + raise ValueError('MF/CSF amplitude and power arrays must have equal length') + + if np.any(np.diff(mass_limits) <= 0): + raise ValueError('mass_limits must be strictly increasing') + + super(MultiplicityPiecewisePowerLaw, self).__init__( + MF_amp=MF_amps[-1], MF_power=MF_powers[-1], + CSF_amp=CSF_amps[-1], CSF_power=CSF_powers[-1], + CSF_max=CSF_max, q_power=q_power, q_min=q_min, + companion_max=companion_max) + + self.mass_limits = mass_limits + self.MF_amps = MF_amps + self.MF_powers = MF_powers + self.CSF_amps = CSF_amps + self.CSF_powers = CSF_powers + + return + + def multiplicity_fraction(self, mass): """ - Generate random eccentricity from the inverse of the CDF where the PDF is f(e) = 2e from Duchene and Kraus 2013 - + Multiplicity fraction as a piecewise power law in primary mass. + Clipped to [0, 1]. + Parameters ---------- - x : float or array_like - Random number between 0 and 1. + mass : float or array_like + Primary mass in solar masses (Msun). Returns ------- - e : float or array_like - companion mass ratio(s) + mf : float or ndarray + Multiplicity fraction, dimensionless, in [0, 1]. + Python float if ``mass`` is scalar, ndarray otherwise. """ - e = np.sqrt(x) - - return e - - def random_keplarian_parameters(self, x, y, z): + mf = _piecewise_powerlaw(mass, self.mass_limits, self.MF_amps, self.MF_powers, clip_min=0.0, clip_max=1.0) + + return mf + + def companion_star_fraction(self, mass): """ - Generate random incliniation and angles of binary system - + Companion star fraction as a piecewise power law in primary mass. + + Clipped to [0, CSF_max] and raised to at least MF. + Parameters ---------- - x : float or array_like - Random number between 0 and 1. - - y : float or array_like - Random number between 0 and 1. - - z : float or array_like - Random number between 0 and 1. + mass : float or array_like + Primary mass in solar masses (Msun). Returns ------- - inclination : float or array_like - Angle of inclination - - Omega : float or array_like - Big Omega: one other angle of the system - - omega : float or array_like - Final angle of the system + csf : float or ndarray + Companion star fraction, dimensionless mean companion + count (not bounded by 1). Python float if ``mass`` is + scalar, ndarray otherwise. """ - sign = np.array([choice([-1,1]) for i in range(len(x))]) - x = sign*x - inclination = np.arccos(x)*180/np.pi #inclination angle in degrees + return_scalar = np.isscalar(mass) + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + csf = _piecewise_powerlaw( + mass_arr, self.mass_limits, self.CSF_amps, self.CSF_powers, + clip_min=0.0, clip_max=self.CSF_max) + mf = _piecewise_powerlaw( + mass_arr, self.mass_limits, self.MF_amps, self.MF_powers, + clip_min=0.0, clip_max=1.0) + csf = np.maximum(csf, mf) + if return_scalar: + return float(csf[0]) + return csf + + +class MultiplicityLogistic(MultiplicityUnresolved): + """ + Generic helper: multiplicity as a logistic in log primary mass. + + Use this for a C-infinity smooth MF/CSF that saturates at low + and high mass. :class:`MultiplicityUnresolvedOffner2023` is this + class with coefficients fitted to Offner et al. (2023) Table 1. + + Notes + ----- + The same 4-parameter logistic is used for MF and CSF with + independent coefficients:: + + f(M) = A + (B - A) / (1 + (M / M0)**(-k)) + + As M → 0, f → A; as M → ∞, f → B. The curve is not a piecewise + interpolation of survey knots. MF is clipped to [0, 1]. CSF is + clipped to [0, CSF_max] and raised to at least MF. Higher-order + multiples are allowed at all masses. If ``companion_max`` is + True, counts are capped at CSF_max at all masses. Mass-ratio + draws use a single ``q_power`` unless a subclass overrides + :meth:`q_power_at_mass`. + + Parameters + ---------- + MF_A, MF_B : float + Low-mass and high-mass asymptotes of the multiplicity-fraction + logistic, dimensionless (MF). + MF_M0 : float + Characteristic mass of the MF logistic, in solar masses (Msun). + MF_k : float + MF logistic slope, dimensionless. + CSF_A, CSF_B : float + Low-mass and high-mass asymptotes of the companion-star-fraction + logistic, dimensionless mean companion count (not bounded by 1). + CSF_M0 : float + Characteristic mass of the CSF logistic, in solar masses (Msun). + CSF_k : float + CSF logistic slope, dimensionless. + CSF_max : float, optional + Maximum companion star fraction, dimensionless mean companion + count (not bounded by 1). Default 3. + q_power : float, optional + Mass-ratio power-law index, dimensionless. Default -0.4. + q_min : float, optional + Minimum mass ratio m_comp/m_prim, dimensionless. Default 0.01. + companion_max : bool, optional + If True, cap companion counts at CSF_max at all masses. + Default False. + """ + def __init__(self, MF_A, MF_B, MF_M0, MF_k, + CSF_A, CSF_B, CSF_M0, CSF_k, + CSF_max=3, q_power=-0.4, q_min=0.01, companion_max=False): + super(MultiplicityLogistic, self).__init__( + MF_amp=1.0, MF_power=0.0, CSF_amp=1.0, CSF_power=0.0, + CSF_max=CSF_max, q_power=q_power, q_min=q_min, + companion_max=companion_max) + self.MF_A = float(MF_A) + self.MF_B = float(MF_B) + self.MF_M0 = float(MF_M0) + self.MF_k = float(MF_k) + self.CSF_A = float(CSF_A) + self.CSF_B = float(CSF_B) + self.CSF_M0 = float(CSF_M0) + self.CSF_k = float(CSF_k) + + return + + def multiplicity_fraction(self, mass): + """ + Multiplicity fraction as a logistic in log primary mass. + Clipped to [0, 1]. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + mf : float or ndarray + Multiplicity fraction, dimensionless, in [0, 1]. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + mf = _logistic_in_logm(mass, self.MF_A, self.MF_B, self.MF_M0, self.MF_k, + clip_min=0.0, clip_max=1.0) + + return mf + + def companion_star_fraction(self, mass): + """ + Companion star fraction as a logistic in log primary mass. + + Clipped to [0, CSF_max] and raised to at least MF. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + csf : float or ndarray + Companion star fraction, dimensionless mean companion + count (not bounded by 1). Python float if ``mass`` is + scalar, ndarray otherwise. + """ + return_scalar = np.isscalar(mass) + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + + # Calculate the multiplicity fraction + mf = _logistic_in_logm(mass_arr, self.MF_A, self.MF_B, self.MF_M0, self.MF_k, + clip_min=0.0, clip_max=1.0) + + # Calculate the companion star fraction + csf = _logistic_in_logm(mass_arr, self.CSF_A, self.CSF_B, self.CSF_M0, self.CSF_k, + clip_min=0.0, clip_max=self.CSF_max) - Omega = 360*y - omega = 360*z + # Ensure the companion star fraction is at least the multiplicity fraction + csf = np.maximum(csf, mf) + + if return_scalar: + return float(csf[0]) + + return csf + + +class MultiplicityUnresolvedOffner2023(MultiplicityLogistic): + """ + Opt-in unresolved multiplicity derived from data in Offner et al. 2023 Table 1, + including brown dwarfs. + + Scientifically preferred over the SPISEA v2.5 default, but **not** + the default (backwards compatibility). Companions only; for orbits + use :class:`MultiplicityResolvedOffner2023`. + + Citation: Offner, S. S. R., Moe, M., Kratter, K. M., Sadavoy, S. I., + Jensen, E. L. N., & Tobin, J. J. 2023, in Protostars and Planets VII, + ASP Conf. Ser. 534, 275 (`arXiv:2203.10066 + `_; ADS + `2023ASPC..534..275O + `_). + Table 1 data: Zenodo `10.5281/zenodo.6628915 + `_. + + Notes + ----- + MF and CSF are a **4-parameter logistic in log-mass**, fitted with + equal weight to the geom-mean MF/CF columns of Table 1:: + + f(M) = A + (B - A) / (1 + (M / M0)**(-k)) + + with (A, B, M0, k) = (0.14, 0.99, 1.41, 1.25) for MF and + (0.12, 2.35, 3.57, 0.96) for CSF/CF. The curve is C-infinity + smooth (not a broken power law and not interpolation of Table 1 + knots), saturates at B ~ 1 for MF so A/B stars stay near the + Raghavan/MDS/Sana points, and has a low-mass floor A ~ 0.14. + Fontanive et al. (2018) 8 ± 6% sits ~0.07 below the curve (~15%), + consistent with the Burgasser/Close BD points and within ~1–2σ + of Fontanive. MF is clipped to [0, 1]. CSF is clipped to + [0, CSF_max] and raised to at least MF. Higher-order multiples + are allowed at all masses, including brown dwarfs. If + ``companion_max`` is True, counts are capped at CSF_max at all + masses. + + **Companion assignment vs Table 1.** Offner et al. 2023 (text + above Table 1): BD primaries include all BD companions; FGKM MS + statistics include only MS companions with M_comp > 0.075 Msun; + OBA include MS companions above q > 0.1. Table 1 stellar MF/CF + therefore exclude BD companions. SPISEA still draws companions + down to ``q_min`` (default 0.01), so brown-dwarf secondaries of + stellar primaries are generated. The solar-type BD-companion + fraction is only ≈ 4% (BD desert at a < 0.5 au), so the + integrated stellar MF is affected very little. Do not interpret + the simulated stellar-primary MF as a stellar-companion-only + statistic. + + **Mass-ratio index.** :meth:`q_power_at_mass` is an error-weighted + logistic in log-mass fitted to Table 1 γ_trunc (1–100 au):: + + γ(M) = A + (B - A) / (1 + (M / M0)**(-k)) + + with (A, B, M0, k) = (6.6, −1.77, 0.0651, 0.629). Call + ``random_q(x, mass=...)``. Without ``mass``, ``random_q(x)`` + keeps the historical stellar-only power law. BD companions are + still more equal-mass than solar-type companions. The + err-weighted fit undershoots Fontanive 4.8 ± 2.2 (~3.3 at + 0.033 Msun); that is the fit, not a bug. + + **Characteristic separation.** :meth:`log_a_mean` / :meth:`a_mean` + are a smooth broken power law in log10 a vs log10 M (FGK-pulled, + s = 0.1 dex), C-infinity via a stable logcosh (not + ``log(cosh x)`` and not a hard ``where`` break):: + + v = log10(M / Mp), yp = log10(μp) + log10 a = yp + 0.5*(αL+αR)*v + 0.5*(αR-αL)*s * logcosh(v/s) + + with μp = 44.46 AU, Mp = 0.819 Msun, αL = 1.005, αR = −0.308, + s = 0.10. Linear-space a is clipped above 0.1 AU. The + implementation uses + ``logcosh x = |x| + log(1 + e**(-2|x|)) - log 2``. + + **Separation scatter.** :meth:`sigma_log_a` is a 2-parameter + logistic pinned at 0.7 / 1.5:: + + σ(M) = 0.7 + 0.8 / (1 + (M / 0.354)**(-6.05)) + + i.e. (A, B, M0, k) = (0.7, 1.5, 0.354, 6.05), clipped to ≥ 0.1. + Resolved draws use these as loc / scale of a truncated lognormal + (see :class:`MultiplicityResolvedOffner2023`). + + Parameters + ---------- + MF_A, MF_B : float, optional + Low- and high-mass MF logistic asymptotes, dimensionless. + Defaults 0.14, 0.99 (Offner et al. 2023 Table 1, + equal-weight geom-mean MF fit). + MF_M0 : float, optional + Characteristic mass of the MF logistic, in solar masses + (Msun). Default 1.41. + MF_k : float, optional + MF logistic slope, dimensionless. Default 1.25. + CSF_A, CSF_B : float, optional + Low- and high-mass CSF logistic asymptotes, dimensionless + mean companion count. Defaults 0.12, 2.35 (Table 1 + equal-weight geom-mean CF fit). + CSF_M0 : float, optional + Characteristic mass of the CSF logistic, in solar masses + (Msun). Default 3.57. + CSF_k : float, optional + CSF logistic slope, dimensionless. Default 0.96. + q_A, q_B : float, optional + Low- and high-mass γ logistic asymptotes, dimensionless. + Defaults 6.6, −1.77 (Table 1 γ_trunc, error-weighted). + q_M0 : float, optional + Characteristic mass of the γ logistic, in solar masses + (Msun). Default 0.0651. + q_k : float, optional + γ logistic slope, dimensionless. Default 0.629. + a_mup : float, optional + Smooth-broken-PL pivot μ(a), in AU. Default 44.46. + a_mp : float, optional + Smooth-broken-PL pivot mass, in solar masses (Msun). + Default 0.819. + a_alphaL, a_alphaR : float, optional + Smooth-broken-PL slopes in log10 a vs log10 M, + dimensionless. Defaults 1.005, −0.308. + a_s : float, optional + Smooth-broken-PL smoothing scale, in dex. Default 0.10. + a_min : float, optional + Minimum characteristic a, in AU. Default 0.1. + sig_A, sig_B : float, optional + Low- and high-mass σ(log10 a) logistic values, in dex. + Defaults 0.7, 1.5 (Table 2 pins). + sig_M0 : float, optional + Characteristic mass of the σ logistic, in solar masses + (Msun). Default 0.354. + sig_k : float, optional + σ logistic slope, dimensionless. Default 6.05. + CSF_max : float, optional + Maximum companion star fraction, dimensionless mean companion + count (not bounded by 1). Default 3. + q_power : float, optional + Fallback mass-ratio power-law index, dimensionless. Ignored for + draws when primary mass is provided (the γ logistic is used); + used by ``random_q(x)`` with no mass. Default 0.2. + q_min : float, optional + Minimum mass ratio m_comp/m_prim, dimensionless, in [q_min, 1]. + Default 0.01. + companion_max : bool, optional + If True, cap companion counts at CSF_max at all masses. + Default False. + """ + def __init__(self, MF_A=0.14, MF_B=0.99, MF_M0=1.41, MF_k=1.25, + CSF_A=0.12, CSF_B=2.35, CSF_M0=3.57, CSF_k=0.96, + q_A=6.6, q_B=-1.77, q_M0=0.0651, q_k=0.629, + a_mup=44.46, a_mp=0.819, a_alphaL=1.005, + a_alphaR=-0.308, a_s=0.10, a_min=0.1, + sig_A=0.7, sig_B=1.5, sig_M0=0.354, sig_k=6.05, + CSF_max=3, q_power=0.2, q_min=0.01, + companion_max=False): + + super(MultiplicityUnresolvedOffner2023, self).__init__( + MF_A=MF_A, MF_B=MF_B, MF_M0=MF_M0, MF_k=MF_k, + CSF_A=CSF_A, CSF_B=CSF_B, CSF_M0=CSF_M0, CSF_k=CSF_k, + CSF_max=CSF_max, q_power=q_power, q_min=q_min, + companion_max=companion_max) + self.q_A = float(q_A) + self.q_B = float(q_B) + self.q_M0 = float(q_M0) + self.q_k = float(q_k) + self.a_mup = float(a_mup) + self.a_mp = float(a_mp) + self.a_alphaL = float(a_alphaL) + self.a_alphaR = float(a_alphaR) + self.a_s = float(a_s) + self.a_min = float(a_min) + self.sig_A = float(sig_A) + self.sig_B = float(sig_B) + self.sig_M0 = float(sig_M0) + self.sig_k = float(sig_k) + + return + + def q_power_at_mass(self, mass): + """ + Mass-ratio power-law index γ(M), P(q) ∝ q^γ on [q_min, 1]. + + Error-weighted logistic in log-mass fitted to Table 1 + γ_trunc (1–100 au). Not an interpolation of the Table 1 + knots:: + + γ(M) = A + (B - A) / (1 + (M / M0)**(-k)) + + with (A, B, M0, k) = (6.6, −1.77, 0.0651, 0.629). + Undershoots Fontanive 4.8 ± 2.2 (~3.3 at 0.033 Msun). + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + gamma : float or ndarray + Mass-ratio power-law index γ, dimensionless. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + gamma = _logistic_in_logm( + mass, self.q_A, self.q_B, self.q_M0, self.q_k) - return inclination, Omega, omega + return gamma + + def log_a_mean(self, mass): + """ + Characteristic log10(a/AU) from the smooth broken power law. + + FGK-pulled, s = 0.1 dex, C-infinity (stable logcosh; not + ``log(cosh x)`` and not a hard ``where`` break):: + + v = log10(M / Mp), yp = log10(μp) + log10 a = yp + 0.5*(αL+αR)*v + 0.5*(αR-αL)*s * logcosh(v/s) + + with μp = 44.46 AU, Mp = 0.819 Msun, αL = 1.005, + αR = −0.308, s = 0.10. Linear-space a is clipped to 0.1 AU. + Uses ``logcosh x = |x| + log(1 + e**(-2|x|)) - log 2``. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + log_a_mean : float or ndarray + Characteristic log10(a / 1 AU) in dex (not ln, not AU). + Python float if ``mass`` is scalar, ndarray otherwise. + """ + # Calculate the characteristic log10(a / 1 AU) using a smooth broken power law + log_a_mean = _smooth_broken_loglog( + mass, self.a_mup, self.a_mp, + self.a_alphaL, self.a_alphaR, self.a_s, + a_min=self.a_min) + + return log_a_mean + + def a_mean(self, mass): + """ + Characteristic μ(a) in AU, ``10 ** log_a_mean(mass)``. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + a_mean : float or ndarray + Characteristic separation μ(a) in AU. Python float if + ``mass`` is scalar, ndarray otherwise. + """ + log_a = self.log_a_mean(mass) + + if np.isscalar(log_a): + a_mean = 10.0 ** log_a + return a_mean + + a_mean = 10.0 ** np.asarray(log_a, dtype=float) + + return a_mean + + def sigma_log_a(self, mass): + """ + σ(log10 a) from a 2-parameter logistic in log-mass. + + Floors/ceilings pinned at 0.7 / 1.5; clipped to ≥ 0.1:: + + σ(M) = 0.7 + 0.8 / (1 + (M / 0.354)**(-6.05)) + + i.e. (A, B, M0, k) = (0.7, 1.5, 0.354, 6.05). + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + + Returns + ------- + sigma_log_a : float or ndarray + Standard deviation of log10(a / 1 AU), in dex. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + + # Calculate the standard deviation of log10(a / 1 AU) using a 2-parameter logistic + sigma_log_a = _logistic_in_logm( + mass, self.sig_A, self.sig_B, self.sig_M0, self.sig_k, + clip_min=0.1) + + return sigma_log_a + + def _q_values_for_primaries(self, prim_subset, n_comp, rng): + """ + Draw mass-dependent q for every primary (BD and stellar). + + Parameters + ---------- + prim_subset : array_like + Primary masses for this companion-count group. Must be + positive, in solar masses (Msun). + n_comp : int + Number of companions per primary, integer count. + rng : numpy.random.Generator + Random generator used for inverse-CDF q draws. + + Returns + ------- + q_values : ndarray + Companion mass ratios m_comp/m_prim, dimensionless, in + [q_min, 1]. Shape (len(prim_subset), n_comp). + """ + # Draw the mass-dependent mass ratios using a logistic in log-mass + q_values = self.random_q(rng.random((len(prim_subset), n_comp)), mass=prim_subset) + + return q_values + + +class MultiplicityResolvedOffner2023(MultiplicityUnresolvedOffner2023, + _ResolvedOrbitalMixin): + """ + Opt-in resolved Offner et al. 2023 multiplicity. + + Same MF/CSF/q as :class:`MultiplicityUnresolvedOffner2023` + (Table 1 logistic MF/CSF, error-weighted γ logistic, Table 1 + companion-cut caveat). Adds mass-dependent separations. + Higher-order multiples are allowed at all masses. + Scientifically preferred over :class:`MultiplicityResolvedDK` + in the brown-dwarf regime, but **not** the default. + + Notes + ----- + :meth:`log_semimajoraxis` draws log10(a/AU) from a truncated + lognormal with ``loc = log_a_mean(mass)`` (smooth broken power + law, s = 0.1 dex, FGK-pulled) and + ``scale = sigma_log_a(mass)`` (2-parameter logistic pinned at + 0.7 / 1.5). Brown-dwarf binaries peak near a few AU + (μ(0.033) ≈ 2 AU). Truncation is 0.01–2000 AU, same limits as + :class:`MultiplicityResolvedDK`. + + Eccentricity and Keplerian angles still follow Duchêne & Kraus + (2013): f(e) = 2e, random inclination and angles. Same mixin as + :class:`MultiplicityResolvedDK`. + + Parameters + ---------- + MF_A, MF_B : float, optional + Low- and high-mass MF logistic asymptotes, dimensionless. + Defaults 0.14, 0.99 (Offner et al. 2023 Table 1, + equal-weight geom-mean MF fit). + MF_M0 : float, optional + Characteristic mass of the MF logistic, in solar masses + (Msun). Default 1.41. + MF_k : float, optional + MF logistic slope, dimensionless. Default 1.25. + CSF_A, CSF_B : float, optional + Low- and high-mass CSF logistic asymptotes, dimensionless + mean companion count. Defaults 0.12, 2.35 (Table 1 + equal-weight geom-mean CF fit). + CSF_M0 : float, optional + Characteristic mass of the CSF logistic, in solar masses + (Msun). Default 3.57. + CSF_k : float, optional + CSF logistic slope, dimensionless. Default 0.96. + q_A, q_B : float, optional + Low- and high-mass γ logistic asymptotes, dimensionless. + Defaults 6.6, −1.77 (Table 1 γ_trunc, error-weighted). + q_M0 : float, optional + Characteristic mass of the γ logistic, in solar masses + (Msun). Default 0.0651. + q_k : float, optional + γ logistic slope, dimensionless. Default 0.629. + a_mup : float, optional + Smooth-broken-PL pivot μ(a), in AU. Default 44.46. + a_mp : float, optional + Smooth-broken-PL pivot mass, in solar masses (Msun). + Default 0.819. + a_alphaL, a_alphaR : float, optional + Smooth-broken-PL slopes in log10 a vs log10 M, + dimensionless. Defaults 1.005, −0.308. + a_s : float, optional + Smooth-broken-PL smoothing scale, in dex. Default 0.10. + a_min : float, optional + Minimum characteristic a, in AU. Default 0.1. + sig_A, sig_B : float, optional + Low- and high-mass σ(log10 a) logistic values, in dex. + Defaults 0.7, 1.5 (Table 2 pins). + sig_M0 : float, optional + Characteristic mass of the σ logistic, in solar masses + (Msun). Default 0.354. + sig_k : float, optional + σ logistic slope, dimensionless. Default 6.05. + CSF_max : float, optional + Maximum companion star fraction, dimensionless mean companion + count (not bounded by 1). Default 3. + q_power : float, optional + Fallback mass-ratio power-law index, dimensionless. Ignored for + draws when primary mass is provided (the γ logistic is used). + Default 0.2. + q_min : float, optional + Minimum mass ratio m_comp/m_prim, dimensionless. Default 0.01. + companion_max : bool, optional + If True, cap companion counts at CSF_max at all masses. + Default False. + sep_sig : array_like, optional + Table 2 σ(log10 a) knots, in dex. Default (0.7, 1.3, 1.5). + sep_sig_mass : array_like, optional + Table 2 primary-mass knots for the σ comparison plot, in + solar masses (Msun). Defaults are geom-mean M1 of the + late-M, early-M, and FGK bins. + """ + def __init__(self, MF_A=0.14, MF_B=0.99, MF_M0=1.41, MF_k=1.25, + CSF_A=0.12, CSF_B=2.35, CSF_M0=3.57, CSF_k=0.96, + q_A=6.6, q_B=-1.77, q_M0=0.0651, q_k=0.629, + a_mup=44.46, a_mp=0.819, a_alphaL=1.005, + a_alphaR=-0.308, a_s=0.10, a_min=0.1, + sig_A=0.7, sig_B=1.5, sig_M0=0.354, sig_k=6.05, + CSF_max=3, q_power=0.2, q_min=0.01, + companion_max=False, + sep_sig=(0.7, 1.3, 1.5), + sep_sig_mass=(np.sqrt(0.075 * 0.15), + np.sqrt(0.3 * 0.6), + np.sqrt(0.75 * 1.25))): + super(MultiplicityResolvedOffner2023, self).__init__( + MF_A=MF_A, MF_B=MF_B, MF_M0=MF_M0, MF_k=MF_k, + CSF_A=CSF_A, CSF_B=CSF_B, CSF_M0=CSF_M0, CSF_k=CSF_k, + q_A=q_A, q_B=q_B, q_M0=q_M0, q_k=q_k, + a_mup=a_mup, a_mp=a_mp, a_alphaL=a_alphaL, + a_alphaR=a_alphaR, a_s=a_s, a_min=a_min, + sig_A=sig_A, sig_B=sig_B, sig_M0=sig_M0, sig_k=sig_k, + CSF_max=CSF_max, q_power=q_power, q_min=q_min, + companion_max=companion_max) + # Table 2 σ knots for the comparison plot; draws use sigma_log_a. + self.sep_sig_mass = np.array(sep_sig_mass, dtype=float) + self.sep_sig = np.array(sep_sig, dtype=float) + + return + + def log_semimajoraxis(self, mass, rng=None): + """ + Draw log10(a/AU) from a mass-dependent truncated lognormal. + + Uses ``loc = log_a_mean(mass)`` and + ``scale = sigma_log_a(mass)``. Truncated so a is between + 0.01 AU and 2000 AU (same limits as + :class:`MultiplicityResolvedDK`). + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + rng : numpy.random.Generator, optional + Random number generator passed to ``truncnorm.rvs``. + Default is a new ``numpy.random.default_rng()``. + + Returns + ------- + log_semimajoraxis : ndarray + Drawn log10(a / 1 AU) in dex (not ln, not AU), truncated + so a is between 0.01 AU and 2000 AU. + """ + if rng is None: + rng = np.random.default_rng() + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + log_a_mean = np.atleast_1d(np.asarray(self.log_a_mean(mass), + dtype=float)) + log_a_std = np.atleast_1d(np.asarray(self.sigma_log_a(mass), + dtype=float)) + + log_a_lower = np.log10(0.01) + log_a_upper = np.log10(2000) + a_lower_std = (log_a_lower - log_a_mean) / log_a_std + a_upper_std = (log_a_upper - log_a_mean) / log_a_std + + # Draw the log10(a / 1 AU) from a truncated normal + log_a = truncnorm.rvs(a_lower_std, a_upper_std, + loc=log_a_mean, scale=log_a_std, + random_state=rng) + + return log_a + + +def _piecewise_powerlaw(mass, mass_limits, amps, powers, clip_min=None, + clip_max=None): + """ + Evaluate y = A_i * M**alpha_i on mass segments. + + Segment i applies for mass_limits[i] <= M < mass_limits[i+1]. + The first segment also covers M below the lowest limit; the last + segment is closed on the right. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + mass_limits : array_like + Segment edges in solar masses (Msun), length N+1, strictly + increasing. + amps : array_like + Length-N amplitudes A_i, in units of y / Msun**alpha_i. + powers : array_like + Length-N power-law indices alpha_i, dimensionless. + clip_min, clip_max : float or None, optional + Optional lower/upper clips on y, in the same units as y. + ``None`` means no clip on that side. + + Returns + ------- + y : float or ndarray + Piecewise power-law value, in the same units as + ``amps * mass**powers``. Python float if ``mass`` is scalar, + ndarray otherwise. + """ + return_scalar = np.isscalar(mass) + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + out = np.empty(mass_arr.shape, dtype=float) + nseg = len(amps) + + # Evaluate the piecewise power-law for each segment + for i in range(nseg): + lo = mass_limits[i] + hi = mass_limits[i + 1] + + # Determine the mask for the current segment + if i == 0: + mask = mass_arr < hi + elif i == nseg - 1: + mask = mass_arr >= lo + else: + mask = (mass_arr >= lo) & (mass_arr < hi) + out[mask] = amps[i] * np.power(mass_arr[mask], powers[i]) + + # Apply the clips + if clip_min is not None: + out = np.maximum(out, clip_min) + + if clip_max is not None: + out = np.minimum(out, clip_max) + + # Return the result + if return_scalar: + return float(out[0]) + + return out + + +def _logistic_in_logm(mass, A, B, M0, k, clip_min=None, clip_max=None): + """ + Evaluate y = A + (B - A) / (1 + (M / M0)**(-k)). + + This is a logistic in log-mass: as M -> 0+, y -> A; as M -> inf, + y -> B. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + A, B : float + Low-mass and high-mass asymptotes, in the same units as y. + Dimensionless for MF and γ; mean companion count for CSF; + dex for σ(log10 a). + M0 : float + Characteristic mass in solar masses (Msun). + k : float + Logistic slope, dimensionless. + clip_min, clip_max : float or None, optional + Optional lower/upper clips on y, in the same units as y. + ``None`` means no clip on that side. + + Returns + ------- + y : float or ndarray + Logistic value in the same units as ``A`` and ``B``. + Python float if ``mass`` is scalar, ndarray otherwise. + """ + return_scalar = np.isscalar(mass) + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + + # Evaluate the logistic in log-mass + out = A + (B - A) / (1.0 + np.power(mass_arr / M0, -k)) + + # Apply the clips + if clip_min is not None: + out = np.maximum(out, clip_min) + if clip_max is not None: + out = np.minimum(out, clip_max) + + # Return the result + if return_scalar: + return float(out[0]) + + return out + + +def _logcosh(x): + """ + Numerically stable log(cosh(x)). + + Uses |x| + log1p(exp(-2|x|)) - log(2) rather than np.log(np.cosh(x)), + which overflows for |x| ≳ 700. + + Parameters + ---------- + x : float or array_like + Argument of cosh, dimensionless (for the smooth broken power + law this is v/s, with v and s in dex). + + Returns + ------- + logcosh_x : float or ndarray + log(cosh(x)), dimensionless. Same shape as ``x``. + """ + ax = np.abs(np.asarray(x, dtype=float)) + + # Calculate the log(cosh(x)) + logcosh_x = ax + np.log1p(np.exp(-2.0 * ax)) - np.log(2.0) + + # Return the result + return logcosh_x + + +def _smooth_broken_loglog(mass, mup, Mp, alpha_L, alpha_R, s, a_min=0.1): + """ + Smooth broken power law in log10(a) vs log10(M). + + v = log10(M / Mp) + yp = log10(mup) + log10(a) = yp + 0.5*(αL+αR)*v + 0.5*(αR-αL)*s * logcosh(v/s) + + ``s`` is the smoothing scale in dex (C-infinity; logcosh). The + linear-space value is clipped to ``a_min``. + + Parameters + ---------- + mass : float or array_like + Primary mass in solar masses (Msun). + mup : float + Characteristic separation at the break mass, in AU. + Mp : float + Break mass in solar masses (Msun). + alpha_L, alpha_R : float + Power-law indices below and above ``Mp``, dimensionless. + s : float + Smoothing scale in dex of log10(M / 1 Msun). + a_min : float, optional + Minimum linear-space separation in AU. Default 0.1 AU. + + Returns + ------- + log_a : float or ndarray + log10(a / 1 AU) in dex (not ln, not AU). Python float if + ``mass`` is scalar, ndarray otherwise. + """ + return_scalar = np.isscalar(mass) + mass_arr = np.atleast_1d(np.asarray(mass, dtype=float)) + v = np.log10(mass_arr / float(Mp)) + yp = np.log10(float(mup)) + + # Calculate the log10(a / 1 AU) using the smooth broken power law + log_a = (yp + + 0.5 * (alpha_L + alpha_R) * v + + 0.5 * (alpha_R - alpha_L) * s * _logcosh(v / s)) + a = np.maximum(10.0 ** log_a, float(a_min)) + log_a = np.log10(a) + + # Return the result + if return_scalar: + return float(log_a[0]) + + return log_a + + +def _q_from_powerlaw(x, q_pow, q_min): + """ + Inverse CDF of P(q) ∝ q**q_pow for q_min <= q <= 1. + + ``q_pow`` may be a scalar or an array broadcastable to ``x``. + The q_pow = -1 (b = 0) limit is q = q_min**(1 - x). + + Parameters + ---------- + x : float or array_like + Uniform random draw, dimensionless, in [0, 1]. + q_pow : float or array_like + Mass-ratio power-law index γ, dimensionless. Broadcastable + to ``x``. + q_min : float + Minimum mass ratio m_comp/m_prim, dimensionless, in (0, 1]. + + Returns + ------- + q : ndarray + Companion mass ratio m_comp/m_prim, dimensionless, in + [q_min, 1]. Same shape as the broadcast of ``x`` and ``q_pow``. + """ + x = np.asarray(x, dtype=float) + q_pow = np.asarray(q_pow, dtype=float) + + # Broadcast the arrays if necessary + if x.ndim > q_pow.ndim: + q_pow = q_pow.reshape(q_pow.shape + (1,) * (x.ndim - q_pow.ndim)) + + b = 1.0 + q_pow + b, x = np.broadcast_arrays(b, x) + + # Create an empty array to store the result + out = np.empty(x.shape, dtype=float) + + # Determine the mask for values near zero + near_zero = np.abs(b) < 1e-12 + + # Determine the mask for values far from zero + ok = ~near_zero + + # Calculate the mass ratio for values near zero + if np.any(near_zero): + out[near_zero] = q_min ** (1.0 - x[near_zero]) + if np.any(ok): + out[ok] = (x[ok] * (1.0 - q_min ** b[ok]) + q_min ** b[ok]) ** (1.0 / b[ok]) + + # Return the result + return out + diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 0fe624a8..092de2df 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -532,14 +532,24 @@ def _make_companions_table_initial(self, star_systems, compMass): companions = Table([system_index], names=['system_idx']) companions.add_column(np.zeros(N_comp_tot, dtype=float), name='mass') - if isinstance(self.imf._multi_props, multiplicity.MultiplicityResolvedDK): - companions.add_column(Column(self.imf._multi_props.log_semimajoraxis(star_systems['mass'][companions['system_idx']]), name='log_a')) - companions.add_column(Column(self.imf._multi_props.random_e(self.rng.random(N_comp_tot)), name='e')) - companions['i'], companions['Omega'], companions['omega'] = self.imf._multi_props.random_keplarian_parameters( - self.rng.random(N_comp_tot), - self.rng.random(N_comp_tot), - self.rng.random(N_comp_tot) - ) + # Duck-type resolved multiplicity: any object with orbital methods + # gets log_a / e / angles, not only MultiplicityResolvedDK. + multi_props = self.imf._multi_props + if (hasattr(multi_props, 'log_semimajoraxis') and + hasattr(multi_props, 'random_e') and + hasattr(multi_props, 'random_keplarian_parameters')): + prim_mass = star_systems['mass'][companions['system_idx']] + companions.add_column(Column( + multi_props.log_semimajoraxis(prim_mass, rng=self.rng), + name='log_a')) + companions.add_column(Column( + multi_props.random_e(self.rng.random(N_comp_tot)), name='e')) + companions['i'], companions['Omega'], companions['omega'] = ( + multi_props.random_keplarian_parameters( + self.rng.random(N_comp_tot), + self.rng.random(N_comp_tot), + self.rng.random(N_comp_tot), + rng=self.rng)) companions['mass'] = compMass.compressed() for key in ['Teff', 'L', 'logg', 'mass_current', 'phase']: diff --git a/spisea/tests/test_imf.py b/spisea/tests/test_imf.py index 943dc0a1..41623a25 100755 --- a/spisea/tests/test_imf.py +++ b/spisea/tests/test_imf.py @@ -30,6 +30,40 @@ def test_generate_cluster(): return + +def test_generate_cluster_offner2023(): + """ + generate_cluster with Offner et al. 2023 multiplicity: vectorized MF, + and Offner q (not Fontanive gamma=6.1) is used for BD companion + masses. + """ + imf_multi = multiplicity.MultiplicityUnresolvedOffner2023() + massLimits = np.array([0.01, 0.05, 0.22, 0.55, 8, 120]) + powers = np.array([-0.6, -0.25, -1.3, -2.3, -2.35]) + my_imf = imf.IMF_broken_powerlaw(massLimits, powers, imf_multi) + my_imf.rng = np.random.default_rng(7) + + M_cl = 2e3 + mass, isMulti, compMass, sysMass = my_imf.generate_cluster(M_cl) + + assert np.abs(M_cl - sysMass.sum()) < M_cl * 0.05 + bd = mass <= 0.08 + assert np.any(isMulti) + + # BD companions should be more equal-mass than SPISEA v2.5 stellar q_power=-0.4 + bd_mult = bd & isMulti + if np.any(bd_mult): + q_bd = [] + for i in np.where(bd_mult)[0]: + comps = compMass[i].compressed() + if len(comps): + q_bd.extend(list(comps / mass[i])) + if len(q_bd) >= 5: + assert np.mean(q_bd) > 0.5 + + return + + def test_prim_power(): #mass_limits = np.array([0.1, 1.0, 100.0]) #powers = np.array([-2.0, -1.8]) diff --git a/spisea/tests/test_multiplicity.py b/spisea/tests/test_multiplicity.py index 53226651..a7ef111f 100755 --- a/spisea/tests/test_multiplicity.py +++ b/spisea/tests/test_multiplicity.py @@ -1,5 +1,6 @@ import numpy as np import time +import os import spisea from spisea.imf import imf, multiplicity @@ -131,13 +132,13 @@ def test_companion_star_fraction(): # csf2_3 = mu1.companion_star_fraction(0.1) # np.testing.assert_almost_equal(csf2_3, 0.159, decimal=2) - # Test brown dwarf csf + # BD CSF follows the stellar power law (no CSF=MF mass cut). csf_bd1 = mu1.companion_star_fraction(0.07) csf_bd2 = mu1.companion_star_fraction(0.04) csf_bd3 = mu1.companion_star_fraction(0.01) - assert np.isclose(csf_bd1, 0.16, atol=0.01) - assert np.isclose(csf_bd2, 0.08, atol=0.01) - assert np.isclose(csf_bd3, 0.0, atol=1e-6) + np.testing.assert_almost_equal(csf_bd1, 0.50 * 0.07 ** 0.45, decimal=4) + np.testing.assert_almost_equal(csf_bd2, 0.50 * 0.04 ** 0.45, decimal=4) + np.testing.assert_almost_equal(csf_bd3, 0.50 * 0.01 ** 0.45, decimal=4) def test_resolvedmult(): @@ -243,3 +244,526 @@ def test_resolvedmult(): f"BD sigma log(a) off: {std_log_a:.2f}" return + + +# --------------------------------------------------------------------------- +# Offner et al. 2023 (Table 1) multiplicity +# --------------------------------------------------------------------------- + +# Published Table 1 MF (%) converted to fraction, CF, and 1-sigma MF error. +# Masses are geometric means of the tabulated M1 intervals. +_OFFNER_TABLE1 = [ + # name, M_lo, M_hi, MF, MF_err, CF + ('Fontanive+2018', 0.019, 0.058, 0.08, 0.06, 0.08), + ('Burgasser 2007', 0.05, 0.08, 0.15, 0.04, 0.16), + ('Close+2003', 0.080, 0.095, 0.19, 0.07, 0.19), + ('Allen+2007', 0.06, 0.15, 0.20, 0.04, 0.20), + ('Winters+2019 late-M', 0.075, 0.15, 0.19, 0.03, 0.21), + ('Winters+2019 mid-M', 0.15, 0.30, 0.23, 0.02, 0.27), + ('Winters+2019 early-M', 0.3, 0.6, 0.30, 0.02, 0.38), + ('Raghavan+2010', 0.75, 1.25, 0.46, 0.03, 0.60), + ('Tokovinin 2014b', 0.85, 1.5, 0.47, 0.03, 0.62), + ('Moe & Kratter 2021', 1.6, 2.4, 0.68, 0.07, 0.99), + ('Moe & Di Stefano 2017 3-5', 3.0, 5.0, 0.81, 0.06, 1.28), + ('Moe & Di Stefano 2017 5-8', 5.0, 8.0, 0.89, 0.05, 1.55), + ('Moe & Di Stefano 2017 8-17', 8.0, 17.0, 0.93, 0.04, 1.80), + ('Sana et al. 17-50', 17.0, 50.0, 0.96, 0.04, 2.10), +] + + +def _table1_mgeom(row): + return np.sqrt(row[1] * row[2]) + + +def test_piecewise_powerlaw_api(): + """Custom piecewise MF/CSF is vectorized and clips MF to [0, 1].""" + mass_limits = np.array([0.1, 1.0, 10.0]) + # First segment: MF = 0.4 * M^0 → 0.4; second: 0.4 * M^1 so MF(10)=4 → clip 1 + mp = multiplicity.MultiplicityPiecewisePowerLaw( + mass_limits, + MF_amps=[0.4, 0.4], MF_powers=[0.0, 1.0], + CSF_amps=[0.4, 0.5], CSF_powers=[0.0, 0.5]) + assert mp.multiplicity_fraction(0.2) == 0.4 + assert mp.multiplicity_fraction(1.0) == 0.4 + np.testing.assert_almost_equal(mp.multiplicity_fraction(10.0), 1.0) + masses = np.array([0.2, 1.0, 10.0]) + mf = mp.multiplicity_fraction(masses) + np.testing.assert_allclose(mf, [mp.multiplicity_fraction(m) for m in masses]) + + +def test_logistic_api(): + """Custom logistic MF/CSF clips, vectorizes, and keeps CSF >= MF.""" + ml = multiplicity.MultiplicityLogistic( + MF_A=0.1, MF_B=1.5, MF_M0=1.0, MF_k=2.0, + CSF_A=0.2, CSF_B=4.0, CSF_M0=2.0, CSF_k=1.0, + CSF_max=2.0) + # Low-mass asymptote A for a very low-mass primary (not a missing mass) + np.testing.assert_almost_equal(ml.multiplicity_fraction(1e-8), 0.1, decimal=4) + # High-mass MF saturates at B then clips to 1 + np.testing.assert_almost_equal(ml.multiplicity_fraction(1e6), 1.0) + # High-mass CSF clips to CSF_max + np.testing.assert_almost_equal(ml.companion_star_fraction(1e6), 2.0) + # No CSF=MF mass cut: low-mass CSF can exceed MF. + assert ml.companion_star_fraction(0.05) > ml.multiplicity_fraction(0.05) + masses = np.array([0.05, 1.0, 100.0]) + mf = ml.multiplicity_fraction(masses) + np.testing.assert_allclose(mf, [ml.multiplicity_fraction(m) for m in masses]) + csf = ml.companion_star_fraction(masses) + np.testing.assert_allclose( + csf, [ml.companion_star_fraction(m) for m in masses]) + assert np.all(csf >= mf - 1e-12) + + +def test_offner2023_logistic_coefficients(): + """Offner stores the equal-weight logistic-in-log-mass coefficients.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + assert isinstance(multi, multiplicity.MultiplicityLogistic) + assert not isinstance(multi, multiplicity.MultiplicityPiecewisePowerLaw) + np.testing.assert_allclose(multi.MF_A, 0.14) + np.testing.assert_allclose(multi.MF_B, 0.99) + np.testing.assert_allclose(multi.MF_M0, 1.41) + np.testing.assert_allclose(multi.MF_k, 1.25) + np.testing.assert_allclose(multi.CSF_A, 0.12) + np.testing.assert_allclose(multi.CSF_B, 2.35) + np.testing.assert_allclose(multi.CSF_M0, 3.57) + np.testing.assert_allclose(multi.CSF_k, 0.96) + np.testing.assert_allclose(multi.q_A, 6.6) + np.testing.assert_allclose(multi.q_B, -1.77) + np.testing.assert_allclose(multi.q_M0, 0.0651) + np.testing.assert_allclose(multi.q_k, 0.629) + np.testing.assert_allclose(multi.a_mup, 44.46) + np.testing.assert_allclose(multi.a_mp, 0.819) + np.testing.assert_allclose(multi.a_alphaL, 1.005) + np.testing.assert_allclose(multi.a_alphaR, -0.308) + np.testing.assert_allclose(multi.a_s, 0.10) + np.testing.assert_allclose(multi.a_min, 0.1) + np.testing.assert_allclose(multi.sig_A, 0.7) + np.testing.assert_allclose(multi.sig_B, 1.5) + np.testing.assert_allclose(multi.sig_M0, 0.354) + np.testing.assert_allclose(multi.sig_k, 6.05) + assert not hasattr(multiplicity, 'FONTANIVE2018_BD_Q_POWER') + assert not any(n.startswith('OFFNER2023_') for n in dir(multiplicity)) + resolved = multiplicity.MultiplicityResolvedOffner2023() + np.testing.assert_allclose(resolved.MF_A, 0.14) + np.testing.assert_allclose(resolved.sig_k, 6.05) + import inspect + res_sig = inspect.signature( + multiplicity.MultiplicityResolvedOffner2023.__init__) + assert res_sig.parameters['MF_A'].default == 0.14 + assert res_sig.parameters['sig_k'].default == 6.05 + + +def test_offner2023_mf_smooth(): + """MF is continuous and nearly C1 around 0.08 and 1.5 Msun.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + eps = 1e-8 + for m in (0.08, 1.5): + mf_left = multi.multiplicity_fraction(m - eps) + mf_right = multi.multiplicity_fraction(m + eps) + mf_at = multi.multiplicity_fraction(m) + np.testing.assert_allclose(mf_left, mf_right, atol=1e-6, rtol=0) + np.testing.assert_allclose(mf_at, mf_right, atol=1e-6, rtol=0) + d_left = (mf_at - mf_left) / eps + d_right = (mf_right - mf_at) / eps + np.testing.assert_allclose(d_left, d_right, atol=1e-3, rtol=0) + for m in (0.04, 0.3, 1.0, 10.0): + expected = multiplicity._logistic_in_logm( + m, multi.MF_A, multi.MF_B, multi.MF_M0, multi.MF_k, + clip_min=0.0, clip_max=1.0) + np.testing.assert_allclose(multi.multiplicity_fraction(m), expected) + + +def test_offner2023_table1_mf(): + """ + Logistic MF matches Offner et al. 2023 Table 1 at geom-mean M1. + + Fontanive (8±6%) sits ~0.07 below the curve (~15%); other rows, + including A/B stars, stay close. + """ + multi = multiplicity.MultiplicityUnresolvedOffner2023() + for row in _OFFNER_TABLE1: + name, mlo, mhi, mf_tab, mf_err, cf_tab = row + m = _table1_mgeom(row) + mf = multi.multiplicity_fraction(m) + tol = max(0.08, 2.0 * mf_err) + assert abs(mf - mf_tab) <= tol, \ + '{0}: MF({1:.3f})={2:.3f} vs Table 1 {3:.2f} ± {4:.2f}'.format( + name, m, mf, mf_tab, mf_err) + assert 0.0 <= mf <= 1.0 + + +def test_offner2023_table1_csf(): + """CSF tracks Table 1 CF; CSF >= MF at all masses (no BD CSF=MF cut).""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + for row in _OFFNER_TABLE1: + name, mlo, mhi, mf_tab, mf_err, cf_tab = row + m = _table1_mgeom(row) + csf = multi.companion_star_fraction(m) + mf = multi.multiplicity_fraction(m) + if mlo >= 1.6: + # Logistic CF tracks A/B; Moe & Kratter residual ~0.1 is ok + tol = max(0.12, 0.12 * cf_tab) + else: + tol = max(0.08, 0.25 * cf_tab) + assert abs(csf - cf_tab) <= tol, \ + '{0}: CSF({1:.3f})={2:.3f} vs Table 1 CF {3:.2f}'.format( + name, m, csf, cf_tab) + assert csf >= mf - 1e-12 + assert csf <= multi.CSF_max + 1e-12 + + +def test_offner2023_array_vs_scalar(): + """Array and scalar MF/CSF evaluations agree.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + masses = np.array([_table1_mgeom(row) for row in _OFFNER_TABLE1]) + mf_arr = multi.multiplicity_fraction(masses) + csf_arr = multi.companion_star_fraction(masses) + for i, m in enumerate(masses): + np.testing.assert_allclose(mf_arr[i], multi.multiplicity_fraction(float(m))) + np.testing.assert_allclose(csf_arr[i], multi.companion_star_fraction(float(m))) + + +def test_higher_order_multiples_at_all_masses(): + """BD and stellar primaries can have n_comp > 1 when CSF/MF allows.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + rng = np.random.default_rng(123) + for m in (0.04, 1.0): + masses = np.full(4000, m) + mf = multi.multiplicity_fraction(masses) + csf = np.full(len(masses), 2.5) + n_comp = multi.draw_n_companions(masses, csf, mf, rng) + assert np.any(n_comp > 1) + assert np.all(n_comp >= 1) + + +def test_companion_max_caps_all_masses(): + """companion_max=True clips n_comp at CSF_max for BD and stellar.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023( + companion_max=True, CSF_max=1) + rng = np.random.default_rng(0) + for m in (0.04, 1.0): + masses = np.full(2000, m) + mf = multi.multiplicity_fraction(masses) + csf = np.full(len(masses), 3.0) + n_comp = multi.draw_n_companions(masses, csf, mf, rng) + assert np.all(n_comp <= 1) + assert np.all(n_comp >= 1) + + +def test_offner2023_q_more_equal_mass_for_bds(): + """BD mass ratios are more equal-mass (higher mean q) than solar-type.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + rng = np.random.default_rng(7) + n = 20000 + q_bd = multi.random_q(rng.random(n), mass=0.04) + q_sun = multi.random_q(rng.random(n), mass=1.0) + assert np.mean(q_bd) > np.mean(q_sun) + 0.1 + # Err-wt logistic undershoots Fontanive 4.8 (~3.3 at 0.033 Msun) + assert multi.q_power_at_mass(0.033) > 2.5 + assert multi.q_power_at_mass(1.0) < 0.5 + + +def test_offner2023_q_sigma_a_closed_form(): + """γ, σ(log a), and log_a_mean match the smooth helpers; not interpolation.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + masses = np.array([0.033, 0.065, 0.3, 1.0, 10.0]) + for m in masses: + np.testing.assert_allclose( + multi.q_power_at_mass(m), + multiplicity._logistic_in_logm( + m, multi.q_A, multi.q_B, multi.q_M0, multi.q_k)) + np.testing.assert_allclose( + multi.sigma_log_a(m), + multiplicity._logistic_in_logm( + m, multi.sig_A, multi.sig_B, multi.sig_M0, multi.sig_k, + clip_min=0.1)) + np.testing.assert_allclose( + multi.log_a_mean(m), + multiplicity._smooth_broken_loglog( + m, multi.a_mup, multi.a_mp, multi.a_alphaL, multi.a_alphaR, + multi.a_s, a_min=multi.a_min)) + # Array vs scalar + g_arr = multi.q_power_at_mass(masses) + sig_arr = multi.sigma_log_a(masses) + loga_arr = multi.log_a_mean(masses) + for i, m in enumerate(masses): + np.testing.assert_allclose(g_arr[i], multi.q_power_at_mass(float(m))) + np.testing.assert_allclose(sig_arr[i], multi.sigma_log_a(float(m))) + np.testing.assert_allclose(loga_arr[i], multi.log_a_mean(float(m))) + # Old L/early-T interpolation knot was 2.5; logistic is not that. + g_knot = multi.q_power_at_mass(0.065) + np.testing.assert_allclose( + g_knot, multiplicity._logistic_in_logm( + 0.065, multi.q_A, multi.q_B, multi.q_M0, multi.q_k)) + assert abs(g_knot - 2.5) > 0.05 + + +def test_offner2023_bd_separations_peak_few_au(): + """BD lognormal separations peak at a few AU (μ(0.033)≈2.1 au).""" + multi = multiplicity.MultiplicityResolvedOffner2023() + rng = np.random.default_rng(0) + log_a = multi.log_semimajoraxis(np.full(5000, 0.04), rng=rng) + med_a = 10 ** np.median(log_a) + assert 1.5 < med_a < 8.0, 'BD median a = {0:.2f} AU'.format(med_a) + # Solar-type should be much wider (smooth-broken μ ~ 44 au) + log_a_s = multi.log_semimajoraxis(np.full(5000, 1.0), rng=rng) + med_a_s = 10 ** np.median(log_a_s) + assert med_a_s > 10.0 + assert med_a_s > med_a + + +def test_offner2023_resolved_methods(): + """Resolved orbital methods exist; no unresolved/resolved alias.""" + assert not hasattr(multiplicity, 'MultiplicityOffner2023') + resolved = multiplicity.MultiplicityResolvedOffner2023() + assert hasattr(resolved, 'log_semimajoraxis') + assert hasattr(resolved, 'log_a_mean') + assert hasattr(resolved, 'sigma_log_a') + np.testing.assert_allclose( + resolved.sep_sig_mass, + [np.sqrt(0.075 * 0.15), np.sqrt(0.3 * 0.6), np.sqrt(0.75 * 1.25)]) + np.testing.assert_allclose(resolved.sep_sig, [0.7, 1.3, 1.5]) + e = resolved.random_e(np.array([0.0, 0.25, 1.0])) + np.testing.assert_allclose(e, [0.0, 0.5, 1.0]) + + +def test_lu2013_defaults_unchanged(): + """SPISEA v2.5 MultiplicityUnresolved defaults and stellar MF unchanged.""" + mu = multiplicity.MultiplicityUnresolved() + assert mu.MF_amp == 0.44 + assert mu.MF_pow == 0.51 + assert mu.CSF_amp == 0.50 + assert mu.CSF_pow == 0.45 + np.testing.assert_almost_equal(mu.multiplicity_fraction(1.0), 0.44, decimal=2) + np.testing.assert_almost_equal(mu.multiplicity_fraction(10.0), 1.0, decimal=2) + np.testing.assert_almost_equal(mu.multiplicity_fraction(0.1), 0.136, decimal=2) + assert mu.bd_q_power == 6.1 + assert np.isclose(mu.q_power_at_mass(0.04), 6.1) + assert np.isclose(mu.q_power_at_mass(1.0), -0.4) + # Scalar BD overrides (SPISEA v2.5 / Fontanive path) + assert np.isclose(mu.multiplicity_fraction(0.07), 0.16, atol=0.01) + assert np.isclose(mu.multiplicity_fraction(0.04), 0.08, atol=0.01) + assert np.isclose(mu.multiplicity_fraction(0.01), 0.0, atol=1e-6) + + +def test_resolveddk_log_a_mean_and_sigma(): + """ + MultiplicityResolvedDK.a_mean / log_a_mean / sigma_log_a match the + Duchêne & Kraus + Fontanive BD interp + sigmoid that + log_semimajoraxis used to inline, and the draw uses those methods. + """ + import astropy.modeling + import inspect + + dk = multiplicity.MultiplicityResolvedDK() + masses = np.array([0.01, 0.04, 0.08, 0.3, 1.0, 5.0, 50.0]) + + def expected_log_a_mean(mass): + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + logm = np.log10(mass) + a_mean_func = astropy.modeling.powerlaws.BrokenPowerLaw1D( + amplitude=dk.a_amp, x_break=dk.a_break, + alpha_1=dk.a_slope1, alpha_2=dk.a_slope2) + log_a_mean_star = np.log10(a_mean_func(mass)) + log_a_mean_bd = np.interp( + logm, + [np.log10(0.01), np.log10(0.08)], + [np.log10(2.5), np.log10(8.0)]) + w = 1.0 / (1.0 + np.exp(-(logm - np.log10(0.08)) / 0.15)) + return (1 - w) * log_a_mean_bd + w * log_a_mean_star + + def expected_sigma(mass): + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + logm = np.log10(mass) + log_a_std_func = astropy.modeling.models.Linear1D( + slope=dk.a_std_slope, intercept=dk.a_std_intercept) + log_a_std_star = log_a_std_func(logm) + log_a_std_star[mass >= 2.9] = log_a_std_func(np.log10(2.9)) + log_a_std_star = np.clip(log_a_std_star, 0.1, None) + log_a_std_bd = np.interp( + logm, + [np.log10(0.01), np.log10(0.08)], + [0.25, 0.5]) + w = 1.0 / (1.0 + np.exp(-(logm - np.log10(0.08)) / 0.15)) + return (1 - w) * log_a_std_bd + w * log_a_std_star + + np.testing.assert_allclose(dk.log_a_mean(masses), expected_log_a_mean(masses)) + np.testing.assert_allclose(dk.sigma_log_a(masses), expected_sigma(masses)) + np.testing.assert_allclose(dk.a_mean(masses), 10.0 ** dk.log_a_mean(masses)) + + for m in masses: + log_a = dk.log_a_mean(float(m)) + sig = dk.sigma_log_a(float(m)) + a_mean = dk.a_mean(float(m)) + assert isinstance(log_a, float) + assert isinstance(sig, float) + assert isinstance(a_mean, float) + np.testing.assert_allclose(log_a, expected_log_a_mean(m)[0]) + np.testing.assert_allclose(sig, expected_sigma(m)[0]) + np.testing.assert_allclose(a_mean, 10.0 ** log_a) + + src = inspect.getsource(dk.log_semimajoraxis) + assert 'self.log_a_mean' in src + assert 'self.sigma_log_a' in src + assert 'BrokenPowerLaw1D' not in src + + log_a_draw = dk.log_semimajoraxis( + np.full(4000, 1.0), rng=np.random.default_rng(0)) + assert np.all(log_a_draw >= np.log10(0.01) - 1e-12) + assert np.all(log_a_draw <= np.log10(2000) + 1e-12) + np.testing.assert_allclose( + np.median(log_a_draw), dk.log_a_mean(1.0), atol=0.15) + + +def test_offner_generate_cluster_companions(): + """IMF cluster generation with Offner multiplicity produces companions.""" + imf_multi = multiplicity.MultiplicityUnresolvedOffner2023() + mass_limits = np.array([0.01, 0.08, 0.5, 120.0]) + powers = np.array([-0.3, -1.3, -2.3]) + my_imf = imf.IMF_broken_powerlaw(mass_limits, powers, imf_multi) + my_imf.rng = np.random.default_rng(42) + mass, is_multi, comp_mass, sys_mass = my_imf.generate_cluster(500.0) + assert np.any(is_multi) + assert np.abs(500.0 - sys_mass.sum()) < 500.0 * 0.05 + + +def test_calc_multi_uses_multiplicity_q_and_counts(): + """ + IMF.calc_multi must not hardcode Fontanive gamma=6.1; q policy + lives on the multiplicity object so Offner γ_trunc (~2–5 for + BDs) actually applies. + """ + import inspect + from spisea.imf import imf as imf_mod + calc_src = inspect.getsource(imf_mod.IMF.calc_multi) + assert '6.1' not in calc_src + assert 'draw_companion_masses' in calc_src + + syn_path = os.path.join(os.path.dirname(spisea.__file__), 'synthetic.py') + with open(syn_path, 'r') as fh: + syn_src = fh.read() + assert 'isinstance(self.imf._multi_props, multiplicity.MultiplicityResolvedDK)' not in syn_src + assert "hasattr(multi_props, 'log_semimajoraxis')" in syn_src + assert "hasattr(multi_props, 'random_e')" in syn_src + assert "hasattr(multi_props, 'random_keplarian_parameters')" in syn_src + assert 'log_semimajoraxis(prim_mass, rng=self.rng)' in syn_src + assert 'random_keplarian_parameters(' in syn_src + assert 'rng=self.rng)' in syn_src + + offner = multiplicity.MultiplicityUnresolvedOffner2023() + lu = multiplicity.MultiplicityUnresolved() + q_off = offner.q_power_at_mass(0.04) + q_lu = lu.q_power_at_mass(0.04) + assert 2.0 <= q_off <= 5.5 + assert np.isclose(q_lu, 6.1) + assert q_off != q_lu + + rng = np.random.default_rng(1) + masses = np.full(3000, 0.04) + is_mult = np.ones(len(masses), dtype=bool) + mf = offner.multiplicity_fraction(masses) + csf = offner.companion_star_fraction(masses) + comp, _, _ = offner.draw_companion_masses( + masses, is_mult, csf, mf, rng, mass_min=0.01) + q = comp.compressed() / 0.04 + q_lu_draw = lu.random_q(np.random.default_rng(1).random(len(q)), mass=0.04) + # Offner BD gamma is shallower than Fontanive 6.1, so mean q is lower. + assert np.mean(q) < np.mean(q_lu_draw) + + +def test_offner2023_mf_is_vectorized(): + """Offner MF is vectorized; SPISEA v2.5 scalar BD bins are not used here.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + masses = np.array([0.03, 0.06, 0.10, 1.0]) + mf = multi.multiplicity_fraction(masses) + for i, m in enumerate(masses): + np.testing.assert_allclose(mf[i], multi.multiplicity_fraction(float(m))) + # Array path is not the SPISEA v2.5 stellar power law 0.44 * M**0.51 + lu_pl = 0.44 * masses ** 0.51 + assert not np.allclose(mf, np.clip(lu_pl, 0, 1), atol=0.02) + + +def _same_seed_rngs(seed=11): + return np.random.default_rng(seed), np.random.default_rng(seed) + + +def test_random_companion_count_same_seed(): + """Same rng seed reproduces random_companion_count.""" + masses = np.full(40, 1.0) + for cls in (multiplicity.MultiplicityUnresolved, + multiplicity.MultiplicityUnresolvedOffner2023): + multi = cls() + mf = multi.multiplicity_fraction(masses) + csf = multi.companion_star_fraction(masses) + rng1, rng2 = _same_seed_rngs() + n1 = multi.random_companion_count( + None, csf, mf, mass=masses, rng=rng1) + n2 = multi.random_companion_count( + None, csf, mf, mass=masses, rng=rng2) + np.testing.assert_array_equal(n1, n2) + + +def test_draw_n_companions_same_seed(): + """Same rng seed reproduces draw_n_companions.""" + masses = np.array([0.04, 0.3, 1.0, 5.0] * 10) + for cls in (multiplicity.MultiplicityUnresolved, + multiplicity.MultiplicityUnresolvedOffner2023): + multi = cls() + mf = multi.multiplicity_fraction(masses) + csf = multi.companion_star_fraction(masses) + rng1, rng2 = _same_seed_rngs() + n1 = multi.draw_n_companions(masses, csf, mf, rng1) + n2 = multi.draw_n_companions(masses, csf, mf, rng2) + np.testing.assert_array_equal(n1, n2) + + +def test_assign_companions_same_seed(): + """Same rng seed reproduces draw_companion_masses assignment.""" + masses = np.array([0.04, 0.3, 1.0, 2.0, 5.0] * 8) + is_mult = np.ones(len(masses), dtype=bool) + for cls in (multiplicity.MultiplicityUnresolved, + multiplicity.MultiplicityUnresolvedOffner2023): + multi = cls() + mf = multi.multiplicity_fraction(masses) + csf = multi.companion_star_fraction(masses) + rng1, rng2 = _same_seed_rngs() + c1, s1, m1 = multi.draw_companion_masses( + masses, is_mult, csf, mf, rng1, mass_min=0.01) + c2, s2, m2 = multi.draw_companion_masses( + masses, is_mult, csf, mf, rng2, mass_min=0.01) + np.testing.assert_array_equal(np.ma.getdata(c1), np.ma.getdata(c2)) + np.testing.assert_array_equal(np.ma.getmaskarray(c1), + np.ma.getmaskarray(c2)) + np.testing.assert_array_equal(s1, s2) + np.testing.assert_array_equal(m1, m2) + + +def test_log_semimajoraxis_same_seed(): + """Same rng seed reproduces log_semimajoraxis.""" + masses = np.full(30, 1.0) + for cls in (multiplicity.MultiplicityResolvedDK, + multiplicity.MultiplicityResolvedOffner2023): + multi = cls() + rng1, rng2 = _same_seed_rngs() + a1 = multi.log_semimajoraxis(masses, rng=rng1) + a2 = multi.log_semimajoraxis(masses, rng=rng2) + np.testing.assert_array_equal(a1, a2) + + +def test_random_keplarian_parameters_same_seed(): + """Same rng seed reproduces random_keplarian_parameters.""" + n = 30 + for cls in (multiplicity.MultiplicityResolvedDK, + multiplicity.MultiplicityResolvedOffner2023): + multi = cls() + rng1, rng2 = _same_seed_rngs() + x1, y1, z1 = rng1.random(n), rng1.random(n), rng1.random(n) + i1, O1, o1 = multi.random_keplarian_parameters( + x1, y1, z1, rng=rng1) + x2, y2, z2 = rng2.random(n), rng2.random(n), rng2.random(n) + i2, O2, o2 = multi.random_keplarian_parameters( + x2, y2, z2, rng=rng2) + np.testing.assert_array_equal(i1, i2) + np.testing.assert_array_equal(O1, O2) + np.testing.assert_array_equal(o1, o2) +