#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# timeutils.py - Waqas Bhatti (wbhatti@astro.princeton.edu) - Sept 2013
'''
Contains various useful tools for dealing with time in astronomical contexts.
'''
#############
## LOGGING ##
#############
import logging
from astrobase import log_sub, log_fmt, log_date_fmt
DEBUG = False
if DEBUG:
level = logging.DEBUG
else:
level = logging.INFO
LOGGER = logging.getLogger(__name__)
logging.basicConfig(
level=level,
style=log_sub,
format=log_fmt,
datefmt=log_date_fmt,
)
LOGDEBUG = LOGGER.debug
LOGINFO = LOGGER.info
LOGWARNING = LOGGER.warning
LOGERROR = LOGGER.error
LOGEXCEPTION = LOGGER.exception
#############
## IMPORTS ##
#############
import os.path
import os
import numpy as np
# we need the astropy.time.Time to convert from UTC to TDB
# this should also add any leap seconds on top of the relevant corrections
import astropy.time as astime
# we need the jplephem package from Brandon Rhodes to import and use the JPL
# ephemerides
from jplephem.spk import SPK
#################
## JPL KERNELS ##
#################
modpath = os.path.abspath(os.path.dirname(__file__))
planetdatafile = os.path.join(modpath,'data/de430.bsp')
# we'll try to load the SPK kernel. if that fails, we'll download it direct from
# JPL so the source distribution is kept small
try:
# load the JPL kernel
jplkernel = SPK.open(planetdatafile)
HAVEKERNEL = True
except Exception:
# this is so we don't download the JPL kernel when this module is imported
# as part of a docs build on RTD
import os
RTD_INVOCATION = os.environ.get('RTD_IGNORE_HTLS_FAIL')
if not RTD_INVOCATION:
# this function is used to check progress of the download
def on_download_chunk(transferred,blocksize,totalsize):
progress = transferred*blocksize/float(totalsize)*100.0
print('{progress:.1f}%'.format(progress=progress),end='\r')
# this is the URL to the SPK
spkurl = (
'https://naif.jpl.nasa.gov/'
'pub/naif/generic_kernels/spk/planets/de430.bsp'
)
LOGINFO('JPL kernel de430.bsp not found. Downloading from:\n\n%s\n' %
spkurl)
from urllib.request import urlretrieve
localf, headerr = urlretrieve(
spkurl,planetdatafile,reporthook=on_download_chunk
)
if os.path.exists(localf):
print('\nDone.')
jplkernel = SPK.open(planetdatafile)
HAVEKERNEL = True
else:
print('failed to download the JPL kernel!')
HAVEKERNEL = False
else:
HAVEKERNEL = False
######################
## USEFUL CONSTANTS ##
######################
# physical constants
CLIGHT_KPS = 299792.458
# various JDs
JD1800 = 2378495.0
JD2000 = 2451545.0
JD2000INT = 2451545
JD2050 = 2469807.5
# conversion factors
MAS_P_YR_TO_RAD_P_DAY = 1.3273475e-11
ARCSEC_TO_RADIANS = 4.84813681109536e-6
KM_P_AU = 1.49597870691e8
SEC_P_DAY = 86400.0
# this is the Earth-Moon mass ratio
# needed for calculating the position of the earth's center wrt SSB.
# the JPL ephemerides provide the position of the Earth-moon barycenter with
# respect to the solar-system barycenter.
EMRAT = 81.30056941599857
#######################
## UTILITY FUNCTIONS ##
#######################
[docs]def precess_coordinates(ra, dec,
epoch_one, epoch_two,
jd=None,
mu_ra=0.0,
mu_dec=0.0,
outscalar=False):
'''Precesses target coordinates `ra`, `dec` from `epoch_one` to `epoch_two`.
This takes into account the jd of the observations, as well as the proper
motion of the target mu_ra, mu_dec. Adapted from J. D. Hartman's
VARTOOLS/converttime.c [coordprecess].
Parameters
----------
ra,dec : float
The equatorial coordinates of the object at `epoch_one` to precess in
decimal degrees.
epoch_one : float
Origin epoch to precess from to target epoch. This is a float, like:
1985.0, 2000.0, etc.
epoch_two : float
Target epoch to precess from origin epoch. This is a float, like:
2000.0, 2018.0, etc.
jd : float
The full Julian date to use along with the propermotions in `mu_ra`, and
`mu_dec` to handle proper motion along with the coordinate frame
precession. If one of `jd`, `mu_ra`, or `mu_dec` is missing, the proper
motion will not be used to calculate the final precessed coordinates.
mu_ra,mu_dec : float
The proper motion in mas/yr in right ascension and declination. If these
are provided along with `jd`, the total proper motion of the object will
be taken into account to calculate the final precessed coordinates.
outscalar : bool
If True, converts the output coordinates from one-element np.arrays to
scalars.
Returns
-------
precessed_ra, precessed_dec : float
A tuple of precessed equatorial coordinates in decimal degrees at
`epoch_two` taking into account proper motion if `jd`, `mu_ra`, and
`mu_dec` are provided.
'''
raproc, decproc = np.radians(ra), np.radians(dec)
if ((mu_ra != 0.0) and (mu_dec != 0.0) and jd):
jd_epoch_one = JD2000 + (epoch_one - epoch_two)*365.25
raproc = (
raproc +
(jd - jd_epoch_one)*mu_ra*MAS_P_YR_TO_RAD_P_DAY/np.cos(decproc)
)
decproc = decproc + (jd - jd_epoch_one)*mu_dec*MAS_P_YR_TO_RAD_P_DAY
ca = np.cos(raproc)
cd = np.cos(decproc)
sa = np.sin(raproc)
sd = np.sin(decproc)
if epoch_one != epoch_two:
t1 = 1.0e-3 * (epoch_two - epoch_one)
t2 = 1.0e-3 * (epoch_one - 2000.0)
a = ( t1*ARCSEC_TO_RADIANS * (23062.181 + t2*(139.656 + 0.0139*t2) +
t1*(30.188 - 0.344*t2+17.998*t1)) )
b = t1*t1*ARCSEC_TO_RADIANS*(79.280 + 0.410*t2 + 0.205*t1) + a
c = (
ARCSEC_TO_RADIANS*t1*(20043.109 - t2*(85.33 + 0.217*t2) +
t1*(-42.665 - 0.217*t2 - 41.833*t2))
)
sina, sinb, sinc = np.sin(a), np.sin(b), np.sin(c)
cosa, cosb, cosc = np.cos(a), np.cos(b), np.cos(c)
precmatrix = np.matrix([[cosa*cosb*cosc - sina*sinb,
sina*cosb + cosa*sinb*cosc,
cosa*sinc],
[-cosa*sinb - sina*cosb*cosc,
cosa*cosb - sina*sinb*cosc,
-sina*sinc],
[-cosb*sinc,
-sinb*sinc,
cosc]])
precmatrix = precmatrix.transpose()
x = (np.matrix([cd*ca, cd*sa, sd])).transpose()
x2 = precmatrix * x
outra = np.arctan2(x2[1],x2[0])
outdec = np.arcsin(x2[2])
outradeg = np.rad2deg(outra)
outdecdeg = np.rad2deg(outdec)
if outradeg < 0.0:
outradeg = outradeg + 360.0
if outscalar:
return float(outradeg), float(outdecdeg)
else:
return outradeg, outdecdeg
else:
# if the epochs are the same and no proper motion, this will be the same
# as the input values. if the epochs are the same, but there IS proper
# motion (and a given JD), then these will be perturbed from the input
# values of ra, dec by the appropriate amount of motion
return np.degrees(raproc), np.degrees(decproc)
############################
## EPOCHS FOR EPHEMERIDES ##
############################
def _single_true(iterable):
'''This returns True if only one True-ish element exists in `iterable`.
Parameters
----------
iterable : iterable
Returns
-------
bool
True if only one True-ish element exists in `iterable`. False otherwise.
'''
# return True if exactly one true found
iterator = iter(iterable)
# consume from "i" until first true or it's exhausted
has_true = any(iterator)
# carry on consuming until another true value / exhausted
has_another_true = any(iterator)
return has_true and not has_another_true
[docs]def get_epochs_given_midtimes_and_period(
t_mid,
period,
err_t_mid=None,
t0_fixed=None,
t0_percentile=None,
verbose=False
):
'''This calculates the future epochs for a transit, given a period and a
starting epoch
The equation used is::
t_mid = period*epoch + t0
Default behavior if no kwargs are used is to define `t0` as the median
finite time of the passed `t_mid` array.
Only one of `err_t_mid` or `t0_fixed` should be passed.
Parameters
----------
t_mid : np.array
A np.array of transit mid-time measurements
period : float
The period used to calculate epochs, per the equation above. For typical
use cases, a period precise to ~1e-5 days is sufficient to get correct
epochs.
err_t_mid : None or np.array
If provided, contains the errors of the transit mid-time
measurements. The zero-point epoch is then set equal to the average of
the transit times, weighted as `1/err_t_mid^2` . This minimizes the
covariance between the transit epoch and the period (e.g., Gibson et
al. 2013). For standard O-C analysis this is the best method.
t0_fixed : None or float:
If provided, use this t0 as the starting epoch. (Overrides all others).
t0_percentile : None or float
If provided, use this percentile of `t_mid` to define `t0`.
Returns
-------
tuple
This is the of the form `(integer_epoch_array, t0)`.
`integer_epoch_array` is an array of integer epochs (float-type),
of length equal to the number of *finite* mid-times passed.
'''
kwargarr = np.array([isinstance(err_t_mid,np.ndarray),
t0_fixed,
t0_percentile])
if not _single_true(kwargarr) and not np.all(~kwargarr.astype(bool)):
raise AssertionError(
'can have at most one of err_t_mid, t0_fixed, t0_percentile')
t_mid = t_mid[np.isfinite(t_mid)]
N_midtimes = len(t_mid)
if t0_fixed:
t0 = t0_fixed
elif isinstance(err_t_mid,np.ndarray):
# get the weighted average. then round it to the nearest transit epoch.
t0_avg = np.average(t_mid, weights=1/err_t_mid**2)
t0_options = np.arange(min(t_mid), max(t_mid)+period, period)
t0 = t0_options[np.argmin(np.abs(t0_options - t0_avg))]
else:
if not t0_percentile:
# if there are an odd number of times, take the median time as
# epoch=0. elif there are an even number of times, take the lower
# of the two middle times as epoch=0.
if N_midtimes % 2 == 1:
t0 = np.median(t_mid)
else:
t0 = t_mid[int(N_midtimes/2)]
else:
t0 = np.sort(t_mid)[int(N_midtimes*t0_percentile/100)]
epoch = (t_mid - t0)/period
# do not convert numpy entries to actual ints, because np.nan is float type
int_epoch = np.round(epoch, 0)
if verbose:
LOGINFO('epochs before rounding')
LOGINFO('\n{:s}'.format(repr(epoch)))
LOGINFO('epochs after rounding')
LOGINFO('\n{:s}'.format(repr(int_epoch)))
return int_epoch, t0
###########################
## JULIAN DATE FUNCTIONS ##
###########################
[docs]def unixtime_to_jd(unix_time):
'''This converts UNIX time in seconds to a Julian date in UTC (JD_UTC).
Parameters
----------
unix_time : float
A UNIX time in decimal seconds since the 1970 UNIX epoch.
Returns
-------
jd : float
The Julian date corresponding to the provided UNIX time.
'''
# use astropy's time module
jdutc = astime.Time(unix_time, format='unix', scale='utc')
return jdutc.jd
[docs]def datetime_to_jd(dt):
'''This converts a Python datetime object (naive, time in UT) to JD_UTC.
Parameters
----------
dt : datetime
A naive Python `datetime` object (e.g. with no tz attribute) measured at
UTC.
Returns
-------
jd : float
The Julian date corresponding to the `datetime` object.
'''
jdutc = astime.Time(dt, format='datetime',scale='utc')
return jdutc.jd
[docs]def jd_to_datetime(jd, returniso=False):
'''This converts a UTC JD to a Python `datetime` object or ISO date string.
Parameters
----------
jd : float
The Julian date measured at UTC.
returniso : bool
If False, returns a naive Python `datetime` object corresponding to
`jd`. If True, returns the ISO format string corresponding to the date
and time at UTC from `jd`.
Returns
-------
datetime or str
Depending on the value of `returniso`.
'''
tt = astime.Time(jd, format='jd', scale='utc')
if returniso:
return tt.iso
else:
return tt.datetime
[docs]def jd_now():
'''
Gets the Julian date at the current time.
Returns
-------
float
The current Julian date in days.
'''
return astime.Time.now().jd
[docs]def jd_to_mjd(jd):
'''
Converts Julian Date to Modified Julian Date.
Parameters
----------
jd : float
The Julian date measured at UTC.
Returns
-------
mjd : float
`mjd = jd - 2400000.5`
'''
return jd - 2400000.5
[docs]def mjd_to_jd(mjd):
'''
Converts Modified Julian date to Julian Date.
Parameters
----------
mjd : float
The Modified Julian date measured at UTC.
Returns
-------
jd : float
`jd = mjd + 2400000.5`
'''
return mjd + 2400000.5
##################################
## CONVERTING JD TO HJD and BJD ##
##################################
[docs]def jd_corr(jd,
ra, dec,
obslon=None,
obslat=None,
obsalt=None,
jd_type='bjd'):
'''Returns BJD_TDB or HJD_TDB for input JD_UTC.
The equation used is::
BJD_TDB = JD_UTC + JD_to_TDB_corr + romer_delay
where:
- JD_to_TDB_corr is the difference between UTC and TDB JDs
- romer_delay is the delay caused by finite speed of light from Earth-Sun
This is based on the code at:
https://mail.scipy.org/pipermail/astropy/2014-April/003148.html
Note that this does not correct for:
1. precession of coordinates if the epoch is not 2000.0
2. precession of coordinates if the target has a proper motion
3. Shapiro delay
4. Einstein delay
Parameters
----------
jd : float or array-like
The Julian date(s) measured at UTC.
ra,dec : float
The equatorial coordinates of the object in decimal degrees.
obslon,obslat,obsalt : float or None
The longitude, latitude of the observatory in decimal degrees and
altitude of the observatory in meters. If these are not provided, the
corrected JD will be calculated with respect to the center of the Earth.
jd_type : {'bjd','hjd'}
Conversion type to perform, either to Baryocentric Julian Date ('bjd')
or to Heliocenter Julian Date ('hjd').
Returns
-------
float or np.array
The converted BJD or HJD.
'''
if not HAVEKERNEL:
LOGERROR('no JPL kernel available, can\'t continue!')
return
# Source unit-vector
## Assume coordinates in ICRS
## Set distance to unit (kilometers)
# convert the angles to degrees
rarad = np.radians(ra)
decrad = np.radians(dec)
cosra = np.cos(rarad)
sinra = np.sin(rarad)
cosdec = np.cos(decrad)
sindec = np.sin(decrad)
# this assumes that the target is very far away
src_unitvector = np.array([cosdec*cosra,cosdec*sinra,sindec])
# Convert epochs to astropy.time.Time
## Assume JD(UTC)
if (obslon is None) or (obslat is None) or (obsalt is None):
t = astime.Time(jd, scale='utc', format='jd')
else:
t = astime.Time(jd, scale='utc', format='jd',
location=('%.5fd' % obslon,
'%.5fd' % obslat,
obsalt))
# Get Earth-Moon barycenter position
## NB: jplephem uses Barycentric Dynamical Time, e.g. JD(TDB)
## and gives positions relative to solar system barycenter
barycenter_earthmoon = jplkernel[0,3].compute(t.tdb.jd)
# Get Moon position vectors from the center of Earth to the Moon
# this means we get the following vectors from the ephemerides
# Earth Barycenter (3) -> Moon (301)
# Earth Barycenter (3) -> Earth (399)
# so the final vector is [3,301] - [3,399]
# units are in km
moonvector = (jplkernel[3,301].compute(t.tdb.jd) -
jplkernel[3,399].compute(t.tdb.jd))
# Compute Earth position vectors (this is for the center of the earth with
# respect to the solar system barycenter)
# all these units are in km
pos_earth = (barycenter_earthmoon - moonvector * 1.0/(1.0+EMRAT))
if jd_type == 'bjd':
# Compute BJD correction
## Assume source vectors parallel at Earth and Solar System
## Barycenter
## i.e. source is at infinity
# the romer_delay correction is (r.dot.n)/c where:
# r is the vector from SSB to earth center
# n is the unit vector from
correction_seconds = np.dot(pos_earth.T, src_unitvector)/CLIGHT_KPS
correction_days = correction_seconds/SEC_P_DAY
elif jd_type == 'hjd':
# Compute HJD correction via Sun ephemeris
# this is the position vector of the center of the sun in km
# Solar System Barycenter (0) -> Sun (10)
pos_sun = jplkernel[0,10].compute(t.tdb.jd)
# this is the vector from the center of the sun to the center of the
# earth
sun_earth_vec = pos_earth - pos_sun
# calculate the heliocentric correction
correction_seconds = np.dot(sun_earth_vec.T, src_unitvector)/CLIGHT_KPS
correction_days = correction_seconds/SEC_P_DAY
# TDB is the appropriate time scale for these ephemerides
new_jd = t.tdb.jd + correction_days
return new_jd