In [ ]:

    

import os
from os import path
from astropy.io import fits
import astropy.units as u
from astropy.table import Table
from astropy.time import Time

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import h5py

from twobody import KeplerOrbit
from thejoker import JokerParams, TheJoker, JokerSamples, RVData
from twoface.sample_prior import make_prior_cache
from twoface.io import load_samples
from twoface.plot import plot_data_orbits


    

In [ ]:

    

t0 = Time('2013-01-01')
baseline = 5 * u.yr # similar to APOGEE2
pars = JokerParams(P_min=1*u.day, P_max=1024*u.day)

K = 1 * u.km/u.s
err = 150 * u.m/u.s


    

In [ ]:

    

N = 3

# t = np.random.uniform(t0.mjd, (t0 + baseline).mjd, size=N)

# t = 10 ** np.random.uniform(np.log10(t0.mjd), np.log10((t0 + baseline).mjd), size=N)
t = t0 + 10 ** np.random.uniform(0, np.log10(baseline.to(u.day).value), size=N) * u.day
t = Time(t, format='mjd')
t.max() - t.min()


    

In [ ]:

    

orb = KeplerOrbit(P=150.*u.day, e=0., omega=0*u.deg, M0=0*u.deg)
rv = K * orb.unscaled_radial_velocity(t)
data = RVData(t, rv, stddev=np.ones_like(rv.value) * err)
_ = data.plot()


    

In [ ]:

    

joker = TheJoker(pars)


    

In [ ]:

    

samples = joker.iterative_rejection_sample(data, n_requested_samples=100, 
                                           n_prior_samples=1000000)


    

In [ ]:

    

_ = plot_data_orbits(data, samples)


    

In [ ]:

    

logP = np.log(samples['P'].to(u.day).value)
plt.scatter(logP, samples['e'], alpha=0.25, marker='.')
plt.axvline(np.log(orb.P.to(u.day).value))


    

In [ ]:

    

from scipy.stats import scoreatpercentile


    

In [ ]:

    

pers = scoreatpercentile(logP, [5, 10, 50, 90, 95])
pers[1]-pers[0], pers[4]-pers[3]


    

In [ ]:

    

pers - np.log(orb.P.to(u.day).value)


    

In [ ]:

    




    

In [ ]: