In [1]:
import emcee
from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [163]:
more_info = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/mdm-targets.txt",
names=True, dtype=None)
more_info
Out[163]:
In [2]:
_phase,_vtemplate = np.loadtxt("/Users/adrian/projects/triand-rrlyrae/data/RVtemplates/Halpha.txt").T
template = interp1d(_phase, _vtemplate) # define an interpolator to get template value at arbitrary phases
In [3]:
# observed parameters
a = 35.6 # ± 2.5
b = 78.2 # ± 2.4
In [4]:
def model(phase, vsys, a, b, AV):
Arv_Ha = a*AV + b
return Arv_Ha*template(phase) + vsys
def ln_likelihood(p, phase, v, sigma, a, a_err, b, b_err, AV):
(vsys,) = p
var = (sigma**2 + template(phase)**2 * (AV**2*a_err**2 + b_err**2))
lnB = -0.5 * var / (sigma*a_err*b_err)
lnL = -0.5 * (v - model(phase, vsys, a, b, AV))**2 / var
return lnB + lnL
def ln_posterior(p, phase, v, sigma, a, a_err, b, b_err, AV):
return ln_likelihood(p, phase, v, sigma, a, a_err, b, b_err, AV).sum()
Test for RRL14 -- from Ally's fit:
$$v_\gamma = -228. \pm 25.$$rv phase:
In [5]:
phase = np.array([0.19, 0.26])
rv = np.array([-256., -210.])
rv_err = np.array([15., 18.])
AV = 1.21*0.453
lls = []
vsyss = np.linspace(-250,-150,100)
for vsys in vsyss:
lls.append(ln_likelihood([vsys], phase, rv, rv_err, a, 2.5, b, 2.4, AV).sum())
plt.plot(vsyss, lls)
Out[5]:
In [30]:
sampler = emcee.EnsembleSampler(nwalkers=32, dim=1, lnpostfn=ln_posterior,
args=(phase, rv, rv_err, a, 2.5, b, 2.4, AV))
In [31]:
p0 = np.random.normal(-200, 10., size=(sampler.k,1))
In [54]:
pos,prob,state = sampler.run_mcmc(p0, 100)
sampler.reset()
pos,prob,state = sampler.run_mcmc(pos, 256)
In [64]:
vhel = model(0.27, sampler.flatchain[:,0], a, b, AV)
v = np.median(vhel)
verr = np.sqrt((a*AV + b + 3.4)**2 * (0.06**2 + (0.1*1.47)**2))
v, verr
Out[64]:
In [15]:
ally_d = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/TriAnd_RRL_26mar15.csv",
skiprows=0, dtype=None, names=True, delimiter=',')
In [6]:
d = np.genfromtxt("/Users/adrian/projects/triand-rrlyrae/data/TriAnd_RRL_RV_intermediate.csv",
delimiter=',', names=True, dtype=None)
d.dtype.names
Out[6]:
In [26]:
def fit_v(phase, rv, rv_err, A_R):
sampler = emcee.EnsembleSampler(nwalkers=32, dim=1, lnpostfn=ln_posterior,
args=(phase, rv, rv_err, a, 2.5, b, 2.4, AV))
p0 = np.random.normal(-200, 10., size=(sampler.k,1))
pos,prob,state = sampler.run_mcmc(p0, 100)
sampler.reset()
pos,prob,state = sampler.run_mcmc(pos, 256)
extra = np.sqrt(np.sum((-(phase-0.5)**2/3 + 0.06)**2))
vsys = sampler.flatchain[:,0]
vhel = model(0.27, vsys, a, b, AV)
v = np.median(vhel)
verr = np.sqrt((a*AV + b + 3.4)**2 * (extra**2 + (0.1*1.47)**2))
return vsys, v, verr
In [61]:
ally_d
Out[61]:
In [33]:
ally_dd.dtype.names
Out[33]:
In [190]:
rows = []
for name in np.unique(d['name']):
dd = d[d['name'] == name]
ally_dd = ally_d[ally_d['name'] == name]
info = more_info[more_info['name'] == name]
print("---------------\n{0}:".format(name))
AV = 1.21 * dd['A_R'][0]
# data to fit
phase = (dd['Phase'] + 0.5/info['period']) % 1.
print("Phase shift: {0}".format((phase - dd['Phase'])[0]))
# phase = dd['phase']
ix = (phase > 0.05) & (phase < 0.85)
if ix.sum() == 0:
continue
phase = phase[ix]
RV_hel = dd['RV_hel'][ix]
RV_err = dd['Err'][ix]
print("\t Phases: {0}".format(phase))
print("{0} measurements initially, {1} with good phase".format(len(dd), len(phase)))
vsys,v,verr = fit_v(phase, RV_hel, RV_err, AV)
plt.figure()
plt.title(name)
plt.errorbar(phase, RV_hel, RV_err, marker='o', ecolor='#666666', linestyle='none')
phase = np.linspace(0.,0.94,100)
for i in range(128):
plt.plot(phase, model(phase, vsys[np.random.randint(len(vsys))], a, b, AV),
marker=None, color='k', alpha=0.05)
rows.append((name, v, verr, ally_dd['Vsys'], ally_dd['Err'], ally_dd['ra'], ally_dd['dec'], ally_dd['dist']))
apw = np.array(rows, dtype=[('name','S12'),('v_apw','f8'),('verr_apw','f8'),('v_ally','f8'),('verr_ally','f8'),
('ra','f8'), ('dec','f8'), ('dist','f8')])
In [191]:
apw.dtype.names
Out[191]:
In [192]:
plt.figure(figsize=(6,6))
plt.plot(apw['v_apw'], apw['v_ally'], linestyle='none')
plt.plot(np.linspace(apw['v_apw'].min(),apw['v_apw'].max(),100),
np.linspace(apw['v_apw'].min(),apw['v_apw'].max(),100), marker=None)
plt.xlabel('Adrian')
plt.ylabel('Ally')
plt.title("Best-fit velocity")
Out[192]:
In [193]:
plt.figure(figsize=(6,6))
plt.plot(apw['verr_apw'], apw['verr_ally'], linestyle='none')
plt.plot(np.linspace(13,27,100), np.linspace(13,27,100), marker=None)
plt.xlim(13,27)
plt.ylim(13,27)
plt.xlabel('Adrian')
plt.ylabel('Ally')
plt.title("Velocity uncertainty")
Out[193]:
In [194]:
from astropy.table import Table
In [195]:
tbl = Table(apw)
tbl.rename_column('v_apw', 'Vsys')
tbl.rename_column('verr_apw', 'Err')
In [196]:
tbl.write("/Users/adrian/projects/triand-rrlyrae/data/apw_velocities.csv", format='ascii', delimiter=',')
In [197]:
!cat /Users/adrian/projects/triand-rrlyrae/data/apw_velocities.csv
In [202]:
pub_tbl = Table(apw)
# fix names
pub_tbl['name'] = np.array([name.replace('l','L') for name in pub_tbl['name']])
pub_tbl['sort_name'] = np.array([int(name[9:]) for name in pub_tbl['name']])
pub_tbl.rename_column('v_apw', 'vgsr')
pub_tbl.rename_column('verr_apw', 'verr')
pub_tbl.remove_column('v_ally')
pub_tbl.remove_column('verr_ally')
pub_tbl.sort('sort_name')
pub_tbl.remove_column('sort_name')
pub_tbl
Out[202]:
In [204]:
pub_tbl.write("/Users/adrian/projects/triand-rrlyrae/data/publication_data.csv", format='ascii', delimiter=',')
In [187]:
# _cache = dict()
In [189]:
fig,axes = plt.subplots(2, 2, figsize=(10,8), sharex=True)
for i,name in enumerate(['TriAndRRl2','TriAndRRl26']):
dd = d[d['name'] == name]
info = more_info[more_info['name'] == name]
print("{0}: {1} measurements".format(name, len(dd)))
print(info['period'])
print(info['hjd0'])
AV = 1.21 * dd['A_R'][0]
phase = np.linspace(0.,0.94,100)
if name not in _cache:
vsys,v,verr = fit_v(dd['Phase'], dd['RV_hel'], dd['Err'], AV)
lines = np.array([model(phase, vsys[np.random.randint(len(vsys))], a, b, AV) for j in range(1000)])
_cache[name] = dict(vsys=vsys, v=v, verr=verr, lines=lines)
else:
vsys = _cache[name]['vsys']
v = _cache[name]['v']
verr = _cache[name]['verr']
lines = _cache[name]['lines']
# fit to RV curve
col = axes[:,i]
col[0].set_title(name.replace('l','L'))
quantiles = np.percentile(lines, [16, 84, 50], axis=0)
col[0].fill_between(phase, quantiles[0], quantiles[1], color='#dddddd')
col[0].plot(phase, quantiles[2], marker=None, lw=2., color='k')
col[0].errorbar(dd['Phase'], dd['RV_hel'], dd['Err'], marker='o', ecolor='#666666', linestyle='none')
# light curve
path = "/Users/adrian/projects/triand-rrlyrae/data/lc_{0}.dat".format(name.replace('l','L'))
lc = np.genfromtxt(path, names=True, dtype=None, skiprows=20)
lc = lc[(lc['filterID'] == 2) & (lc['sextractorFlags'] == 0) & (lc['magerr'] < 0.15)]
col[1].set_xlabel('Phase')
phase = ((lc['epoch'] - info['hjd0'] + 0.5) % info['period']) / info['period']
# if '10' in name:
# phase = ((lc['epoch'] - info['hjd0'] + 0.08) % info['period']) / info['period']
col[1].errorbar(phase, lc['mag'], lc['magerr'], linestyle='none',
ecolor='#666666', marker='o', alpha=0.75, ms=4)
axes[1,1].set_xlim(-0.02, 1.02)
axes[0,0].set_ylabel(r'$v_r$ [${\rm km~s}^{-1}$]')
axes[1,0].set_ylabel(r'$R$ [${\rm mag}$]')
# left column, top
axes[0,0].set_ylim(-300, -150)
axes[0,0].set_yticks([-160, -200, -240, -280])
# left column, bottom
axes[1,0].set_ylim(17.65, 16.75)
axes[1,0].set_yticks([17.5, 17.3, 17.1, 16.9])
# right column, top
axes[0,1].set_ylim(-275, -90)
# axes[0,1].set_yticks([-200, -240, -280, -320])
axes[1,1].set_ylim(18.4, 17.2)
axes[1,1].set_yticks([18.2, 17.9, 17.6, 17.3])
fig.tight_layout()
fig.subplots_adjust(hspace=0.)
fig.savefig("/Users/adrian/projects/triand-rrlyrae/plots/rv_lc.pdf")
In [ ]: