In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os
from isochrones.dartmouth import Dartmouth_Isochrone # Get using pip install isochrones
from isochrones import StarModel
from astropy import units as u, constants
G = constants.G
%matplotlib inline
In [2]:
# Read in MCMC samples. You can make this with the Do_MCMC_Fit notebook
mcmc_fname = 'data/SB2_samples.npy'
samples = np.load(mcmc_fname)
In [3]:
# Use Dartmouth isochrones to get the primary star mass from its parameters (from Simbad and Brugameyer et al (in prep) )
dar = Dartmouth_Isochrone()
Teff = (6546, 42)
logg = (3.9, 0.11)
feh = (-0.1, 0.05)
J = (3.803, 0.264)
H = (3.648, 0.248)
K = (3.502, 0.362)
model = StarModel(dar, Teff=Teff, logg=logg, feh=feh, J=J, H=H, K=K)
model.fit() # This will take a while, especially if using emcee instead of MultiNest
Now that we have everything read in and calculated, we can start deriving new quantities. We will first measure the mass ratio (q) from the semiamplitudes K1 and K2. For all variables, the values I quote come from percentiles of the samples. For large numbers, that is the same as doing
$v = \int_{-\infty }^{y} P(x)dx$
where y = 0.5 for the central value, and y=0.16 and 0.84 give the $1\sigma$ confidence intervals.
Once we have those measured, we need to get all the samples the same size so that we can do algebra to them and derive things like the companion mass, companion temperature, inclination, and semi-major axis.
In [4]:
# Get percentiles of the stuff in samples
K1_samples = samples[:, 0]
K2_samples = samples[:, 1]
q_samples = K1_samples / K2_samples
In [5]:
# Make numpy arrays with all the parameters I care about. Make sure they are the same length!
idx = np.random.randint(0, samples.shape[0], model.samples.shape[0])
q = q_samples[idx]
M1 = model.samples['mass'].values * u.M_sun
P = samples[idx, 2] * u.day
K1 = K1_samples[idx] * u.km/u.s
In [6]:
# Calculate the companion mass
M2 = M1*q
l, m, h = np.percentile(M2, [16, 50, 84])
print ('M2 = {:.3f} +{:.3f} / -{:.3f} Msun'.format(m, h-m, m-l))
In [7]:
# Use the dartmouth isochrones to get samples of the companion temperature
T2 = dar.Teff(M2, model.samples.age, model.samples.feh)
l, m, h = np.percentile(T2, [16, 50, 84])
print ('T2 = {:.0f} +{:.0f} / -{:.0f} Kelvin'.format(m, h-m, m-l))
In [8]:
# Get the inclination
sin3i = P/(2*np.pi*G) * K1**3 / (q*M1)
sini = (sin3i**(1./3.)).decompose()
i = np.arcsin(sini).to(u.degree)
l, m, h = np.percentile(i, [16, 50, 84])
print ('i = {:.0f} +{:.0f} / -{:.0f} degrees'.format(m, h-m, m-l))
In [9]:
# Get the semi-major axis
a3 = G*M1*(1+q)*P**2 / (4*np.pi**2)
a = (a3**(1./3.)).to(u.AU)
l, m, h = np.percentile(a, [16, 50, 84])
print ('a = {:.1f} +{:.1f} / -{:.1f} AU'.format(m, h-m, m-l))
In [ ]: