<jupyter_start><jupyter_text>Let's test that it works by choosing a random SNR..<jupyter_code>snr = 1000

asinis = []
lc = SyntheticBinary(time, freqs, amps, period, eccen, asini, varpi, tref, tau=False)
s_rms = lc.amplitude.max() / (np.sqrt(np.pi / len(lc.time)) * snr)
for i in tqdm(range(5000)):
    noisy_flux = lc.flux_true + np.random.normal(loc=0.0, scale=s_rms, size=len(lc.time))
    x0 = np.array([freqs[0], amps[0], 0]).flatten()
    popt, pcov = curve_fit(model, lc.time, noisy_flux, p0=x0)
    asinis.append(popt[-1])
plt.hist(asinis);<jupyter_output><empty_output><jupyter_text>Now let's do it for a range of SNRs. The vals below are the d Sct SNR percentiles obtained in previous notebook<jupyter_code>vals = [14.6757, 1.4668, 0.72857]
from tqdm import tqdm

np.random.seed(42)
snrs = np.geomspace(10,1000,50)

asinis_sd = []
for snr in tqdm(snrs):
    asinis = []
    lc = SyntheticBinary(time, freqs, amps, period, eccen, asini, varpi, tref, tau=False)
    s_rms = lc.amplitude.max() / (np.sqrt(np.pi / len(lc.time)) * snr)
    for i in range(1000):
        noisy_flux = lc.flux_true + np.random.normal(loc=0.0, scale=s_rms, size=len(lc.time))
        x0 = np.array([freqs[0], amps[0], 0]).flatten()
        popt, pcov = curve_fit(model, lc.time, noisy_flux, p0=x0)
        asinis.append(popt[-1])
    asinis_sd.append(np.std(asinis))
from scipy import stats
import astropy.units as u

slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(snrs), np.log(asinis_sd))
vals = [14.6757, 1.4668, 0.72857]
fig, axes = plt.subplots(1,2,
                         figsize=[7,2.5], 
                         constrained_layout=True)

linewidth=0.8
ax = axes[0]
ax.plot(snrs, asinis_sd, '.k', markersize=3, linewidth=0.7)
# ax.plot(snrs, np.exp((slope*np.log(snrs) + intercept)), '-.', c='black', linewidth=0.7)

ax.fill_between(snrs, 10**2,np.exp((slope*np.log(snrs) + intercept)), alpha=0.2)
ax.set_xscale('log')
ax.set_yscale('log')


ax.plot([0,92.1],[8,8],c=red, linewidth=linewidth)
ax.plot([0,358.4],[2,2], c=blue, linewidth=linewidth)
ax.plot([0,749.4],[0.96680069,0.96680069], c=green, linewidth=linewidth)

ax.plot([92.1,92.1],[0, 8],  c=red, linewidth=linewidth)
ax.plot([358.4,358.4],[0,2],  c=blue, linewidth=linewidth)
ax.plot([749.4,749.4],[0,0.96680069], c=green, linewidth=linewidth)

ax.set_xlabel('SNR of primary frequency')
ax.set_ylabel(r'$a\sin{i}/c$ scatter (s)')

ax.set_xlim(snrs[0], snrs[-1])
ax.set_ylim(None,10**2)


ax = axes[1]
periods = np.geomspace(10,1470,500)
stellar_mass = 1.8

worst_case = (mass_function(periods*u.day, 8*u.s).to(u.M_jup))**(1/3) * ((stellar_mass*u.M_sun).to(u.M_jup))**(2/3)
typical_case = (mass_function(periods*u.day, 2*u.s).to(u.M_jup))**(1/3) * ((stellar_mass*u.M_sun).to(u.M_jup))**(2/3)
best_case = (mass_function(periods*u.day, 1*u.s).to(u.M_jup))**(1/3) * ((stellar_mass*u.M_sun).to(u.M_jup))**(2/3)

ax.plot(periods, worst_case, label='SNR=50', c=red, linewidth=linewidth)
ax.plot(periods, typical_case, label='SNR=500', c=blue, linewidth=linewidth)
ax.plot(periods, best_case, label='SNR=1000', c=green, linewidth=linewidth)

ax.fill_between(periods, 75, 13, alpha=0.2, color='orange')
ax.set_xlabel(r'P$_{\rm orb}$ (d)')
ax.set_ylabel(r'Companion mass (M$_{\rm Jup}$)')

ax.set_xscale('log')
ax.set_yscale('log')

ax.set_yticks([1,10,100])
ax.set_xlim(10,1470)

per_med = len(periods) // 2
ax.text(periods[per_med], worst_case[per_med].value, '8 s', rotation=-18, c=red)
ax.text(periods[per_med], typical_case[per_med].value-0.5, '2 s', rotation=-18, c=blue)
ax.text(periods[per_med], best_case[per_med].value-0.5, '1 s', rotation=-18, c=green)

axes[0].annotate('Detectable', (0.65,0.7), xycoords='axes fraction')
plt.savefig(overleaf_path + 'best_mass_range.pdf', dpi=300, bbox_inches='tight', pad_inches=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2006793833289783, 'fraction_numerical': 0.06454141625293965, 'mean_word_length': 5.2042139384116695, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '488203', 'license': 'MIT', 'n_tokens_mistral': 1851, 'n_tokens_neox': 1685, 'n_words': 271, 'path': 'paper/notebooks/(5.1) SNR vs asini relations.ipynb', 'repo_name': 'danielhey/Maelstrom', 'size': '3827'}