Skip to content

Instantly share code, notes, and snippets.

@megbedell

megbedell/barycentric.ipynb

Created January 14, 2019 04:55
Show Gist options
  • Save megbedell/5f3746ebe4332e82bf46a65fde929751 to your computer and use it in GitHub Desktop.
Save megbedell/5f3746ebe4332e82bf46a65fde929751 to your computer and use it in GitHub Desktop.
Download ZIP
barycentric
Display the source blob
Display the rendered blob
Raw
barycentric.ipynb
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Raw
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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment