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_BH2 = np.array([-0.28826183, 0.7105017 , -0.52574706, -0.21933624, 0.08650991,
-0.14516878, 0.6302555 , -0.3771263 , 0.5500957 , -0.73675364,
-0.16729414, 0.30398202, -0.20453878, -0.8524959 , 0.44096807,
0.93901193, -0.483403 , 0.7650907 , -0.40073904, 0.760147 ,
-0.01839051, -0.89376473, 0.5308219 , -0.76396966, 0.40606955,
-0.8563013 , 0.0148804 , -0.9522201 , -0.09896384, 0.21685918,
-0.17396761, -0.91208214, 0.5061753 , 0.94459516, 0.0334551 ,
-0.06152862, -0.25910425, 0.15509658, -0.2196958 , -0.46669292,
0.19075292, 0.5711617 , -0.15632574, 0.12374099, 0.5730349 ,
-0.6190153 , 0.3727271 , -0.5269465 , 0.511369 , -0.69945323,
-0.20429476, -0.72100574, 0.716509 , -0.25523147, -0.3031498 ,
0.18343681, -0.07696132, 0.14444804, -0.7800489 , 0.65127426,
0.67564595, 0.33760983, -0.36696544, 0.71416605, 0.66939396,
-0.4554749 , 0.76010406, -0.08568522, 0.51160055, -0.71235156,
0.8925269 , 0.42780364, 0.7911277 , -0.8147842 , 0.5150208 ,
0.14617956, -0.67924684, 0.6091381 , 0.5922516 , -0.23008141,
0.45837632, -0.8792923 , 0.9372642 , 0.6043152 , 0.73485035,
-0.7568019 , 0.6737804 , 0.27924234, -0.70115626, 0.7368394 ,
0.93496126, 0.8697892 , -0.5762039 , 0.7589055 , 0.08591607,
0.41520476, -0.5027001 , 0.85735095, -0.8110985 , -0.46167427,
-0.4300237 , -0.5127689 , -0.06291884, 0.42949262, 0.30797017])
# ra, dec, parallax, pmra, pmdec, A_pred, B_pred, F_pred, G_pred, C_pred, H_pred, v_com, ecc, period, t_peri
mu_vec_BH2 = np.array([ 2.07569716e+02, -5.92390050e+01, 8.59191858e-01, -1.04778102e+01,
-4.61109332e+00, 2.47823048e+00, 2.45044157e+00, -1.98274465e+00,
3.00783355e+00, 2.06147425e+00, -1.95257220e+00, -3.20952270e+00,
5.32327853e-01, 1.35228894e+03, 5.27881716e+01])
err_vec_BH2 = np.array([0.06482769, 0.027924482, 0.018293705, 0.09809765, 0.062035598, 0.046954542,
0.096332304, 0.10347517, 0.10905123, 0.112943575, 0.12393499, 0.6448445,
0.015242959, 45.469105, 3.8733842])
# the RV data
jds_BH2 = np.array([2459814.4989, 2459916.8408, 2459917.8428, 2459918.8523, 2459919.8432, 2459920.8374, 2459921.8421,
2459923.8343, 2459924.8366, 2459932.8561, 2459935.8419, 2459940.8282, 2459947.8362, 2459957.8404,
2459957.8637, 2459958.8580, 2459959.8584, 2459960.8611, 2459961.8795, 2459962.8667, 2459965.8774,
2459966.8686, 2459967.8619, 2459972.7330, 2459984.8437, 2459985.8316, 2459986.8430, 2459987.8336,
2459989.7996, 2459990.7861, 2459992.8014, 2459994.7108, 2459997.8192, 2459999.8278, 2460000.8586,
2460001.8520, 2460002.8348, 2460003.8399, 2460004.8448, 2460005.8399, 2460006.84823, 2460012.7997,
2460013.8823, 2460014.8624, 2460015.8741, 2460017.8436, 2460018.6749, 2460018.8868])
rvs_BH2 = np.array([10.99, -2.05, -2.24, -2.60, -2.84, -3.06, -3.40, -4.08, -4.21, -6.69, -7.67, -9.43, -12.02,
-15.95, -15.89, -16.29, -16.70, -17.16, -17.52, -17.95, -19.26, -19.58, -19.99, -22.10, -26.82,
-27.18, -27.56, -27.93, -28.60, -28.96, -29.66, -30.23, -31.28, -31.88, -32.15, -32.50, -32.76,
-33.05, -33.31, -33.55, -33.82, -35.14, -35.38, -35.55, -35.71, -36.06, -36.22, -36.27])
rv_errs_BH2 = np.array([0.04, 0.04, 0.04, 0.03, 0.04, 0.08, 0.03, 0.09, 0.03, 0.03, 0.10, 0.03, 0.04, 0.1, 0.05,
0.05, 0.03, 0.04, 0.04, 0.04, 0.03, 0.03, 0.03, 0.1, 0.03, 0.03, 0.04, 0.03, 0.03, 0.03,
0.03, 0.1, 0.04, 0.03, 0.04, 0.04, 0.04, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03,
0.03, 0.1, 0.04])
rv_names_BH2 = np.array(['FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS',
'UVES', 'FEROS', 'FEROS', 'UVES', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS',
'FEROS', 'FEROS', 'FEROS', 'UVES', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS',
'FEROS', 'UVES', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS',
'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'FEROS', 'UVES', 'FEROS'])
# prior on the mass of the luminous star
m1_obs, m1_err = 1.07, 0.19
# the reference epoch of the gaia t_periastron; i.e., JD 2016.0
jd_ref = 2457389.0
In [3]:
# convert the correlation matrix to a covariance matrix
Nparam = 15
triangle_nums = (np.arange(1, Nparam)*(np.arange(1, Nparam) - 1)//2)
cov_mat_BH2 = []
count = 0
for i, sigx in enumerate(err_vec_BH2):
this_row = []
for j, sigy in enumerate(err_vec_BH2):
if i == j:
this_row.append(err_vec_BH2[i]**2)
elif j < i:
this_row.append(cov_mat_BH2[j][i])
else:
this_idx = triangle_nums[j-1] + i
this_row.append(corr_vec_BH2[this_idx]*sigx*sigy)
cov_mat_BH2.append(this_row)
cov_mat_BH2 = np.array(cov_mat_BH2)
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_and_a1_au(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)
a1_au = a_au*(q/(1+q))
return a0_mas, a1_au
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 * np.sin(inc_deg*np.pi/180)**3 / \
(period*86400 * (1 - ecc**2)**(3/2)))**(1/3)
return Kstar/1000 # km/s
def lnL_with_covariance(y_vec, mu_vec, Sigma_mat):
'''
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
'''
return -0.5*np.dot((y_vec - mu_vec).T, np.dot(np.linalg.inv(Sigma_mat), 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, q = inc_deg*np.pi/180, omega_deg*np.pi/180, w_deg*np.pi/180, m2/m1
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_BH2-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_BH2)**2/rv_errs_BH2**2)
else:
lnL = 0
a0_mas, a1_au = get_a0_mas_and_a1_au(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) )
C_pred = a1_au*np.sin(w)*np.sin(inc)
H_pred = a1_au*np.cos(w)*np.sin(inc)
y_vec = np.array([ra, dec, parallax, pmra, pmdec, A_pred, B_pred, F_pred, G_pred, C_pred, H_pred, v_com, ecc, period, t_peri])
lnL += lnL_with_covariance(y_vec = y_vec, mu_vec = mu_vec_BH2, Sigma_mat = cov_mat_BH2)
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([207.569716, -59.239005, 0.859191858, -10.4778102,
-4.61109332, 1352.28894, 0.532327853, 35.1, 268.7, 130.2, 52.8, -3.2, 9.0, 1.07])
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)
# 3000 more steps
sampler_norvs.reset()
state = sampler_norvs.run_mcmc(state, 3000, progress = True)
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [00:11<00:00, 425.19it/s] 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3000/3000 [00:07<00:00, 424.97it/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, 3000, progress = True)
4%|█████▎ | 177/5000 [00:06<02:43, 29.46it/s]/var/folders/_4/hh3dgzmd13l0wvjyf6tkz1cm0000gn/T/ipykernel_68425/3181202268.py:26: RuntimeWarning: invalid value encountered in double_scalars Kstar = (2*np.pi*G*(m2*Msun) * (m2/(m1 + m2))**2 * np.sin(inc_deg*np.pi/180)**3 / \ /var/folders/_4/hh3dgzmd13l0wvjyf6tkz1cm0000gn/T/ipykernel_68425/3181202268.py:8: RuntimeWarning: invalid value encountered in double_scalars a_au = (((period*86400)**2 * G * (m1*Msun + m2*Msun)/(4*np.pi**2))**(1/3.))/AU 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 5000/5000 [02:50<00:00, 29.24it/s] 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 3000/3000 [01:42<00:00, 29.15it/s]
In [9]:
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)
# fix the formatting
fig.axes[10].set_title('$\\rm ecc$ = ${%.3f}_{-%.3f}^{+%.3f}$' % (np.median(half_sampler.T[1]), np.median(half_sampler.T[1])-np.percentile(half_sampler.T[1], 16), np.percentile(half_sampler.T[1], 86) - np.median(half_sampler.T[1])))
#fig.axes[20].set_title(' '+fig.axes[20].get_title())
#fig.axes[30].set_title(' '+fig.axes[30].get_title())
#fig.axes[40].set_title(' '+fig.axes[40].get_title())
#fig.axes[50].set_title(' '+fig.axes[50].get_title())
#fig.axes[60].set_title(' '+fig.axes[60].get_title())
#fig.axes[70].set_title(' '+fig.axes[70].get_title())
#fig.axes[80].set_title(' '+fig.axes[80].get_title())
Out[9]:
Text(0.5, 1.0, '$\\rm ecc$ = ${0.517}_{-0.001}^{+0.001}$')
In [10]:
# posterior predictive check for the "with RVs" solution.
from astropy.time import Time
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_BH2) - 600, np.max(jds_BH2) + 600, 3000)
colors = ['#e41a1c', '#feb24c', '#e7298a', '#66a61e', '#a6761d', '#377eb8', '#756bb1']
f = plt.figure(figsize = (12, 11))
ax0 = f.add_axes([0.05, 0.7, 0.9, 0.25])
ax1 = f.add_axes([0.05, 0.375, 0.9, 0.25])
ax2 = f.add_axes([0.05, 0.05, 0.4, 0.25])
ax3 = f.add_axes([0.55, 0.05, 0.4, 0.25])
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])
ax0.plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'c', alpha=0.1, rasterized=True)
alpha_ = 1
for i, name in enumerate(['FEROS', 'UVES']):
m = rv_names_BH2 == name
ax0.errorbar(Time(jds_BH2[m], format='jd').decimalyear, rvs_BH2[m], yerr= rv_errs_BH2[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax1.errorbar(Time(jds_BH2[m], format='jd').decimalyear, rvs_BH2[m], yerr= rv_errs_BH2[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax2.errorbar(Time(jds_BH2[m], format='jd').decimalyear, rvs_BH2[m], yerr= rv_errs_BH2[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax3.errorbar(Time(jds_BH2[m], format='jd').decimalyear, rvs_BH2[m], yerr= rv_errs_BH2[m], fmt='o', color = colors[i], label = name, alpha=alpha_)
ax0.plot([], [], 'c', label = r'$Gaia\,\,{\rm solution\,\,predictions}$')
ax0.set_xlabel(r'$\rm year$', fontsize=20)
ax0.set_ylabel(r'$\rm RV\,[km\,s^{-1}]$', fontsize=20)
ax0.legend(loc = 'lower left', frameon = False)
ax0.set_xlim(np.min(Time(t_grid, format= 'jd').decimalyear), np.max(Time(t_grid, format= 'jd').decimalyear))
ra, dec, parallax, pmra, pmdec, period, ecc, inc_deg, omega_deg, w_deg, t_peri, v_com, m2, m1 = sampler_withrvs.flatchain.T
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)
ax2.plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'k', alpha=0.1, rasterized=True)
ax3.plot(Time(t_grid, format= 'jd').decimalyear, rv_pred, 'k', alpha=0.1, rasterized=True)
ax2.set_xlim(2022.9, 2022.96)
ax2.set_ylim(-5, -1)
ax3.set_xlim(2023.1, 2023.25)
ax3.set_ylim(-40, -25)
ax1.plot([], [], 'k', label = r'${\rm Joint}\,\,Gaia\,\,+ {\rm RV\,\,fit}$')
ax1.set_xlabel(r'$\rm year$', fontsize=20)
ax1.set_ylabel(r'$\rm RV\,[km\,s^{-1}]$', fontsize=20)
ax1.legend(loc = 'lower left', frameon = False)
ax1.set_xlim(np.min(Time(t_grid, format= 'jd').decimalyear), np.max(Time(t_grid, format= 'jd').decimalyear))
ax2.set_xlabel(r'$\rm year - 2022$', fontsize=20)
ax2.set_ylabel(r'$\rm RV\,[km\,s^{-1}]$', fontsize=20)
ax3.set_xlabel(r'$\rm year - 2023$', fontsize=20)
ax3.set_ylabel(r'$\rm RV\,[km\,s^{-1}]$', fontsize=20)
import matplotlib.patches as patches
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
def custom_mark_inset(axA, axB, fc='None', ec='k'):
xx = axB.get_xlim()
yy = axB.get_ylim()
xy = (xx[0], yy[0])
width = xx[1] - xx[0]
height = yy[1] - yy[0]
pp = axA.add_patch(patches.Rectangle(xy, width, height, fc=fc, ec=ec))
p1 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[0], yy[0]), xyB=(xx[0], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB, linewidth=1, color='0.2'))
p2 = axA.add_patch(patches.ConnectionPatch(
xyA=(xx[1], yy[0]), xyB=(xx[1], yy[1]),
coordsA='data', coordsB='data',
axesA=axA, axesB=axB,linewidth=1, color='0.2'))
return pp, p1, p2
pp, p1, p2 = custom_mark_inset(ax1, ax2)
ppb, p1b, p2b = custom_mark_inset(ax1, ax3)
from matplotlib import ticker
def format_fn(value, tick_number):
str_ = f"{value:.2f}"
return str_.replace( str_.split('.')[0], '0')
ax2.xaxis.set_major_formatter(ticker.FuncFormatter(format_fn))
ax3.xaxis.set_major_formatter(ticker.FuncFormatter(format_fn))
In [ ]: