#!/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_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) # 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) # 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()) # 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[ ]: