#!/usr/bin/env python # coding: utf-8 # In[1]: import numpy as np import matplotlib.pyplot as plt import emcee import corner get_ipython().run_line_magic('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) # 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) # In[9]: get_ipython().run_line_magic('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) # 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) # 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!') # In[ ]: