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 0000000..b00e27a 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 0000000..20da517 Binary files /dev/null and b/docs/figures/mf_offner_vs_spisea2.5.png differ 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 0000000..be66ced --- /dev/null +++ b/docs/figures/plot_mf_offner_vs_spisea2.5.py @@ -0,0 +1,148 @@ +#!/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 array power law, v2.5 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 _spisea25_array_mf(mass): + """SPISEA v2.5 array path: stellar power law, no BD staircase.""" + mf = 0.44 * np.asarray(mass, dtype=float) ** 0.51 + return np.clip(mf, 0.0, 1.0) + + +def _spisea25_scalar_bd_staircase(mass): + """SPISEA v2.5 scalar-only BD bins in MultiplicityUnresolved.multiplicity_fraction.""" + mass = np.asarray(mass, dtype=float) + mf = _spisea25_array_mf(mass) + mf = np.where(mass < 0.02, 0.0, mf) + mf = np.where((mass > 0.02) & (mass <= 0.06), 0.08, mf) + mf = np.where((mass > 0.06) & (mass <= 0.08), 0.16, mf) + return mf + + +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() + + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + mf_off = offner.multiplicity_fraction(m_wide) + mf_lu = _spisea25_array_mf(m_wide) + + # Staircase sampled densely so the steps render as vertical jumps. + 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 = _spisea25_scalar_bd_staircase(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_sep_offner_vs_spisea2.5.py b/docs/figures/plot_q_sep_offner_vs_spisea2.5.py new file mode 100644 index 0000000..8a1dd62 --- /dev/null +++ b/docs/figures/plot_q_sep_offner_vs_spisea2.5.py @@ -0,0 +1,362 @@ +#!/usr/bin/env python +""" +Generate q / separation comparison figures vs SPISEA v2.5 defaults: + + docs/figures/q_offner_vs_spisea2.5.png + docs/figures/sep_offner_vs_spisea2.5.png + docs/figures/sig_loga_offner_vs_spisea2.5.png + docs/figures/meanq_offner_vs_spisea2.5.png + +Two-panel layout matching ``plot_mf_offner_vs_spisea2.5.py``: brown-dwarf +zoom and BD through O. Offner curves are evaluated from the +multiplicity objects (γ logistic, smooth-broken μ(a), σ logistic) +so they cannot drift from the code. + +Run from the repository root:: + + python docs/figures/plot_q_sep_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 + + +# 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), +] + +_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) + + +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 _spisea25_dk_mean_a_au(dk, mass): + """ + SPISEA v2.5 ``MultiplicityResolvedDK`` mean a in AU (Duchêne & Kraus + coefficients, plus the BD log-a interpolation and sigmoid blend in + ``log_semimajoraxis``). + """ + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + logm = np.log10(mass) + x = mass / dk.a_break + a_star = np.where( + mass < dk.a_break, + dk.a_amp * np.power(x, -dk.a_slope1), + dk.a_amp * np.power(x, -dk.a_slope2), + ) + a_star = np.maximum(a_star, 1e-30) + log_a_mean_star = np.log10(a_star) + 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.0 - w) * log_a_mean_bd + w * log_a_mean_star + return 10.0 ** log_a_mean + + +def _spisea25_dk_sig_loga(dk, mass): + """ + SPISEA v2.5 ``MultiplicityResolvedDK`` σ(log10 a) (Duchêne & Kraus + coefficients, plus the BD interpolation and sigmoid blend in + ``log_semimajoraxis``). + """ + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + logm = np.log10(mass) + sig_star = dk.a_std_slope * logm + dk.a_std_intercept + sig_star = np.array(sig_star, dtype=float, copy=True) + sig_star[mass >= 2.9] = dk.a_std_slope * np.log10(2.9) + dk.a_std_intercept + sig_star = np.clip(sig_star, 0.1, None) + sig_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.0 - w) * sig_bd + w * sig_star + + +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) + out = np.empty(np.shape(g), dtype=float) + 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 _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(__file__), filename) + fig.savefig(out, dpi=_DPI, bbox_inches='tight', facecolor='white') + plt.close(fig) + print('Wrote', out) + + +def _spisea25_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 plot_gamma(offner, lu): + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + g_off = offner.q_power_at_mass(m_wide) + m_step = _spisea25_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') + + +def plot_sep(resolved, dk): + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + a_off = resolved.a_mean(m_wide) + a_lu = _spisea25_dk_mean_a_au(dk, 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') + + +def plot_sig_loga(resolved, dk): + m_wide = np.logspace(np.log10(0.012), np.log10(40.0), 800) + sig_off = resolved.sigma_log_a(m_wide) + sig_lu = _spisea25_dk_sig_loga(dk, 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') + + +def plot_meanq(offner, lu): + 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 = _spisea25_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') + + +def main(): + offner = multiplicity.MultiplicityUnresolvedOffner2023() + lu = multiplicity.MultiplicityUnresolved() + resolved = multiplicity.MultiplicityResolvedOffner2023() + dk = multiplicity.MultiplicityResolvedDK() + plot_gamma(offner, lu) + plot_sep(resolved, dk) + plot_sig_loga(resolved, dk) + plot_meanq(offner, lu) + + +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 0000000..059743b 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 0000000..e34f2ba 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 0000000..792d243 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 ed60177..91c2946 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 9a0f740..fdfe005 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 5417738..2ba3df9 100644 --- a/docs/multiplicity.rst +++ b/docs/multiplicity.rst @@ -7,13 +7,26 @@ 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 keeps the + historical stellar-only power law. +* ``random_companion_count(x, CSF, MF, mass=None, rng=None)`` — + companion-count policy, including the binaries-only BD cap when + ``mass`` is given. +* attributes ``companion_max``, ``CSF_max``, ``q_min`` The user can choose either an unresolved or a resolved multiplicity object. If a resolved @@ -27,8 +40,8 @@ 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 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. @@ -40,8 +53,231 @@ Unresolved Multiplicity Classes :members: companion_star_fraction, multiplicity_fraction, random_q +.. 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 + + +Offner et al. 2023 multiplicity +------------------------------------------ +The recommended multiplicity class to use is that derived from +data summarized in Offner et al. (2023) (`arXiv:2203.10066 +`_; ADS +`2023ASPC..534..275O +`_). Table 1 +data: Zenodo `10.5281/zenodo.6628915 +`_. +This class is not the default (for backwards compatability) but +is strongly preferred. + + +The SPISEA v2.5 :class:`~imf.multiplicity.MultiplicityUnresolved` / +:class:`~imf.multiplicity.MultiplicityResolvedDK` objects remain the +default; but has known limitations in the brown dwarf regime. + +Unresolved (companions, no orbits):: + + from spisea.imf import imf, multiplicity + from spisea import synthetic + + multi = multiplicity.MultiplicityUnresolvedOffner2023() + # alias: + # multi = multiplicity.MultiplicityOffner2023() + imf_obj = imf.Kroupa_2001(multiplicity=multi) + cluster = synthetic.ResolvedCluster(iso, imf_obj, Mcl) + +Resolved (same MF/CSF/q, plus mass-dependent separations):: + + multi = multiplicity.MultiplicityResolvedOffner2023() + +The generic helpers :class:`~imf.multiplicity.MultiplicityLogistic` and +:class:`~imf.multiplicity.MultiplicityPiecewisePowerLaw` are available +for other surveys. Offner does **not** evaluate MF/CSF, :math:`\gamma`, +or separations as a piecewise interpolation of Table 1/2 knots. + +Multiplicity and companion-star fractions +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +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: + +.. math:: + + f(M) = A + \frac{B - A}{1 + (M / M_0)^{-k}} + +with :math:`(A, B, M_0, k) = (0.14, 0.99, 1.41, 1.25)` for MF and +:math:`(0.12, 2.35, 3.57, 0.96)` for CSF. The curve is C-infinity +smooth and saturates at :math:`B \approx 1` for MF. MF is clipped to +:math:`[0, 1]`. CSF is clipped to :math:`[0, \mathrm{CSF_{max}}]`, +raised to at least MF, and forced equal to MF for +:math:`M \le 0.08\,M_\odot` (binaries only). + +SPISEA v2.5 uses :math:`\mathrm{MF} = 0.44\,M^{0.51}` (clipped +to 1) for arrays, plus a scalar-only brown-dwarf staircase +(0 / 8% / 16%). Cluster generation on that class still uses the +stellar power law for brown-dwarf primaries. + +.. 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. + Left: brown-dwarf zoom. Right: BD through early B. Solid: Offner + logistic in log-mass. Dashed: SPISEA v2.5 :math:`0.44\,M^{0.51}`. + Dotted: SPISEA v2.5 scalar BD staircase. Blue points: Offner + Table 1. Shaded: :math:`M \le 0.08\,M_\odot`. + +Mass-ratio index +~~~~~~~~~~~~~~~~ +Companion mass ratios follow :math:`P(q) \propto q^{\gamma}` on +:math:`q_{\min} \le q \le 1` (default :math:`q_{\min} = 0.01`). +:math:`\gamma(M)` is an error-weighted logistic in log-mass fitted to +Table 1 :math:`\gamma_{\mathrm{trunc}}` (1–100 au): + +.. math:: + + \gamma(M) = A + \frac{B - A}{1 + (M / M_0)^{-k}} + +with :math:`(A, B, M_0, k) = (6.6, -1.77, 0.0651, 0.629)`. Call +``q_power_at_mass(mass)`` or ``random_q(x, mass=...)``. Without +``mass``, ``random_q(x)`` keeps the historical stellar-only power law. + +The err-weighted fit undershoots Fontanive et al. (2018) +:math:`8\pm 6\%` MF and :math:`\gamma = 4.8\pm 2.2`: at +:math:`0.033\,M_\odot`, :math:`\gamma \approx 3.3`. That is the +fit, not a bug. SPISEA v2.5 is a step: :math:`\gamma = 6.1` for +:math:`M \le 0.08\,M_\odot` (Fontanive) and :math:`-0.4` above. + +.. 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. Solid: Offner error-weighted logistic. Dashed: + SPISEA v2.5 step (6.1 below :math:`0.08\,M_\odot`, :math:`-0.4` + above). Blue points: Table 1 :math:`\gamma_{\mathrm{trunc}}`. + +.. 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` on :math:`[0.01, 1]` implied by + :math:`P(q)\propto q^{\gamma}`. Offner brown-dwarf companions + are more equal-mass than SPISEA v2.5 stellar :math:`q`. + +Characteristic separation +~~~~~~~~~~~~~~~~~~~~~~~~~ +The characteristic :math:`\mu(a)` is a smooth broken power law in +:math:`\log_{10} a` vs :math:`\log_{10} M`, FGK-pulled, with smoothing +scale :math:`s = 0.1` dex. It is C-infinity (stable +:math:`\log\cosh`; not :math:`\log(\cosh x)` and not a hard +``where`` break): + +.. math:: + + v = \log_{10}(M / M_p),\quad + y_p = \log_{10}(\mu_p) + +.. math:: + + \log_{10} a = y_p + \tfrac{1}{2}(\alpha_L+\alpha_R)\,v + + \tfrac{1}{2}(\alpha_R-\alpha_L)\,s\,\log\cosh(v/s) + +with :math:`\mu_p = 44.46` AU, :math:`M_p = 0.819\,M_\odot`, +:math:`\alpha_L = 1.005`, :math:`\alpha_R = -0.308`, :math:`s = 0.10`. +Linear-space :math:`a` is clipped above 0.1 AU. The implementation +uses the stable form +:math:`\log\cosh x = |x| + \log(1+e^{-2|x|}) - \log 2`. +``log_a_mean(mass)`` returns :math:`\log_{10}(a/\mathrm{AU})`; +``a_mean(mass)`` returns :math:`a` in AU. + +The SPISEA v2.5 :class:`~imf.multiplicity.MultiplicityResolvedDK` +uses a Duchêne & Kraus (2013) broken power law in :math:`a` with a +brown-dwarf blend. That law is not meant for the BD regime. + +.. 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. Solid: Offner smooth broken power law. Dashed: + SPISEA v2.5 mean :math:`a`. Open circles: Table 1 + :math:`\tilde{a}_{\mathrm{all}}`. Filled squares: Table 2 + :math:`\mu` knots (4, 25, 40 AU). Offner BD binaries peak at a + few AU. + +Separation scatter +~~~~~~~~~~~~~~~~~~ +:math:`\sigma(\log_{10} a)` is a 2-parameter logistic pinned at +0.7 / 1.5: + +.. math:: + + \sigma(M) = 0.7 + \frac{0.8}{1 + (M / 0.354)^{-6.05}} + +i.e. :math:`(A, B, M_0, k) = (0.7, 1.5, 0.354, 6.05)`, clipped to +:math:`\ge 0.1`. Call ``sigma_log_a(mass)``. SPISEA v2.5 DK is a +linear fit in :math:`\log M` that saturates above +:math:`2.9\,M_\odot`; the dip near :math:`0.08\,M_\odot` is the +BD/stellar blend. + +.. 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. Solid: Offner 2-parameter logistic. Dashed: + SPISEA v2.5 DK. Filled squares: Table 2 knots (0.7, 1.3, 1.5). + +Resolved draws +~~~~~~~~~~~~~~ +:class:`~imf.multiplicity.MultiplicityResolvedOffner2023` draws +:math:`\log_{10}(a/\mathrm{AU})` from a truncated lognormal with +``loc = log_a_mean(mass)`` and ``scale = sigma_log_a(mass)``, +truncated to 0.01–2000 AU (same limits as +:class:`~imf.multiplicity.MultiplicityResolvedDK`). Eccentricity +and Keplerian angles still follow Duchêne & Kraus (2013) +(:math:`f(e)=2e`, random inclination and angles). + +Brown-dwarf policy and Table 1 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Primaries at or below :math:`0.08\,M_\odot` are binaries only: +CSF = MF and the companion count is capped at 1. + +Table 1 FGKM MF/CF exclude brown-dwarf companions +(:math:`M_{\mathrm{comp}} > 0.075\,M_\odot` for FGKM; +OBA use :math:`q > 0.1`). SPISEA still draws companions down to +``q_min`` (default 0.01), so some stellar primaries get BD +secondaries (~4% for solar-type). Do not read the simulated +stellar-primary MF as a stellar-companion-only statistic. + +Reproducing the figures +~~~~~~~~~~~~~~~~~~~~~~~ +The comparison figures are generated from the multiplicity object +methods (``multiplicity_fraction``, ``q_power_at_mass``, +``a_mean``, ``sigma_log_a``) so they cannot drift from the code. +From the repository root:: + + python docs/figures/plot_mf_offner_vs_spisea2.5.py + python docs/figures/plot_q_sep_offner_vs_spisea2.5.py + Resolved Multiplicity Classes ------------------------------------------ .. autoclass:: imf.multiplicity.MultiplicityResolvedDK :show-inheritance: + +.. autoclass:: imf.multiplicity.MultiplicityResolvedOffner2023 + :show-inheritance: + :members: log_semimajoraxis, log_a_mean, a_mean, sigma_log_a + diff --git a/spisea/imf/imf.py b/spisea/imf/imf.py index 241e04d..fbc3a44 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 4258cb4..198dc8b 100755 --- a/spisea/imf/multiplicity.py +++ b/spisea/imf/multiplicity.py @@ -13,8 +13,108 @@ 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 + +# Fontanive et al. (2018) mass-ratio power-law index used for BD primaries +# in the SPISEA v2.5 MultiplicityUnresolved implementation. +FONTANIVE2018_BD_Q_POWER = 6.1 + # Eventually we should add in separation properties. (a_mean, a_sigma) +# Equal-weight logistic-in-log-mass fit to Offner et al. 2023 Table 1 +# geom-mean (M, MF) and (M, CF) points: +# y(M) = A + (B - A) / (1 + (M / M0)**(-k)) +OFFNER2023_MF_A = 0.14 +OFFNER2023_MF_B = 0.99 +OFFNER2023_MF_M0 = 1.41 +OFFNER2023_MF_K = 1.25 +OFFNER2023_CSF_A = 0.12 +OFFNER2023_CSF_B = 2.35 +OFFNER2023_CSF_M0 = 3.57 +OFFNER2023_CSF_K = 0.96 + + +# Error-weighted logistic in log-mass for Table 1 γ_trunc (1–100 au). +# The model is this logistic, not interpolation of the arrays below. +OFFNER2023_Q_A = 6.6 +OFFNER2023_Q_B = -1.77 +OFFNER2023_Q_M0 = 0.0651 +OFFNER2023_Q_K = 0.629 + +# Smooth broken power law in log10(a) vs log10(M), s=0.1 dex, FGK-pulled. +OFFNER2023_A_MUP = 44.46 +OFFNER2023_A_MP = 0.819 +OFFNER2023_A_ALPHAL = 1.005 +OFFNER2023_A_ALPHAR = -0.308 +OFFNER2023_A_S = 0.10 +OFFNER2023_A_MIN = 0.1 + +# 2-parameter logistic for σ(log10 a); floors/ceilings pinned at 0.7 / 1.5. +OFFNER2023_SIG_A = 0.7 +OFFNER2023_SIG_B = 1.5 +OFFNER2023_SIG_M0 = 0.354 +OFFNER2023_SIG_K = 6.05 + + +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): + """ + Generate random incliniation and angles of binary system + + 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. + + 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 + """ + 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 + + Omega = 360*y + omega = 360*z + + return inclination, Omega, omega + class MultiplicityUnresolved(object): """ The properties of stellar companions (see notes below). @@ -104,12 +204,18 @@ class MultiplicityUnresolved(object): companion_max : bool, optional Sets CSF_max is the max as the max number of companions. Default False. + + binary_only_mass_max : float, optional + Primary mass in solar masses (Msun) at and below which systems + are restricted to at most one companion (CSF = MF). Default is + 0.08 Msun. """ 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, + binary_only_mass_max=H_BURNING_MASS): self.MF_amp = MF_amp self.MF_pow = MF_power @@ -119,6 +225,7 @@ def __init__(self, self.q_pow = q_power self.q_min = q_min self.companion_max = companion_max + self.binary_only_mass_max = binary_only_mass_max def multiplicity_fraction(self, mass): """ @@ -130,14 +237,16 @@ def multiplicity_fraction(self, mass): Parameters ---------- - mass : float or numpy array - Mass of primary star. + mass : float or array_like + Primary mass must be positive, in 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 @@ -165,14 +274,16 @@ def companion_star_fraction(self, mass): Parameters ---------- - mass : float or numpy array - Mass of primary star + mass : float or array_like + Primary mass must be positive, 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 @@ -180,37 +291,73 @@ def companion_star_fraction(self, mass): if np.isscalar(csf): if csf > self.CSF_max: csf = self.CSF_max - if (mass <= 0.08): + if (mass <= self.binary_only_mass_max): csf = self.multiplicity_fraction(mass) else: csf[csf > self.CSF_max] = self.CSF_max - bd = mass <= 0.08 + bd = mass <= self.binary_only_mass_max 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 <= binary_only_mass_max) use + gamma = 6.1 from Fontanive et al. (2018), matching the + companion-mass draw previously special-cased in ``imf.calc_multi``. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, 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 <= self.binary_only_mass_max] = FONTANIVE2018_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 must be positive, 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,23 +365,210 @@ def random_is_multiple(self, x, MF): """ return x < MF - def random_companion_count(self, x, CSF, MF): - """ - Helper function: calculate number of companions. + def random_companion_count(self, x, CSF, MF, mass=None, rng=None): """ - # 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 must be positive, in solar masses (Msun). If given, primaries + at or below ``binary_only_mass_max`` are limited to one + companion. Cluster generation always passes mass so + subclasses can override the BD companion-count policy here. + 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 - + if (mass is not None) and (np.asarray(mass, dtype=float).reshape(-1)[0] <= self.binary_only_mass_max): + n_comp = min(int(n_comp), 1) + 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 mass is not None: + mass = np.atleast_1d(np.asarray(mass, dtype=float)) + bd = mass <= self.binary_only_mass_max + n_comp[bd & (n_comp > 1)] = 1 + + 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` + with ``mass`` so BD companion-count policy lives on this object. + + 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 <= self.binary_only_mass_max + 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 @@ -270,6 +604,8 @@ 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 + + return def log_semimajoraxis(self, mass): """ @@ -328,55 +664,907 @@ def log_semimajoraxis(self, mass): log_semimajoraxis = truncnorm.rvs(a_lower_std, a_upper_std, loc=log_a_mean, scale=log_a_std) return log_semimajoraxis - def random_e(self, x): + +class MultiplicityPiecewisePowerLaw(MultiplicityUnresolved): + """ + Multiplicity described by a piecewise power law in primary mass. + + 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 forced + equal to MF for primaries at or below ``binary_only_mass_max`` + (binaries only). CSF is also raised to at least MF so the Poisson + companion-count draw is well defined. + + 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. Default False. + binary_only_mass_max : float, optional + Primary mass in solar masses (Msun) at and below which systems + are binaries only. Default 0.08 Msun. + """ + 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, + binary_only_mass_max=H_BURNING_MASS): + 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, + binary_only_mass_max=binary_only_mass_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 must be positive, 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], raised to at least MF, and set equal + to MF for primaries at or below ``binary_only_mass_max``. + 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 must be positive, 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) + bd = mass_arr <= self.binary_only_mass_max + csf[bd] = mf[bd] + if return_scalar: + return float(csf[0]) + return csf + + +class MultiplicityLogistic(MultiplicityUnresolved): + """ + Multiplicity described by a logistic in log primary mass. + + f(M) = A + (B - A) / (1 + (M / M0)**(-k)) + + As M → 0, f → A; as M → ∞, f → B. The same functional form is + used for MF and CSF with independent coefficients. + ``MultiplicityUnresolvedOffner2023`` uses this class with + coefficients fitted to Offner et al. (2023) Table 1. + + MF is clipped to [0, 1]. CSF is clipped to [0, CSF_max], raised + to at least MF, and forced equal to MF for primaries at or below + ``binary_only_mass_max`` (binaries only). + + 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. Default False. + binary_only_mass_max : float, optional + Primary mass in solar masses (Msun) at and below which systems + are binaries only. Default 0.08 Msun. + """ + 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, + binary_only_mass_max=H_BURNING_MASS): + 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, + binary_only_mass_max=binary_only_mass_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 must be positive, 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], raised to at least MF, and set equal + to MF for primaries at or below ``binary_only_mass_max``. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, 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) + + # Fix all brown dwarf binaries so they only have one companion + bd = mass_arr <= self.binary_only_mass_max + csf[bd] = mf[bd] + + if return_scalar: + return float(csf[0]) + + return csf + + +class MultiplicityUnresolvedOffner2023(MultiplicityLogistic): + """ + Unresolved multiplicity from Offner et al. 2023 Table 1, including + brown dwarfs. + + 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 + `_. + + The multiplicity fraction and companion frequency 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), 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%), which is consistent with the + Burgasser/Close BD points and within ~1–2σ of Fontanive. MF is + clipped to [0, 1]. Below 0.08 Msun, CSF = MF (binaries only; + THF is tiny). + + 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 draws use 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). 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). + + Characteristic separation μ(a) is a **smooth broken power law** + in log10(a) vs log10(M) (s = 0.1 dex, FGK-pulled), C-infinity + via a stable logcosh. σ(log10 a) is a 2-parameter logistic + pinned at 0.7 / 1.5. See :meth:`log_a_mean` and + :meth:`sigma_log_a`. Resolved draws use those as loc / scale + of a truncated lognormal. + + This class is opt-in; it does not change the SPISEA v2.5 + :class:`MultiplicityUnresolved` default. + + Parameters + ---------- + 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. Default False. + binary_only_mass_max : float, optional + Primary mass in solar masses (Msun) at and below which systems + are binaries only (CSF = MF, at most one companion). Default + 0.08 Msun. + """ + def __init__(self, CSF_max=3, q_power=0.2, q_min=0.01, + companion_max=False, binary_only_mass_max=H_BURNING_MASS): + + super(MultiplicityUnresolvedOffner2023, self).__init__( + MF_A=OFFNER2023_MF_A, + MF_B=OFFNER2023_MF_B, + MF_M0=OFFNER2023_MF_M0, + MF_k=OFFNER2023_MF_K, + CSF_A=OFFNER2023_CSF_A, + CSF_B=OFFNER2023_CSF_B, + CSF_M0=OFFNER2023_CSF_M0, + CSF_k=OFFNER2023_CSF_K, + CSF_max=CSF_max, q_power=q_power, q_min=q_min, + companion_max=companion_max, + binary_only_mass_max=binary_only_mass_max) + + # Table 1/2 data that was fit; evaluation uses a smooth function. + self.q_mass = np.array(OFFNER2023_Q_MASS, dtype=float) + self.q_gamma = np.array(OFFNER2023_Q_GAMMA, dtype=float) + + return + + def q_power_at_mass(self, mass): + """ + Mass-ratio power-law index γ(M), P(q) ∝ q^γ. + + Error-weighted logistic in log-mass fitted to Table 1 + γ_trunc. Not an interpolation of the Table 1 knots. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, 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, OFFNER2023_Q_A, OFFNER2023_Q_B, OFFNER2023_Q_M0, OFFNER2023_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). + Linear-space a is clipped to 0.1 AU. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, 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, OFFNER2023_A_MUP, OFFNER2023_A_MP, + OFFNER2023_A_ALPHAL, OFFNER2023_A_ALPHAR, OFFNER2023_A_S, + a_min=OFFNER2023_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 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): + 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. + + 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. + """ + + # Calculate the standard deviation of log10(a / 1 AU) using a 2-parameter logistic + sigma_log_a = _logistic_in_logm(mass, OFFNER2023_SIG_A, OFFNER2023_SIG_B, + OFFNER2023_SIG_M0, OFFNER2023_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): + """ + Resolved Offner et al. 2023 multiplicity: Table 1 MF/CF plus + mass-dependent separations. + + Separations are drawn from a truncated lognormal in log10(a/AU) + with loc = :meth:`log_a_mean` (smooth broken power law, s = 0.1 + dex, FGK-pulled) and scale = :meth:`sigma_log_a` (2-parameter + logistic pinned at 0.7 / 1.5). Brown-dwarf binaries peak near a + few AU. Truncation is 0.01–2000 AU, same as + :class:`MultiplicityResolvedDK`. + + Eccentricity and Keplerian angles follow Duchêne & Kraus (2013), + same as :class:`MultiplicityResolvedDK`. + + Opt-in; does not replace :class:`MultiplicityResolvedDK`. + + Parameters + ---------- + 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. Default False. + binary_only_mass_max : float, optional + Primary mass in solar masses (Msun) at and below which systems + are binaries only. Default 0.08 Msun. + """ + def __init__(self, **kwargs): + super(MultiplicityResolvedOffner2023, self).__init__(**kwargs) + # Table 1/2 data that was fit; draws use log_a_mean / sigma_log_a. + self.sep_mass = np.array(OFFNER2023_SEP_MASS, dtype=float) + self.sep_mu_au = np.array(OFFNER2023_SEP_MU_AU, dtype=float) + self.sep_sig_mass = np.array(OFFNER2023_SEP_SIG_MASS, dtype=float) + self.sep_sig = np.array(OFFNER2023_SEP_SIG, dtype=float) + + return + + def log_semimajoraxis(self, mass): + """ + Draw log10(a/AU) from a mass-dependent truncated lognormal. + + Parameters + ---------- + mass : float or array_like + Primary mass must be positive, in solar masses (Msun). + + 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. + """ + 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) + + return log_a + +# Convenience alias; unresolved Table 1 model is the usual opt-in object. +MultiplicityOffner2023 = MultiplicityUnresolvedOffner2023 + + +def _two_point_powerlaw(mass_1, y_1, mass_2, y_2): + """ + Amplitude and power for y = A * M**alpha through two (M, y) points. + + Parameters + ---------- + mass_1, mass_2 : float + Primary masses of the two anchor points, in solar masses + (Msun). Must be positive and distinct. + y_1, y_2 : float + Ordinate values at ``mass_1`` and ``mass_2``. Units match the + fitted quantity (dimensionless for MF/γ, mean companion count + for CSF, AU for characteristic a, dex for σ). + + Returns + ------- + amp : float + Power-law amplitude A, in units of y / Msun**alpha. + power : float + Power-law index alpha, dimensionless. + """ + power = np.log(y_2 / y_1) / np.log(mass_2 / mass_1) + amp = y_1 / (mass_1 ** power) + + return amp, power + + +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 must be positive, 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 must be positive, 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 must be positive, 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 + + +def _offner2023_table1_geom_mass(m_lo, m_hi): + """ + Geometric-mean primary mass of a Table 1 M1 interval. + + Parameters + ---------- + m_lo, m_hi : float + Low and high edges of the Table 1 primary-mass interval, in + solar masses (Msun). Must be positive. + + Returns + ------- + mass : float + Geometric mean sqrt(m_lo * m_hi) in solar masses (Msun). + """ + # Calculate the geometric mean of the two masses + return float(np.sqrt(m_lo * m_hi)) + +# Table 1/2 data that was fit; the model is a smooth function, +# not interpolation of these arrays. +OFFNER2023_Q_MASS = np.array([ + _offner2023_table1_geom_mass(0.019, 0.058), # Fontanive+2018: 4.8 + 0.065, # L/early-T: 2-3 (text) + _offner2023_table1_geom_mass(0.080, 0.095), # Close+2003: 3.3 + _offner2023_table1_geom_mass(0.06, 0.15), # Allen+2007: 1.7 + _offner2023_table1_geom_mass(0.15, 0.30), # Winters mid-M: 0.7 + _offner2023_table1_geom_mass(0.3, 0.6), # Winters early-M: 0.1 + _offner2023_table1_geom_mass(0.75, 1.25), # Raghavan FGK: 0.2 + _offner2023_table1_geom_mass(1.6, 2.4), # De Rosa A: -1.3 + _offner2023_table1_geom_mass(3.0, 5.0), # MDS 3-5: -1.0 + _offner2023_table1_geom_mass(5.0, 8.0), # MDS 5-8: -1.7 + _offner2023_table1_geom_mass(8.0, 17.0), # MDS 8-17: -1.6 + _offner2023_table1_geom_mass(17.0, 50.0), # Sana O: -1.4 +]) +OFFNER2023_Q_GAMMA = np.array([ + 4.8, 2.5, 3.3, 1.7, 0.7, 0.1, 0.2, -1.3, -1.0, -1.7, -1.6, -1.4 +]) + +# Table 1 ã_all (au) and Table 2 lognormal μ (au) vs geom-mean M1. +# Table 2 μ is listed where both exist (late-M, early-M, FGK). +# Fit data only; evaluation uses the smooth broken power law. +OFFNER2023_SEP_MASS = np.array([ + _offner2023_table1_geom_mass(0.019, 0.058), # Fontanive ã_all=2.9 + _offner2023_table1_geom_mass(0.080, 0.095), # Close ã_all=3.7 + _offner2023_table1_geom_mass(0.075, 0.15), # Table 2 late-M μ=4 + _offner2023_table1_geom_mass(0.15, 0.30), # Winters mid-M ã_all=10 + _offner2023_table1_geom_mass(0.3, 0.6), # Table 2 early-M μ=25 + _offner2023_table1_geom_mass(0.75, 1.25), # Table 2 FGK μ=40 + _offner2023_table1_geom_mass(1.6, 2.4), # Moe & Kratter ã_all=32 + _offner2023_table1_geom_mass(3.0, 5.0), # MDS ã_all=28 + _offner2023_table1_geom_mass(5.0, 8.0), # MDS ã_all=25 + _offner2023_table1_geom_mass(8.0, 17.0), # MDS ã_all=23 + _offner2023_table1_geom_mass(17.0, 50.0), # Sana ã_all=19 +]) +OFFNER2023_SEP_MU_AU = np.array([ + 2.9, 3.7, 4.0, 10.0, 25.0, 40.0, 32.0, 28.0, 25.0, 23.0, 19.0 +]) +# Table 2 σ_log a at the three published bins. Fit data only; +# evaluation uses the 2-parameter logistic. +OFFNER2023_SEP_SIG_MASS = np.array([ + _offner2023_table1_geom_mass(0.075, 0.15), + _offner2023_table1_geom_mass(0.3, 0.6), + _offner2023_table1_geom_mass(0.75, 1.25), +]) +OFFNER2023_SEP_SIG = np.array([0.7, 1.3, 1.5]) diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 0fe624a..bc80be5 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -532,10 +532,15 @@ 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( + # 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')): + companions.add_column(Column(multi_props.log_semimajoraxis(star_systems['mass'][companions['system_idx']]), 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) diff --git a/spisea/tests/test_imf.py b/spisea/tests/test_imf.py index 943dc0a..984125f 100755 --- a/spisea/tests/test_imf.py +++ b/spisea/tests/test_imf.py @@ -30,6 +30,42 @@ def test_generate_cluster(): return + +def test_generate_cluster_offner2023(): + """ + generate_cluster with Offner et al. 2023 multiplicity: vectorized MF, + BD primaries have at most one companion, 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 + n_comp = np.sum(~compMass.mask, axis=1) + bd = mass <= 0.08 + assert np.all(n_comp[bd] <= 1) + 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 5322665..f0d6782 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 @@ -243,3 +244,339 @@ 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], + binary_only_mass_max=0.05) + 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 sets BD CSF = MF.""" + ml = multiplicity.MultiplicityLogistic( + MF_A=0.1, MF_B=1.5, MF_M0=1.0, MF_k=2.0, + CSF_A=0.05, CSF_B=4.0, CSF_M0=2.0, CSF_k=1.0, + CSF_max=2.0, binary_only_mass_max=0.08) + # 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) + # BD CSF = 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) + + +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, 0.14, 0.99, 1.41, 1.25, 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 matches Table 1 CF for stellar primaries; CSF = MF for BDs.""" + 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 m <= multiplicity.H_BURNING_MASS: + assert np.isclose(csf, mf, atol=1e-12), \ + '{0}: BD CSF should equal MF'.format(name) + elif mlo >= 1.6: + # Logistic CF tracks A/B; Moe & Kratter residual ~0.1 is ok + tol = max(0.12, 0.12 * 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) + 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_offner2023_bd_binaries_only(): + """BD primaries have CSF = MF and companion counts of 0 or 1.""" + multi = multiplicity.MultiplicityUnresolvedOffner2023() + rng = np.random.default_rng(123) + masses = np.array([0.02, 0.04, 0.07, 0.08]) + mf = multi.multiplicity_fraction(masses) + csf = multi.companion_star_fraction(masses) + np.testing.assert_allclose(csf, mf) + # Force multiples so we test the count draw, not the MF coin flip. + n_comp = multi.draw_n_companions(masses, csf, mf, rng) + assert np.all(n_comp <= 1) + assert np.all(n_comp >= 1) + + # Full companion-mass assignment: never more than one companion column + is_mult = np.ones(len(masses), dtype=bool) + comp, sys_mass, is_mult_out = multi.draw_companion_masses( + masses, is_mult, csf, mf, rng, mass_min=0.01) + assert comp.shape[1] == 1 + assert np.all(np.sum(~comp.mask, axis=1) <= 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, 6.6, -1.77, 0.0651, 0.629)) + np.testing.assert_allclose( + multi.sigma_log_a(m), + multiplicity._logistic_in_logm( + m, 0.7, 1.5, 0.354, 6.05, clip_min=0.1)) + np.testing.assert_allclose( + multi.log_a_mean(m), + multiplicity._smooth_broken_loglog( + m, 44.46, 0.819, 1.005, -0.308, 0.10, a_min=0.1)) + # 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, 6.6, -1.77, 0.0651, 0.629)) + 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() + np.random.seed(0) + log_a = multi.log_semimajoraxis(np.full(5000, 0.04)) + 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)) + med_a_s = 10 ** np.median(log_a_s) + assert med_a_s > 10.0 + assert med_a_s > med_a + + +def test_offner2023_alias_and_resolved_methods(): + """Public names and resolved orbital methods exist.""" + assert multiplicity.MultiplicityOffner2023 is \ + multiplicity.MultiplicityUnresolvedOffner2023 + resolved = multiplicity.MultiplicityResolvedOffner2023() + assert hasattr(resolved, 'log_semimajoraxis') + assert hasattr(resolved, 'log_a_mean') + assert hasattr(resolved, 'sigma_log_a') + 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) + # 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_offner_generate_cluster_companions(): + """IMF cluster generation with Offner multiplicity produces BD binaries only.""" + 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) + bd = mass <= 0.08 + n_comp = np.sum(~comp_mass.mask, axis=1) + assert np.all(n_comp[bd] <= 1) + 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 or the BD + companion cap; those policies live 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 + + 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) +