Commit ee116ecc authored by Monica Rainer's avatar Monica Rainer
Browse files

Change barycentric correction

parent a90d1650
......@@ -16,12 +16,15 @@
"""
from drslib.config import CONFIG
from drslib.berv import baryvel
#from drslib.berv import baryvel
from astropy import constants as const
from astropy import units as u
from astropy import coordinates as coord
from astropy.io import fits
from astropy.time import Time
import barycorrpy as bc
import numpy as np
import numpy.polynomial.polynomial as poly
......@@ -1095,7 +1098,7 @@ def check_keywords(header,filename):
#--------------------- Compute barycentric correction -------------------
def berv_corr(hdr):
def berv_corr_old(hdr):
# CORREZIONE BARICENTRICA + BJD
# TNG coordinates
latitude = 28.754
......@@ -1129,24 +1132,96 @@ def berv_corr(hdr):
pmd = 0.0
# JD (if MJD then convert to JD) + half exposure time
try:
expt = float(hdr[CONFIG['KEYS']['EXPTIME']])/2.0
except:
expt = 0.0
try:
mjd = float(hdr[CONFIG['DRS_MJD'][0]])
except:
mjd = float(hdr[CONFIG['KEYS']['MJD']])
mjd = mjd + (expt/(86400.0))
if mjd < 100000:
mjd = mjd + 2400000.5
barycorr, hjd, bjd, vbar, vdiurnal = baryvel.helcorr(longitude,latitude,elevation,ra,dec,mjd,equinox,pma,pmd)
return barycorr, hjd, bjd
#--------------------- Compute barycentric correction NEW -------------------
def berv_corr(hdr):
# CORREZIONE BARICENTRICA + BJD
# TNG coordinates
latitude = 28.754
longitude = -17.889056
elevation = 2387.2
tng = coord.EarthLocation.from_geodetic(latitude,-longitude,elevation)
try:
ra = hdr[CONFIG['KEYS']['RA']]
dec = hdr[CONFIG['KEYS']['DEC']]
except:
barycorr = 0.0
hjd = hdr[CONFIG['DRS_MJD'][0]]
return barycorr, hjd
# proper motions
try:
pma = float(hdr[CONFIG['KEYS']['PMA']])
pmd = float(hdr[CONFIG['KEYS']['PMD']])
except:
pma = 0.0
pmd = 0.0
if abs(pma) > 100:
pma = 0.0
if abs(pmd) > 100:
pmd = 0.0
# target coordinates
radec = coord.SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
ra = radec.ra.value
dec = radec.dec.value
equinox = float(hdr[CONFIG['KEYS']['EQUINOX']])
t = Time(equinox,format='jyear')
epoch = t.jd
# JD (if MJD then convert to JD) + half exposure time
try:
expt = float(hdr[CONFIG['KEYS']['EXPTIME']])/2.0
except:
expt = 0.0
mjd = mjd + (expt/(86400.0))
try:
mjd = float(hdr[CONFIG['DRS_MJD'][0]])
except:
mjd = float(hdr[CONFIG['KEYS']['MJD']])
mjd = mjd + (expt/(86400.0))
if mjd < 100000:
mjd = mjd + 2400000.5
times = Time(mjd, format='jd', scale='utc', location=tng)
barycorr, hjd, bjd, vbar, vdiurnal = baryvel.helcorr(longitude,latitude,elevation,ra,dec,mjd,equinox,pma,pmd)
#ltt_bary = times.light_travel_time(radec)
#bjd = times.tdb + ltt_bary
#bjd = bjd.value
ltt_helio = times.light_travel_time(radec, 'heliocentric')
hjd = times.utc + ltt_helio
hjd = hjd.value
return barycorr, hjd, bjd
#berv = bc.get_BC_vel(JDUTC=mjd, ra=ra, dec=dec, lat = latitude, longi = longitude, alt = elevation, pmra = pma, pmdec = pmd, epoch=epoch, ephemeris = 'de430', leap_update=True)
berv = bc.get_BC_vel(JDUTC=mjd, ra=ra, dec=dec, lat = latitude, longi = longitude, alt = elevation, pmra = pma, pmdec = pmd, epoch=epoch, leap_update=True)
berv = berv[0][0]
corr_time = bc.utc_tdb.JDUTC_to_BJDTDB(mjd, ra=ra, dec=dec, lat = latitude, longi = longitude, alt = elevation, pmra = pma, pmdec = pmd, epoch=epoch, leap_update=True)
corr_time = corr_time[0][0]
return berv, hjd, corr_time
#--------------------- Create s1d output -------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment