barycentric

barycentric.ipynb

barycentric.py

# coding: utf-8

# In[1]:


import wobble
import numpy as np
import matplotlib.pyplot as plt


# #### Get the Gaia coordinates:

# In[2]:


from astropy.time import Time
from astroquery.gaia import Gaia
import astropy.units as u
from astropy import coordinates
coordinates.solar_system_ephemeris.set('jpl')
from astropy.coordinates import SkyCoord, EarthLocation


# In[3]:


coord = SkyCoord(ra=269.4486, dec=4.7379807, unit=(u.degree, u.degree), frame='icrs')
width = u.Quantity(0.01, u.degree)
height = u.Quantity(0.01, u.degree)
r = Gaia.query_object_async(coordinate=coord, width=width, height=height);


# In[4]:


star = r[r['source_id'] == 4472832130942575872]


# In[5]:


sc = SkyCoord(ra=star['ra'][0]*u.deg, dec=star['dec'][0]*u.deg)


# #### Try a naive BERV calculation neglecting proper motion effects:

# In[6]:


loc = EarthLocation.of_site('lasilla')


# In[7]:


data = wobble.Data('barnards_e2ds.hdf5', filepath='/Users/mbedell/python/wobble/data/', orders=[60])


# In[8]:


dates = Time(data.dates, format='jd')
pipeline_bervs = data.bervs * u.m / u.s


# In[9]:


naive_bervs = sc.radial_velocity_correction(obstime=dates, location=loc)


# In[10]:


plt.plot_date(dates.plot_date, naive_bervs - pipeline_bervs, 'r.')
plt.ylabel('Resids w.r.t \n HARPS pipeline BERVs (m/s)', fontsize=12);


# #### Try it using coordinates calculated with Gaia proper motions:

# In[11]:


sc = SkyCoord(ra=star['ra'][0]*u.deg, dec=star['dec'][0]*u.deg,
             pm_ra_cosdec=star['pmra'][0]*u.mas/u.year,
             pm_dec=star['pmdec'][0]*u.mas/u.year,
             obstime=Time(2015.5, format='decimalyear'))
#bervs = sc.radial_velocity_correction(obstime=dates, location=loc)


# In[12]:


def calc_new_coords(sc, date):
    new_sc = SkyCoord(ra=sc.ra + sc.pm_ra_cosdec / np.cos(sc.dec.radian) * (date - sc.obstime),
                      dec=sc.dec + sc.pm_dec * (date - sc.obstime),
                      obstime=date)
    return new_sc


# In[13]:


sc.ra, sc.dec


# In[14]:


new_sc = sc.apply_space_motion(new_obstime=Time(2005., format='decimalyear'))
new_sc.ra, new_sc.dec


# In[15]:


new_sc = calc_new_coords(sc, Time(2005., format='decimalyear'))
new_sc.ra, new_sc.dec


# In[16]:


scs = [calc_new_coords(sc, Time(d, format='jd')) for d in dates]


# In[17]:


bervs = [s.radial_velocity_correction(location=loc).value for s in scs]


# In[18]:


plt.plot_date(dates.plot_date, bervs - data.bervs, 'r.')
plt.ylabel('Resids w.r.t \n HARPS pipeline BERVs (m/s)', fontsize=12);


# In[19]:


plt.plot(data.bervs, bervs - data.bervs, 'r.')


# In[20]:


plt.plot_date(dates.plot_date, bervs - naive_bervs.value, 'r.')


# In[21]:


plt.plot(bervs, bervs - naive_bervs.value, 'r.')


# #### Try to replicate pipeline with HARPS header coords:

# In[22]:


from astropy.io import fits
from tqdm import tqdm
scs = []
harps_bervs = []
for f in tqdm(data.filelist):
    sp = fits.open(f)
    sc = SkyCoord(ra=sp[0].header['RA'] * u.deg,
                  dec=sp[0].header['DEC'] * u.deg,
                  pm_ra_cosdec = sp[0].header['HIERARCH ESO TEL TARG PMA'] * u.arcsec / u.year,
                  pm_dec = sp[0].header['HIERARCH ESO TEL TARG PMD'] * u.arcsec / u.year,
                  obstime=Time('J2000'))
    scs.append(sc)
    date = Time(sp[0].header['HIERARCH ESO DRS BJD'], format='jd')
    new_sc = calc_new_coords(sc, date)
    harps_bervs.append(new_sc.radial_velocity_correction(location=loc).value)


# In[43]:


plt.hist([sc.pm_dec.value for sc in scs])
plt.axvline(star['pmdec'][0], c='k')
plt.xlabel('PM Dec (arcsec/yr)', fontsize=16);


# In[45]:


plt.hist([sc.pm_ra_cosdec.value for sc in scs])
plt.axvline(star['pmra'][0], c='k')
plt.xlabel('PM RA (arcsec/yr)', fontsize=16);


# In[25]:


plt.plot_date(dates.plot_date, harps_bervs - data.bervs, 'r.')
plt.ylabel('Resids w.r.t \n HARPS pipeline BERVs (m/s)', fontsize=12);


# In[26]:


plt.plot_date(dates.plot_date, np.array(harps_bervs) - np.array(bervs), 'r.')


# In[27]:


plt.plot(data.bervs, harps_bervs - data.bervs, 'r.')


# In[28]:


adc_ra = sp[0].header['HIERARCH ESO INS ADC1 RA']
adc_dec = sp[0].header['HIERARCH ESO INS ADC1 DEC']


# In[29]:


sp[0].header['MJD-OBS'] # MJD at observation start


# In[30]:


sp[0].header['MJD-END']


# In[31]:


sp[0].header['HIERARCH ESO DRS BJD'] - 2400000.5


# In[32]:


(sp[0].header['HIERARCH ESO DRS BJD'] - 2400000.5 - sp[0].header['MJD-OBS'])*24.*60.


# In[33]:


sp[0].header['EXPTIME']/60.


# ### how much does flux weighting matter?

# In[34]:


sc = SkyCoord(ra=star['ra'][0]*u.deg, dec=star['dec'][0]*u.deg,
             pm_ra_cosdec=star['pmra'][0]*u.mas/u.year,
             pm_dec=star['pmdec'][0]*u.mas/u.year,
             obstime=Time(2015.5, format='decimalyear'))


# In[ ]:


dates2 = dates + np.random.normal(0., 1./24./60., len(dates)) # 1-minute randomness in midpoint


# In[ ]:


scs = [calc_new_coords(sc, Time(d, format='jd')) for d in dates2]


# In[ ]:


bervs2 = [s.radial_velocity_correction(location=loc).value for s in scs]


# In[ ]:


plt.scatter(bervs, np.array(bervs2) - np.array(bervs))


# In[ ]:


np.std(np.array(bervs2) - np.array(bervs))


# In[ ]:


mid = 56966.24165454
start = 56966.23644620
end = 56966.24686289


# In[ ]:


(end - start)/2. + start # naively calculated midpoint time


# In[ ]:


2456966.74612107 - 2400000.5 # BJD


# In[ ]:


end - 22.6/60./60./24.


# In[36]:


sc


# In[39]:


with coordinates.solar_system_ephemeris.set('jpl'):
    print(sc.radial_velocity_correction(location=loc))


# In[40]:


with coordinates.solar_system_ephemeris.set('builtin'):
    print(sc.radial_velocity_correction(location=loc))