'''
.. module:: timing
.. moduleauthor:: Anna-Lea Lesage
timing functions used for MASCARA
Includes:
- Timing routines:
- check_discontinuity
- astrotiming
- calendar2jd
- jd2calendar
- jd2epoch
'''
__version__ = '15.10.02'
### ----- TIMING ROUTINES ----- ###
# To verify that the time cadence of the image is regular...
import numpy as np
from datetime import datetime, timedelta
[docs]def checkdiscontinuity(obstime, texp):
''' Given an exposure time -texp- and an array containing at least the last 2
dates of observation, checktimescale verifies that the latest exposure is
made in the continuity of the previous one. This is relevant in cases of
a serie of observations during one night. If the 2 exposures diverge by
more than twice the exposure time, than a discontinuity is found, and the
function returns True.
Input:
- obstime, an array containing the dates of observation.
- texp, a scalar for the exposure time, and including the read-
out time of the detector.
Output:
- a boolean value, True if a discontinuity is found, else False
'''
try:
dim = len(obstime)
except AttributeError:
dim = 1
if dim < 2:
return False
else:
to = jd2calendar(obstime.pop(-2))
ti = jd2calendar(obstime[-1])
tmp = ti-to
if tmp.seconds < 2*texp:
return False
else:
return True
###### ------------------------------#####
### ------ TIME CONVERSIONS -----------###
[docs]def calendar2jd(*args, **kwargs):
''' Calendar to Jd converts a datetime object into a Julian date object
Input:
- date, a datetime object or a timedelta object or an array of datetime
objects
Output:
- JD, the Julian date in np.float64
'''
if len(args) == 1:
date = args[0]
# Case 1 it is a datetime object:
if type(date) == datetime:
utcoffset = date.utcoffset()
### No datetime awareness -> maybe in the keyword UTC?
if utcoffset is None:
UTC = kwargs.get("UTC", True)
if not UTC:
print "Not enough information for the conversion to JD"
return
utcoffset = timedelta(seconds=0)
date = date - utcoffset
yr = date.year
mth = date.month
day = date.day
hr = np.float64(date.hour)
mn = np.float64(date.minute)
sec = np.float64(date.second + date.microsecond * 1e-6)
# Case 2: it's an array of datetime objects:
elif type(date) == np.ndarray:
getyears = np.vectorize(lambda x:x.year)
getmonth = np.vectorize(lambda x:x.month)
getday = np.vectorize(lambda x:x.day)
gethour = np.vectorize(lambda x:np.float64(x.hour))
getmin = np.vectorize(lambda x:np.float64(x.minute))
getsec = np.vectorize(lambda x:np.float64(x.second + x.microsecond*1e-6))
yr = getyears(date)
mth = getmonth(date)
day = getday(date)
hr = gethour(date)
mn = getmin(date)
sec = getsec(date)
# Case 3: it is a list/tuple with year, month, day, hour, minutes, seconds format
elif isinstance(date, timedelta):
from mascara.constants import JDsec
sec = date.total_seconds()
return sec*JDsec
elif len(args) == 6:
yr, mth, day, hr, mn, sec = args
hr = np.float64(hr)
mn = np.float64(mn)
sec = np.float64(sec)
# Condition for the leap years:
m3 = mth < 3
L = np.zeros(np.size(mth), dtype=int)
if np.size(L) == 1:
if m3:
L = -1
else:
L[m3] = -1
#if m3:
# L = -1
#else: L = 0
thour = hr + mn/60 + sec/3600
jd = np.float64(day - 32075 + 1461*(yr + 4800 +L)/4 + 367*(mth -2 -L*12)/12 - \
3 *(yr + 4900 + L)/400) + thour/24 - 0.5
return jd
[docs]def jd2calendar(jd, tz = None):
''' Convert a Julian day date (JD) into a calendar date, using a datetime format.
Input:
- JD, a scalar or a vector, in float64 precision
Output:
- date, a datetime object in UTC timezone
'''
try:
import pytz
utc = pytz.utc
except ImportError:
print 'No timezone information available'
if type(jd) != np.float64 and type(jd) != np.ndarray:
raise SyntaxError('JD is a float64 scalar...')
frac = jd%1 + 0.5
xjd = np.array(jd - jd%1, dtype=long)
if np.isscalar(frac) and frac>1:
xjd = xjd + 1
frac = frac -1
elif not np.isscalar(frac):
xjd[frac>1] = xjd[frac>1]+1
frac[frac>1] = frac[frac>1]-1
hr = frac*24
l = xjd + 68569
n = (4*l / 146097)
l = l - (146097*n + 3) / 4
yr = 4000*(l+1) / 1461001
l = l - 1461*yr / 4 + 31
mth = 80*l / 2447
day = l - 2447*mth / 80
l = mth/11
mth = mth + 2 - 12*l
yr = 100*(n-49) + yr + l
mn = (hr%1)*60
sec = (mn%1)*60
usec = np.array((sec%1)*1e6, dtype=long)
hr = np.array(hr, dtype = long)
mn = np.array(mn, dtype=long)
sec = np.array(sec, dtype=long)
if (hr == 24.0).any():
hr[hr==24] = hr[hr==24] - 24.0
day[hr==24] = day[hr==24] + 1
def todate(year, month, day, hour, mn, sec, usec):
return datetime(year, month, day, hour, mn, sec, usec, tzinfo=utc)
indate = np.vectorize(todate)
date = indate(yr, mth, day, hr, mn, sec, usec)
if len(date.shape) > 1:
date.shape = date.shape[1]
return date
else:
return date[()]
[docs]def jd2epoch(jd, Julian=True): # TESTED WORKING!!!
''' Conversion of a given Julian date to an epoch expressed in decimal
years.
Input: jd, the Julian date to convert, if jd is None, it takes the
current date
Keyword: Julian. If not set, the reference epoch is 1900
Output: epoch. The date converted in decimal years.
'''
from mascara.constants import JD2000, JDyr
jd = np.array(jd, copy=False)
if Julian:
epoch = 2000.00 + (jd - JD2000)/(JDyr/100)
else:
epoch = 1900.00 + (jd-2415020.31352)/365.242198781
return epoch
[docs]def epoch2jd(epoch):
'''Convert an epoch into a Julian Date'''
print 'do something'
return
