Commit 3840fc17 authored by Monica Rainer's avatar Monica Rainer
Browse files

Fix BERV computations

parent 4dfd170e
...@@ -4,8 +4,8 @@ import datetime ...@@ -4,8 +4,8 @@ import datetime
import numpy import numpy
from numpy import sin, cos, tan, sqrt, arcsin from numpy import sin, cos, tan, sqrt, arcsin
#from PyAstronomy.pyaC import pyaErrors as PE #from PyAstronomy.pyaC import pyaErrors as PE
import pyaErrors as PE #import pyaErrors as PE
import six.moves as smo #import six.moves as smo
def daycnv(xjd, mode="idl"): def daycnv(xjd, mode="idl"):
...@@ -76,10 +76,10 @@ def daycnv(xjd, mode="idl"): ...@@ -76,10 +76,10 @@ def daycnv(xjd, mode="idl"):
Converted to IDL V5.0 W. Landsman September 1997 Converted to IDL V5.0 W. Landsman September 1997
""" """
if not mode in ('idl', 'dtlist', 'dt'): #if not mode in ('idl', 'dtlist', 'dt'):
raise(PE.PyAValError("Unknown mode: " + str(mode), \ # raise(PE.PyAValError("Unknown mode: " + str(mode), \
where="daycnv", \ # where="daycnv", \
solution="Use any of 'idl', 'dtlist', or 'dt'.")) # solution="Use any of 'idl', 'dtlist', or 'dt'."))
# Adjustment needed because Julian day starts at noon, calendar day at midnight # Adjustment needed because Julian day starts at noon, calendar day at midnight
...@@ -115,15 +115,15 @@ def daycnv(xjd, mode="idl"): ...@@ -115,15 +115,15 @@ def daycnv(xjd, mode="idl"):
if mode == 'dtlist': if mode == 'dtlist':
if not iterable: if not iterable:
return [yr[0], mn[0], day[0], hour[0], minute[0], sec[0], msec[0]] return [yr[0], mn[0], day[0], hour[0], minute[0], sec[0], msec[0]]
return [[yr[i], mn[i], day[i], hour[i], minute[i], sec[i], msec[i]] for i in smo.range(len(yr))] return [[yr[i], mn[i], day[i], hour[i], minute[i], sec[i], msec[i]] for i in xrange(len(yr))]
# Return datetime object # Return datetime object
dts = [datetime.datetime(*(yr[i], mn[i], day[i], hour[i], minute[i], sec[i], msec[i])) for i in smo.range(len(yr))] dts = [datetime.datetime(*(yr[i], mn[i], day[i], hour[i], minute[i], sec[i], msec[i])) for i in xrange(len(yr))]
if not iterable: if not iterable:
return dts[0] return dts[0]
return dts return dts
if not iterable: if not iterable:
return [yr[0], mn[0], day[0], hr[0]] return [yr[0], mn[0], day[0], hr[0]]
return [[yr[i], mn[i], day[i], hr[i]] for i in smo.range(len(yr))] return [[yr[i], mn[i], day[i], hr[i]] for i in xrange(len(yr))]
...@@ -284,7 +284,7 @@ def bprecess(ra, dec, mu_radec = None, \ ...@@ -284,7 +284,7 @@ def bprecess(ra, dec, mu_radec = None, \
s1 = r1/rmag ; s1_dot = r1_dot/rmag s1 = r1/rmag ; s1_dot = r1_dot/rmag
s = s1 s = s1
for j in smo.range(2): for j in xrange(2):
r = s1 + A - ((s * A).sum())*s r = s1 + A - ((s * A).sum())*s
s = r/rmag s = r/rmag
...@@ -901,10 +901,10 @@ def helio_jd(date, ra, dec, B1950 = False, TIME_DIFF = False): ...@@ -901,10 +901,10 @@ def helio_jd(date, ra, dec, B1950 = False, TIME_DIFF = False):
# Because XYZ uses default B1950 coordinates, we'll convert everything to B1950 # Because XYZ uses default B1950 coordinates, we'll convert everything to B1950
if date > 2.4e6: #if date > 2.4e6:
PE.warn(PE.PyAValError("The given Julian Date ( " + str(date) + ") is exceedingly large far a reduced JD.", # PE.warn(PE.PyAValError("The given Julian Date ( " + str(date) + ") is exceedingly large far a reduced JD.",
solution="Did you forget to subtract 2.4e6?", # solution="Did you forget to subtract 2.4e6?",
where="helio_jd")) # where="helio_jd"))
if not B1950: if not B1950:
bpresult = bprecess(ra,dec) bpresult = bprecess(ra,dec)
......
from __future__ import print_function, division from __future__ import print_function, division
import numpy as np import numpy as np
from astroTimeLegacy import premat, daycnv, precess, helio_jd from astroTimeLegacy import premat, daycnv, precess, helio_jd
from idlMod import idlMod
#from PyAstronomy.pyaC import pyaErrors as PE def idlMod(a, b):
import pyaErrors as PE """
import six Emulate 'modulo' behavior of IDL.
import six.moves as smo
Parameters
----------
a : float or array
Numerator
b : float
Denominator
Returns
-------
IDL modulo : float or array
The result of IDL modulo operation.
"""
if isinstance(a, np.ndarray):
s = np.sign(a)
m = np.mod(a, b)
m[(s < 0)] -= b
else:
m = a % b
if a < 0: m -= b
return m
def baryvel(dje, deq): def baryvel(dje, deq):
""" """
...@@ -234,7 +255,7 @@ def baryvel(dje, deq): ...@@ -234,7 +255,7 @@ def baryvel(dje, deq):
pertld = 0.0 pertld = 0.0
pertr = 0.0 pertr = 0.0
pertrd = 0.0 pertrd = 0.0
for k in smo.range(15): for k in xrange(15):
a = idlMod((dcargs[k,0] + dt*dcargs[k,1]), dc2pi) a = idlMod((dcargs[k,0] + dt*dcargs[k,1]), dc2pi)
cosa = np.cos(a) cosa = np.cos(a)
sina = np.sin(a) sina = np.sin(a)
...@@ -269,7 +290,7 @@ def baryvel(dje, deq): ...@@ -269,7 +290,7 @@ def baryvel(dje, deq):
pertld = 0.0 pertld = 0.0
pertp = 0.0 pertp = 0.0
pertpd = 0.0 pertpd = 0.0
for k in smo.range(3): for k in xrange(3):
a = idlMod((dcargm[k,0] + dt*dcargm[k,1]), dc2pi) a = idlMod((dcargm[k,0] + dt*dcargm[k,1]), dc2pi)
sina = np.sin(a) sina = np.sin(a)
cosa = np.cos(a) cosa = np.cos(a)
...@@ -293,7 +314,7 @@ def baryvel(dje, deq): ...@@ -293,7 +314,7 @@ def baryvel(dje, deq):
dxbd = dxhd*dc1mme dxbd = dxhd*dc1mme
dybd = dyhd*dc1mme dybd = dyhd*dc1mme
dzbd = dzhd*dc1mme dzbd = dzhd*dc1mme
for k in smo.range(4): for k in xrange(4):
plon = forbel[k+3] plon = forbel[k+3]
pomg = sorbel[k+1] pomg = sorbel[k+1]
pecc = sorbel[k+9] pecc = sorbel[k+9]
...@@ -453,8 +474,8 @@ def helcorr(obs_long, obs_lat, obs_alt, ra2000, dec2000, jd, equinox=2000.0, pma ...@@ -453,8 +474,8 @@ def helcorr(obs_long, obs_lat, obs_alt, ra2000, dec2000, jd, equinox=2000.0, pma
# East longitudes are positive # East longitudes are positive
obs_long = -obs_long obs_long = -obs_long
if jd < 2.4e6: # if jd < 2.4e6:
PE.warn(PE.PyAValError("The given Julian Date (" + str(jd) + ") is exceedingly small. Did you subtract 2.4e6?")) # PE.warn(PE.PyAValError("The given Julian Date (" + str(jd) + ") is exceedingly small. Did you subtract 2.4e6?"))
# Covert JD to Gregorian calendar date # Covert JD to Gregorian calendar date
xjd = jd xjd = jd
......
...@@ -875,15 +875,18 @@ def berv_corr(hdr): ...@@ -875,15 +875,18 @@ def berv_corr(hdr):
hjd = hdr[CONFIG['DRS_MJD'][0]] hjd = hdr[CONFIG['DRS_MJD'][0]]
return barycorr, hjd return barycorr, hjd
# target coordinates
radec = coord.SkyCoord(ra,dec, unit=(u.hourangle, u.deg)) radec = coord.SkyCoord(ra,dec, unit=(u.hourangle, u.deg))
ra = radec.ra.value ra = radec.ra.value
dec = radec.dec.value dec = radec.dec.value
equinox = float(hdr[CONFIG['KEYS']['EQUINOX']]) equinox = float(hdr[CONFIG['KEYS']['EQUINOX']])
# proper motions
try: try:
pma = float(hdr[CONFIG['KEYS']['PMA']]) pma = float(hdr[CONFIG['KEYS']['PMA']])
pmd = float(hdr[CONFIG['KEYS']['PMD']]) pmd = float(hdr[CONFIG['KEYS']['PMD']])
except: except:
print CONFIG['KEYS']['PMA']
pma = 0.0 pma = 0.0
pmd = 0.0 pmd = 0.0
...@@ -892,10 +895,21 @@ def berv_corr(hdr): ...@@ -892,10 +895,21 @@ def berv_corr(hdr):
if abs(pmd) > 100: if abs(pmd) > 100:
pmd = 0.0 pmd = 0.0
# JD (if MJD then convert to JD) + half exposure time
try: try:
mjd = hdr[CONFIG['DRS_MJD'][0]] mjd = float(hdr[CONFIG['DRS_MJD'][0]])
except: except:
mjd = hdr[CONFIG['KEYS']['MJD']] mjd = float(hdr[CONFIG['KEYS']['MJD']])
if mjd < 100000:
mjd = mjd + 2400000.5
try:
expt = float(hdr[CONFIG['KEYS']['EXPTIME']])/2.0
except:
expt = 0.0
mjd = mjd + (expt/(86400.0))
barycorr, hjd, vbar, vdiurnal = baryvel.helcorr(longitude,latitude,elevation,ra,dec,mjd,equinox,pma,pmd) barycorr, hjd, vbar, vdiurnal = baryvel.helcorr(longitude,latitude,elevation,ra,dec,mjd,equinox,pma,pmd)
return barycorr, hjd return barycorr, hjd
......
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