<jupyter_start><jupyter_text>This is an example to infer eccentricities for individual wide binaries, reproducing the result in Fig. 8 in [Hwang et al. 2022](https://arxiv.org/abs/2111.01789)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from astropy.table import Table
import astropy.units as u
%matplotlib inline 
from matplotlib import rcParams
rcParams.update({'font.size': '18'})
rcParams.update({'font.family': 'sans-serif'})
rcParams.update({'axes.facecolor': 'white'})
rcParams.update({'axes.edgecolor': 'black'})
rcParams.update({'axes.labelcolor': 'black'})
rcParams.update({'xtick.top': 'True'})
rcParams.update({'xtick.major.pad': '6.0'})
rcParams.update({'xtick.color': 'k'})
rcParams.update({'xtick.direction': 'in'})
rcParams.update({'ytick.right': 'True'})
rcParams.update({'ytick.major.pad': '6.0'})
rcParams.update({'ytick.color': 'k'})
rcParams.update({'ytick.direction': 'in'})
rcParams.update({'legend.fontsize': '18'})
rcParams.update({'figure.figsize': '8.0, 6.0'})
rcParams.update({'figure.dpi': '100'})
rcParams.update({'figure.facecolor': 'white'})
rcParams.update({'figure.edgecolor': 'white'})
rcParams.update({'image.cmap': 'rainbow'})<jupyter_output><empty_output><jupyter_text>Read wide binary table, which contains two rows from the eDR3 wide binary catalog ([El-Badry et al. 2021](https://ui.adsabs.harvard.edu/abs/2021MNRAS.506.2269E/abstract))<jupyter_code>wb_table = Table.read('example_binary.fits')<jupyter_output><empty_output><jupyter_text>Take a look at the content:<jupyter_code>wb_table<jupyter_output><empty_output><jupyter_text>Calculate their v-r angles:<jupyter_code>#difference in ra-dec coordinates
wb_table['dra'] = (wb_table['ra2'] - wb_table['ra1']) * np.cos((wb_table['dec1']+wb_table['dec2'])/2. * np.pi / 180.)
wb_table['ddec'] = (wb_table['dec2'] - wb_table['dec1'])
wb_table['dradec'] = np.sqrt(wb_table['dra']**2 + wb_table['ddec']**2)

#difference in proper motions
wb_table['dpmra'] = wb_table['pmra2'] - wb_table['pmra1']
wb_table['dpmdec'] = wb_table['pmdec2'] - wb_table['pmdec1']
wb_table['dpmradec'] = np.sqrt(wb_table['dpmra']**2 + wb_table['dpmdec']**2)

#v-r angles
wb_table['vr_angle'] = 180. / np.pi * np.arccos((wb_table['dra'] * wb_table['dpmra'] + wb_table['ddec'] * wb_table['dpmdec']) / wb_table['dpmradec'] / wb_table['dradec'])


#uncertainties of proper motion differences and v-r angles
wb_table['dpm_error'] = 1. / wb_table['dpmradec'] * np.sqrt(
        wb_table['dpmra']**2 * (wb_table['pmra_error1']**2 + wb_table['pmra_error2']**2) +
        wb_table['dpmdec']**2 * (wb_table['pmdec_error1']**2 + wb_table['pmdec_error2']**2) 
    )

wb_table['dpm_over_error'] = wb_table['dpmradec'] / wb_table['dpm_error']

wb_table['vr_angle_error'] = 180. / np.pi / wb_table['dpm_over_error']

wb_table['vr_angle_error'].unit = u.deg<jupyter_output><empty_output><jupyter_text>These two wide binaries have very different v-r angles: one is 3.6 deg and one is 86.5 deg.<jupyter_code>wb_table['vr_angle', 'vr_angle_error']<jupyter_output><empty_output><jupyter_text>Below we infer the probability distribution function (PDF) of their individual eccentricties, $p(e|\gamma_{obs})$ where $\gamma_{obs}$ is the observed v-r angle. The Bayes' theorem says$$p(e|\gamma_{obs}) \propto p(\gamma_{obs}|e) p(e)$$And we marginalize over the uncertainty of the v-r angle:$$p(e|\gamma_{obs}) \propto \int p(\gamma_{obs}|\gamma_{true}) p(\gamma_{true}|e) p(e) d\gamma_{true},$$where $p(\gamma_{obs}|\gamma_{true})$ is a Gaussian (truncated at 0 and 180 deg) with a width of v-r angle uncertainty. The term $p(\gamma_{true}|e)$ is pre-computed from simulations. Here we load the result. See more details in "Demo 1 - Using p(gamma|e) grid" notebook.<jupyter_code>from scipy import interpolate

e_list, vr_angle_list_center, hist_list = np.load('./grid_data/data.npy', allow_pickle=True)
p_gamma_bar_e = interpolate.interp2d(e_list, vr_angle_list_center, hist_list.T, kind='cubic')<jupyter_output><empty_output><jupyter_text>First we define a function to compute the truncated Gaussian for $ p(\gamma_{obs}|\gamma_{true}) $. scipy 1.7 is needed here.<jupyter_code>import scipy
from scipy import stats

def vr_angle_pdf_gaussian(vr_angle, vr_angle_error, vr_angle_list=np.arange(0., 181., 1.)):
    
    gaussian = scipy.stats.norm(vr_angle, vr_angle_error).pdf(vr_angle_list) 
    return gaussian / np.sum(gaussian)
scipy.__version__<jupyter_output><empty_output><jupyter_text>Then we compute the integration $\int p(\gamma_{obs}|\gamma_{true}) p(\gamma_{true}|e) d\gamma_{true}$ for these two wide binaries. We do this integration for each $e$ from 0 to 1.<jupyter_code>p_gamma_bar_e_list = []

#grid for integration
e_dense_list = np.arange(0., 1., 0.01)
vr_angle_list_to_integrate = np.arange(0., 181., 1)

N_vr_angle_list = len(vr_angle_list_to_integrate)


for pair in wb_table:
    
    #initialize the array as zero
    p_e_bar_gamma = np.zeros_like(e_dense_list)
    
    vr_angle_obs = pair['vr_angle']
    vr_angle_error_obs = pair['vr_angle_error']
    
    #compute the Gaussian
    vr_angle_pdf = vr_angle_pdf_gaussian(vr_angle_obs, vr_angle_error_obs)

    #integrate (i.e. summation) over vr angle from 0 to 180 
    p_e_bar_gamma = np.sum(p_gamma_bar_e(e_dense_list,vr_angle_list_to_integrate)* vr_angle_pdf[:,None], axis=0)

    p_gamma_bar_e_list.append(p_e_bar_gamma)<jupyter_output><empty_output><jupyter_text>We then need to determine the prior $p(e)$. In [Hwang et al. 2022](https://arxiv.org/abs/2111.01789), we measure the eccentricity distribution (parametrized by a power law, $p(e)\propto e^\alpha$) as a function of binary separations. Then we describe the relation between $\alpha$ (eccentricity distribution) and binary separation using a tanh() function. Below we calculate this $\alpha$ for these two binaries based on their separations.<jupyter_code>A = 0.91
B = 2.29
C = 0.61
D = 0.44

wb_table['alpha'] = A * np.tanh( (np.log10(wb_table['sep_AU']) - B) / C ) + D<jupyter_output><empty_output><jupyter_text>Check their alpha:<jupyter_code>wb_table['sep_AU', 'alpha']<jupyter_output><empty_output><jupyter_text>For the purpse of this jupyter notebook (and Fig. 8 in [Hwang et al. 2022](https://arxiv.org/abs/2111.01789) ), we explore how different priors affect the eccentricity inference using $\alpha$=0, 0.5, and 1. Below, for each $\alpha$, we calculate the $p(e|\gamma_{obs})$ by multiplying the $\int p(\gamma_{obs}|\gamma_{true}) p(\gamma_{true}|e) d\gamma_{true}$ computed above with the prior $p(e)=(1+\alpha)e^\alpha$.<jupyter_code>e_dense_list = np.arange(0., 1., 0.01)

for i, pair in enumerate(wb_table):

    alpha= 0.5
    
    #load the pre-computed integration
    p_gamma_bar_e = p_gamma_bar_e_list[i]
    lw=3
    
    #alpha=0.5 prior
    plt.plot(
        e_dense_list,
        p_gamma_bar_e * (1. + alpha) * e_dense_list**alpha,
        c='C%d' %(i), lw=lw
    )
    
    #alpha=0 prior
    plt.plot(
        e_dense_list,
        p_gamma_bar_e,
        c='C%d' %(i), ls=':', lw=lw
    )
    
    #alpha=1 prior
    plt.plot(
        e_dense_list,
        p_gamma_bar_e * (2. * e_dense_list),
        c='C%d' %(i), ls='--', lw=lw
    )
    
    
    
    plt.text(
        [0.8, 0.1][i],
        [0.016, 0.016][i],
        r'$\gamma=%.1f^{\circ}$' %(pair['vr_angle'])
    )

    i+=1

    
    
plt.plot(
    [],[], ls=':',
    label=r'$\alpha=0$ prior',
    c='k', lw=lw
)

plt.plot(
    [],[], 
    label=r'$\alpha=0.5$ prior',
    c='k', lw=lw
)


plt.plot(
    [],[], ls='--',
    label=r'$\alpha=1$ prior',
    c='k', lw=lw
)

plt.legend()

plt.ylim(0, 0.028)

plt.xlabel('Eccentricity')
plt.ylabel(r'Posterior $p(e|\gamma_{obs})$')
plt.yticks([0, 0.01, 0.02])

plt.show()<jupyter_output>Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
Substituting with a symbol from Computer Modern.
{'fraction_non_alphanumeric': 0.16064603624707188, 'fraction_numerical': 0.02527431882628529, 'mean_word_length': 4.960323291697281, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 4, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '649084', 'license': 'MIT', 'n_tokens_mistral': 3275, 'n_tokens_neox': 3071, 'n_words': 655, 'path': 'Demo 2 - Infer eccentricities for individual wide binaries.ipynb', 'repo_name': 'HC-Hwang/Eccentricity-of-wide-binaries', 'size': '8111'}