In [1]:
import numpy as np
import matplotlib.pyplot as plt
import emcee
import corner
%matplotlib inline
In [2]:
# the correlation matrix from Gaia, and the vector of best-fit parameters.
# these can be retrieved from the Gaia archive; I just copy them here so
# this notebook can be run without an internet connection or TAP instance.
corr_vec_NS1 = np.array([0.62615186, -0.84156424, -0.75608426, -0.18818806, -0.04869297, 0.18428592,
0.65590525, 0.27018258, -0.40238973, -0.44358525, 0.8284085, 0.48271143,
-0.72155017, -0.0024235346, 0.30181977, 0.98108965, 0.5335061, -0.7680332,
-0.21360154, 0.74536324, 0.7958441, -0.9792101, -0.5300254, 0.7683197,
0.15293673, -0.7035678, -0.8137645, -0.9937877, -0.9042292, -0.35511398,
0.66397995, 0.32691798, -0.77383995, -0.7256531, -0.94355947, 0.92267823,
0.5648844, 0.9491089, -0.6664481, 0.052533828, 0.23379754, 0.39255613, 0.4690493,
-0.4670044, -0.26245797, -0.55593956, -0.10765043, 0.19196706, 0.17331037,
-0.8728608, -0.25485885, -0.6702889, 0.65313816, 0.6671808, -0.11448689, 0.9874635,
0.55741996, -0.794174, -0.19639102, 0.7141803, 0.8134602, 0.9979075, -0.99479944,
-0.9351806, 0.49180174, -0.63134485])
mu_vec_NS1 = np.array([ 2.18086202e+02, -1.03663474e+01, 1.23158367e+00, -2.28644068e-02,
-9.23538943e+01, 7.29328579e-01, 2.04250825e-01, 1.48166628e-01,
2.05054356e+00, 1.35366326e-01, 7.36017850e+02, 2.28220981e+02])
# ra, dec, parallax, pmra, pmdec, A, B, F, G, ecc, period, t_peri
# we're going to inflate the astrometric uncertainties by a factor of 3
# we also need to convert mas to degrees
mas_to_deg = 3600000
err_vec_NS1 = 3*np.array([0.12653045/mas_to_deg, 0.04505481/mas_to_deg, 0.040705957, 0.024724633, 0.056686666,
0.039688792, 0.69762796, 0.23487046, 0.07399644, 0.035554353, 11.611804, 37.795406])
Nparam = 12
triangle_nums = (np.arange(1, Nparam)*(np.arange(1, Nparam) - 1)//2)
cov_mat_NS1 = []
count = 0
for i, sigx in enumerate(err_vec_NS1):
this_row = []
for j, sigy in enumerate(err_vec_NS1):
if i == j:
this_row.append(err_vec_NS1[i]**2)
elif j < i:
this_row.append(cov_mat_NS1[j][i])
else:
this_idx = triangle_nums[j-1] + i
this_row.append(corr_vec_NS1[this_idx]*sigx*sigy)
cov_mat_NS1.append(this_row)
cov_mat_NS1 = np.array(cov_mat_NS1)
cov_mat_inv = np.linalg.inv(cov_mat_NS1)
In [3]:
# let's gather the RVs
jds = np.array([2459795.5592, 2459767.5238, 2459813.4918, 2459831.5038,
2459832.4866, 2459929.0525, 2459981.0076, 2459990.8528,
2459992.8178, 2460012.816 , 2460026.7002, 2460036.8989,
2460039.7097, 2460046.9207, 2460050.6922, 2460059.8472,
2460072.7468, 2460078.79 , 2460085.6962, 2460094.8017,
2460098.6573, 2460109.7473, 2460110.6418, 2460133.6779,
2460139.5731, 2460185.4878, 2460186.4889, 2460306.0444,
2460044.7941, 2460044.8121, 2460054.8803, 2460054.8601,
2460061.9128, 2460061.8931, 2460072.8152, 2460072.8354,
2460086.8885, 2460086.9084, 2460106.8061, 2460122.781 ,
2460122.7615, 2460303.0427, 2459977.8497, 2460037.6527,
2460057.7359, 2460095.6516, 2460132.5799, 2460168.5282,
2460076.611 , 2460038.7825, 2460327.0182, 2460338.8541,
2460340.8396, 2460345.8491, 2460352.9597, 2460374.9193,
2460386.9141, 2460400.8157, 2460406.6343])
rvs = np.array([136.19, 130. , 139.74, 143.61, 143.6 , 154.65, 153.69, 153.22,
153.16, 151.78, 150.29, 149.37, 149.14, 148.36, 147.89, 146.89,
145.13, 144.35, 143.28, 142.6 , 141.94, 140.43, 140.22, 136.95,
136.22, 128.99, 128.99, 113.88, 148.47, 148.7 , 147.34, 147.23,
146.64, 146.48, 145.07, 145.4 , 143.2 , 143.62, 140.97, 138.64,
138.51, 114.39, 154.02, 149.41, 147.01, 142.37, 136.99, 131.85,
144.73, 149.18, 112.34, 111.97, 111.67, 111.44, 111.21, 111.29,
111.63, 112.42, 112.81])
rv_errs = np.array([3. , 3. , 0.1 , 0.41 , 0.25 , 0.5 , 0.085, 0.15 , 0.15 ,
0.12 , 0.15 , 0.068, 0.09 , 0.067, 0.1 , 0.105, 0.2 , 0.16 ,
0.22 , 0.082, 0.249, 0.065, 0.073, 0.066, 0.132, 0.212, 0.148,
0.116, 0.081, 0.093, 0.135, 0.112, 0.082, 0.091, 0.227, 0.13 ,
0.156, 0.149, 0.17 , 0.105, 0.088, 0.321, 0.053, 0.04 , 0.053,
0.032, 0.084, 0.082, 0.04 , 0.022, 0.067, 0.105, 0.084, 0.104,
0.08 , 0.094, 0.093, 0.06 , 0.12])
rv_names = np.array(['MagE', 'MagE', 'FEROS', 'FEROS', 'FEROS', 'PEPSI', 'TRES',
'FEROS', 'FEROS', 'FEROS', 'FEROS', 'TRES', 'FEROS', 'TRES',
'FEROS', 'TRES', 'FEROS', 'TRES', 'FEROS', 'TRES', 'FEROS', 'TRES',
'FEROS', 'TRES', 'FEROS', 'FEROS', 'FEROS', 'TRES', 'APF', 'APF',
'APF', 'APF', 'APF', 'APF', 'APF', 'APF', 'APF', 'APF', 'APF',
'APF', 'APF', 'APF', 'MIKE', 'MIKE', 'MIKE', 'MIKE', 'MIKE',
'MIKE', 'MIKE', 'MIKE', 'TRES', 'FEROS', 'FEROS', 'MIKE', 'TRES',
'TRES', 'TRES', 'FEROS', 'MIKE'])
# prior on the mass of the luminous star
m1_obs, m1_err = 0.79, 0.01
# the reference epoch of the gaia t_periastron; i.e., JD 2016.0
jd_ref = 2457389.0
In [4]:
# functions to predict RVs at a given time
def fsolve_newton(Mi, ecc, xtol = 1e-10):
'''
numerically solve with newton's method.
Mi: float, 2*pi/P*(t - T_p)
ecc: float, eccentricity
xtol: float, tolerance
'''
eps = 1
EE = Mi + ecc*np.sin(Mi) + ecc**2/2*np.sin(2*Mi) + ecc**3/8*(3*np.sin(3*Mi) - np.sin(Mi))
while eps > xtol:
EE0 = EE
EE = EE0 + (Mi + ecc*np.sin(EE0) - EE0)/(1 - ecc*np.cos(EE0))
eps = np.abs(EE0 - EE)
return EE
def solve_kep_eqn(t, T_p, P, ecc):
'''
Solve Keplers equation E - e*sin(E) = M for E
Here, M = 2*pi/P*(t - T_p)
t: array of times
T_p: float; periastron time
P: float, period
ecc: float, eccentricity
'''
M = 2*np.pi/P * (t - T_p)
E = np.zeros(t.shape)
for i,Mi in enumerate(M):
E[i] = fsolve_newton(Mi = Mi, ecc = ecc, xtol = 1e-10)
return E
def get_radial_velocities(t, P, T_p, ecc, K, omega, gamma):
'''
t: array of times
P: float, period
T_p: float, periastron time
ecc: float, eccentricity
K: float, RV semi-amplitude
omega: float, longitude of periastron (radians, not degrees)
gamma: float, center-of-mass RV
'''
E = solve_kep_eqn(t = t, T_p = T_p, P = P, ecc = ecc)
f = 2*np.arctan2(np.sqrt(1+ecc)*np.sin(E*.5),np.sqrt(1-ecc)*np.cos(E*.5))
vr = K*(np.cos(f + omega) + ecc*np.cos(omega)) + gamma
return vr
In [5]:
def get_a0_mas(period, m1, m2, parallax):
'''
so we can avoid astropy units
'''
G = 6.6743e-11
Msun = 1.98840987069805e+30
AU = 1.4959787e+11
a_au = (((period*86400)**2 * G * (m1*Msun + m2*Msun)/(4*np.pi**2))**(1/3.))/AU
a_mas = a_au*parallax
q = m2/m1
a0_mas = a_mas*q/(1+q)
return a0_mas
def get_Kstar_kms(period, inc_deg, m1, m2, ecc):
'''
RV semi-amplitude of the luminous star
period: float, days
inc_deg: float, inclination in degrees
m1: float, mass of luminous star in Msun
m2: float, mass of companion in Msun
ecc: float, eccentricity
'''
G = 6.6743e-11 # SI units
Msun = 1.98840987069805e+30
Kstar = (2*np.pi*G*(m2*Msun) * (m2/(m1 + m2))**2 / (period*86400 * (1 - ecc**2)**(3/2)))**(1/3) * np.sin(inc_deg*np.pi/180)
return Kstar/1000 # km/s
def lnL_with_covariance(y_vec, mu_vec, Sigma_mat, Sigma_mat_inv):
'''
calculates the likelihood of a vector y_vec, given a multivariate Gaussian likelihood characterized by
a mean mu_vec and a covariance matrix Sigma_ma
Sigma_mat_inv is just just the inverse of the covariance matrix (faster to precompute)
'''
return -0.5*np.dot((y_vec - mu_vec).T, np.dot(Sigma_mat_inv, y_vec - mu_vec)) \
-0.5*np.log(np.linalg.det(Sigma_mat)) - len(mu_vec)/2*np.log(2*np.pi)
In [6]:
def ln_likelihood(theta, fit_rvs = False):
'''
theta = ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1
fit_rvs: bool; whether to include the RVs or not
'''
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1 = theta
inc, omega, w = inc_deg*np.pi/180, omega_deg*np.pi/180, w_deg*np.pi/180
if v_com < 0 or v_com > 200 or m1 < 0 or m2 < 0: # broad prior to prevent numerical problems
return -np.inf
if fit_rvs:
Kstar_kms = get_Kstar_kms(period = period, inc_deg = inc_deg, m1 = m1, m2 = m2, ecc = ecc)
rv_pred = get_radial_velocities(t = jds-jd_ref, P = period, T_p = t_peri, ecc = ecc, K = Kstar_kms,
omega = w, gamma = v_com)
lnL = -0.5*np.sum((rv_pred - rvs)**2/rv_errs**2)
else:
lnL = 0
a0_mas = get_a0_mas(period, m1, m2, parallax)
A_pred = a0_mas*( np.cos(w)*np.cos(omega) - np.sin(w)*np.sin(omega)*np.cos(inc) )
B_pred = a0_mas*( np.cos(w)*np.sin(omega) + np.sin(w)*np.cos(omega)*np.cos(inc) )
F_pred = -a0_mas*( np.sin(w)*np.cos(omega) + np.cos(w)*np.sin(omega)*np.cos(inc) )
G_pred = -a0_mas*( np.sin(w)*np.sin(omega) - np.cos(w)*np.cos(omega)*np.cos(inc) )
y_vec = np.array([ra, dec, parallax, pmra, pmdec, A_pred, B_pred, F_pred, G_pred, ecc, period, t_peri])
lnL += lnL_with_covariance(y_vec = y_vec, mu_vec = mu_vec_NS1, Sigma_mat = cov_mat_NS1, Sigma_mat_inv=cov_mat_inv)
lnL += -0.5*(m1 - m1_obs)**2/m1_err**2
if np.isfinite(lnL):
return lnL
else:
return -np.inf
In [7]:
ndim = 14
nwalkers = 64
p0 = np.array([2.18086202e+02, -1.03663474e+01, 1.23, -0.1,
-92.3, 731, 0.124, 68.9, 82.7, 259.6, 187.7, 133.0, 1.90, 0.79])
p0 = np.tile(p0, (nwalkers, 1))
p0 += p0 * 1.0e-10 * np.random.normal(size=(nwalkers, ndim))
sampler_norvs = emcee.EnsembleSampler(nwalkers, ndim, ln_likelihood, args=[False])
# 5000 step burn-in
state = sampler_norvs.run_mcmc(p0, 5000, progress = True)
# 5000 more steps
sampler_norvs.reset()
state = sampler_norvs.run_mcmc(state, 5000, progress = True)
100%|████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [00:05<00:00, 949.07it/s] 100%|████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [00:05<00:00, 957.12it/s]
In [8]:
# next, the "with-rvs" fit. This will take a few minutes because predicting RVs is moderately expensive
# (in production runs we use a compiled C function to speed up the RV prediction)
sampler_withrvs = emcee.EnsembleSampler(nwalkers, ndim, ln_likelihood, args=[True])
# 5000 step burn-in
state = sampler_withrvs.run_mcmc(p0, 5000, progress = True)
# 5000 more steps
sampler_withrvs.reset()
state = sampler_withrvs.run_mcmc(state, 5000, progress = True)
100%|█████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [01:41<00:00, 49.21it/s] 100%|█████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [01:43<00:00, 48.21it/s]
In [9]:
%matplotlib inline
labels = [r'$P_{\rm orb}\,\,[\rm days]$', r'$\rm ecc$', r'$\rm inc\,\,[\rm deg]$', r'$\Omega\,\,[\rm deg]$',
r'$\omega\,\,[\rm deg]$', r'$T_{p}\,\,[\rm days]$', r'$\gamma\,\,[\rm km\,s^{-1}]$',
r'$M_2\,\,[M_{\odot}]$', r'$M_{\star}\,\,[M_{\odot}]$']
# now compare the fits
half_sampler = sampler_withrvs.flatchain.T[5:].T
half_sampler_norvs = sampler_norvs.flatchain.T[5:].T
fig = corner.corner(half_sampler_norvs, labels=labels, show_titles=True, plot_datapoints = False,
plot_density = False, color = 'c')
fig = corner.corner(half_sampler, labels = labels, show_titles=True, plot_datapoints = False,
plot_density = False, color = 'k', fig = fig, label_kwargs = {'fontsize': 20})
fig.axes[4].plot([], [], 'c', label = r'$\rm astrometry\,\,only$')
fig.axes[4].plot([], [], 'k', label = r'$\rm astrometry+RVs$')
fig.axes[4].legend(loc = 'upper right', frameon = False, fontsize=30)
Out[9]:
<matplotlib.legend.Legend at 0x11921b070>
In [10]:
from astropy.time import Time
f = plt.figure(figsize = (12, 11))
ax0 = f.add_axes([0.05, 0.65, 0.9, 0.3])
ax1 = f.add_axes([0.05, 0.35, 0.9, 0.3])
ax2 = f.add_axes([0.05, 0.05, 0.9, 0.3])
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1 = sampler_norvs.flatchain.T
randints = np.random.randint(0, len(period), 50)
t_grid = np.linspace(np.min(jds) - 600, np.max(jds) + 600, 3000)
for i in randints:
Kstar_kms = get_Kstar_kms(period = period[i], inc_deg = inc_deg[i], m1 = m1[i], m2 = m2[i], ecc = ecc[i])
rv_pred = get_radial_velocities(t = t_grid-jd_ref, P = period[i], T_p = t_peri[i], ecc = ecc[i], K = Kstar_kms,
omega = w_deg[i]*np.pi/180, gamma = 133.47)
ax0.plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'c', alpha=0.1, rasterized=True)
colors = ['#e41a1c', '#feb24c', '#e7298a', '#66a61e', '#a6761d', '#377eb8', '#756bb1']
alpha_ = 1
for i, name in enumerate(['FEROS', 'MagE', 'PEPSI', 'TRES', 'MIKE', 'APF']):
m = rv_names == name
ax0.errorbar(Time(jds[m], format='jd').decimalyear, rvs[m], yerr= rv_errs[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax1.errorbar(Time(jds[m], format='jd').decimalyear, rvs[m], yerr= rv_errs[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax0.plot([], [], 'c', label = r'$Gaia\,\,{\rm predictions}$')
ax0.set_xlabel(r'$\rm year$', fontsize=20)
ax0.set_ylabel(r'$\rm RV\,[km\,s^{-1}]$', fontsize=20)
ax0.legend(loc = 'upper left', frameon = False, ncol =2)
ax0.set_xlim(np.min(Time(t_grid, format= 'jd').decimalyear), np.max(Time(t_grid, format= 'jd').decimalyear))
ax1.plot([], [], 'k', label = r'$\rm joint\,\,fit$')
ax1.legend(loc = 'upper left', frameon = False, ncol =2)
best = np.argmax(sampler_withrvs.flatlnprobability)
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1 = sampler_withrvs.flatchain.T
randints = np.random.randint(0, len(period), 50)
t_grid = np.linspace(np.min(jds) - 600, np.max(jds) + 600, 3000)
for i in randints:
Kstar_kms = get_Kstar_kms(period = period[i], inc_deg = inc_deg[i], m1 = m1[i], m2 = m2[i], ecc = ecc[i])
rv_pred = get_radial_velocities(t = t_grid-jd_ref, P = period[i], T_p = t_peri[i], ecc = ecc[i], K = Kstar_kms,
omega = w_deg[i]*np.pi/180, gamma = v_com[i])
ax1.plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'k', alpha=0.1, rasterized=True)
Kstar_kms = get_Kstar_kms(period = period[best], inc_deg = inc_deg[best], m1 = m1[best], m2 = m2[best], ecc = ecc[best])
rv_pred = get_radial_velocities(t = t_grid-jd_ref, P = period[best], T_p = t_peri[best], ecc = ecc[best], K = Kstar_kms,
omega = w_deg[best]*np.pi/180, gamma = v_com[best])
resid = rvs - np.interp(jds, t_grid, rv_pred)
for i, name in enumerate(['FEROS', 'MagE', 'PEPSI', 'TRES', 'MIKE', 'APF']):
m = rv_names == name
ax2.errorbar(Time(jds[m], format='jd').decimalyear, resid[m], yerr= rv_errs[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax2.set_xlim(2021.5, 2024.5)
ax1.set_xlim(2021.5, 2024.5)
ax0.set_xlim(2021.5, 2024.5)
ax2.set_ylim(-0.8, 0.8)
ax2.plot([2021.5, 2024.5], [0, 0], 'k--', lw=0.5)
ax0.set_xticklabels([])
ax1.set_xticklabels([])
ax0.set_ylim(103, 162)
ax1.set_ylim(103, 162)
ax2.set_xlabel(r'$\rm year$', fontsize=20)
ax0.set_ylabel(r'$\rm RV\,\,[km\,s^{-1}]$', fontsize=20)
ax1.set_ylabel(r'$\rm RV\,\,[km\,s^{-1}]$', fontsize=20)
ax2.set_ylabel(r'$\rm residual\,\,[km\,s^{-1}]$', fontsize=20)
Out[10]:
Text(0, 0.5, '$\\rm residual\\,\\,[km\\,s^{-1}]$')
In [11]:
# note that the Gaia-only RV curve here is more noisy than what is shown in Fig 2 of the paper, because here we
# inflated the Gaia uncertainties by a factor of 3. In the paper we did that too for the mass constraints reported
# in the tables.
In [12]:
print('done!')
done!
In [ ]: