diff --git a/plio/date/__init__.py b/plio/date/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/plio/date/astrodate.py b/plio/date/astrodate.py new file mode 100644 index 0000000000000000000000000000000000000000..69a03867b9a4808df761e84aaccbacaee8f34d4a --- /dev/null +++ b/plio/date/astrodate.py @@ -0,0 +1,314 @@ +""" For more information about astronomical date specifications, +consult a reference source such as +U{this page <http://tycho.usno.navy.mil/systime.html>} provided by +the US Naval Observatory. + +Constants and formulae in this module were taken from the times.h +include file of the tpm package by Jeff Percival, to ensure compatibility. + +""" + +import types +import datetime + +#first define some globally-useful constants +B1950 = 2433282.42345905 +J2000 = 2451545.0 +CB = 36524.21987817305 +CJ = 36525.0 +MJD_0 = 2400000.5 + +def jyear2jd(jyear): + """ + @param jyear: decimal Julian year + @type jyear: float + @return: Julian date + @rtype: float + """ + return (J2000 + ((jyear)-2000.0)*(CJ/100.0)) + +def jd2jyear(jd): + """ + @return: decimal Julian year + @rtype: float + @param jd: Julian date + @type jd: float + """ + return (2000.0 + ((jd)-J2000)*(100.0/CJ)) + +def byear2jd(byear): + """ + @param byear: decimal Besselian year + @type byear: float + @return: Julian date + @rtype: float + """ + return (B1950 + ((byear)-1950.0)*(CB/100.0)) + +def utc2jd(utc): + """ + Convert UTC to Julian date. + + Conversion translated from TPM modules utcnow.c and gcal2j.c, which + notes that the algorithm to convert from a gregorian proleptic calendar + date onto a julian day number is taken from + The Explanatory Supplement to the Astronomical Almanac (1992), + section 12.92, equation 12.92-1, page 604. + + @param utc: UTC (Universal Civil Time) + @type utc: U{datetime<http://docs.python.org/lib/datetime.html>} object + @return: Julian date (to the nearest second) + @rtype: float + + + + """ + +## But, beware of different rounding behavior between C and Python! +## - integer arithmetic truncates -1.07 to -2 in Python; to -1 in C +## - to reproduce C-like behavior in Python, do the math with float +## arithmetic, then explicitly cast to int. + + + y=float(utc.year) + m=float(utc.month) + d=float(utc.day) + hr=utc.hour + min=utc.minute + sec=utc.second + + #Address differences between python and C time conventions + # C: Python datetime + # 0 <= mon <= 11 1 <= month <= 12 + # + + #C code to get the julian date of the start of the day */ + #takes as input 1900+ptm->tm_year, ptm->tm_mon+1, ptm->tm_mday + # So we can use just (year, month, mday) + + mterm=int((m-14)/12) + aterm=int((1461*(y+4800+mterm))/4) + + bterm=int((367*(m-2-12*mterm))/12) + + cterm=int((3*int((y+4900+mterm)/100))/4) + + j=aterm+bterm-cterm+d + j -= 32075 + #offset to start of day + j -= 0.5 + + +# print "h/m/s: %f/%f/%f"%(hr,min,sec) + + #Apply the time + jd = j + (hr + (min + (sec/60.0))/60.0)/24.0 + + return jd + +def AstroDate(datespec=None): + """AstroDate can be used as a class for managing astronomical + date specifications (despite the fact that it was implemented + as a factory function) that returns either a BesselDate or a + JulianDate, depending on the properties of the datespec. + + AstroDate was originally conceived as a Helper class for the + Position function for use with pytpm functionality, but also as a + generally useful class for managing astronomical date specifications. + + The philosophy is the same as Position: to enable the user to specify + the date once and for all, and access it in a variety of styles. + + @param datespec: Date specification as entered by the user. Permissible + specifications include: + - Julian year: 'J1997', 'J1997.325', 1997.325: return a JulianDate + - Besselian year: 'B1950','B1958.432': return a BesselDate + - Julian date: 'JD2437241.81', '2437241.81', 2437241.81: return a JulianDate + - Modified Julian date: 'MJD37241.31': returns a JulianDate + - A U{datetime <http://docs.python.org/lib/datetime.html>} object: return a JulianDate (assumes input time is UTC) + - None: returns the current time as a JulianDate + + @type datespec: string, float, integer, U{datetime <http://docs.python.org/lib/datetime.html>}, or None + + @rtype: L{JulianDate} or L{BesselDate} + + @raise ValueError: Raises an exception if the date specification is a + string, but begins with a letter that is not 'B','J','JD', or 'MJD' + (case insensitive). + + @todo: Add math functions! Addition, subtraction. + @todo: Is there a need to support other date specifications? + eg FITS-style dates? + """ + + if datespec is None: + return JulianDate(datetime.datetime.utcnow()) + + try: + dstring=datespec.upper() + if dstring.startswith('B'): + #it's a Besselian date + return BesselDate(datespec) + elif ( dstring.startswith('JD') or + dstring.startswith('MJD') or + dstring.startswith('J') ): + #it's Julian, one way or another + return JulianDate(datespec) + elif dstring[0].isalpha(): + raise ValueError, "Invalid system specification: must be B, J, JD, or MJD" + else: #it must be a numeric string: assume Julian + return JulianDate(datespec) + + except AttributeError: + #if no letter is specified, assume julian + return JulianDate(datespec) + + + +class JulianDate: + """ + @ivar year: Decimal Julian year + @type year: float + + @ivar jd: Julian date + @type jd: float + + @ivar mjd: Modified Julian Date + @type mjd: float + + @ivar datespec: Date specification as entered by the user + """ + + def __init__(self,datespec): + self.datespec=datespec + + if isinstance(datespec,datetime.datetime): + #Note assumption that datespec is already in UTC + self.jd=utc2jd(datespec) + self.mjd=self.jd-MJD_0 + self.year=jd2jyear(self.jd) + + elif type(datespec) is types.StringType: + if datespec.upper().startswith('JD'): + #it's a julian date + self.jd=float(datespec[2:]) + self.mjd=self.jd-MJD_0 + self.year=jd2jyear(self.jd) + elif datespec.upper().startswith('MJD'): + #it's a modified julian date + self.mjd=float(datespec[3:]) + self.jd=self.mjd+MJD_0 + self.year=jd2jyear(self.jd) + elif datespec.upper().startswith('J'): + #it's a julian decimal year + self.year=float(datespec[1:]) + self.jd=jyear2jd(self.year) + self.mjd=self.jd-MJD_0 + + elif not datespec[0].isalpha(): + #somebody put a numeric date in quotes + datespec=float(datespec) + if datespec < 10000: + #it's a year. Assume julian. + self.year=float(datespec) + self.jd=jyear2jd(self.year) + self.mjd=self.jd-MJD_0 + else: + #it's a date. Assume JD not MJD. + self.jd=float(datespec) + self.mjd=self.jd+MJD_0 + self.year=jd2jyear(self.jd) + else: + print "help, we are confused" + + else: #it's a number + if datespec < 10000: + #it's a year. Assume julian. + self.year=float(datespec) + self.jd=jyear2jd(self.year) + self.mjd=self.jd-MJD_0 + else: + #it's a date. Assume JD not MJD. + self.jd=datespec + self.mjd=self.jd+MJD_0 + self.year=jd2jyear(self.jd) + + def __repr__(self): + return str(self.datespec) + + def __equals__(self,other): + """ All comparisons will be done based on jd attribute """ + return self.jd == other.jd + + def __lt__(self,other): + return self.jd < other.jd + + def __gt__(self,other): + return self.jd > other.jd + + def __le__(self,other): + return self.jd <= other.jd + + def __ge__(self,other): + return self.jd >= other.jd + + def byear(self): + """ Return Besselian year based on previously calculated + julian date. + @return: decimal Besselian year + @rtype: float + """ + ans=(1950.0 + ((x)-B1950)*(100.0/CB)) + return ans + + +class BesselDate: + """ + @ivar year: Decimal Besselian year + @type year: float + + @ivar jd: Julian date + @type jd: float + + @ivar mjd: Modified Julian Date + @type mjd: float + + @ivar datespec: Date specification as entered by the user + + """ + + def __init__(self,datespec): + self.datespec=datespec + try: + self.year=float(datespec) + except ValueError: + self.year=float(datespec[1:]) + self.jd=byear2jd(self.year) + self.mjd=self.jd-MJD_0 + + def __equals__(self,other): + """ All comparisons will be done based on jd attribute """ + return self.jd == other.jd + + def __lt__(self,other): + return self.jd < other.jd + + def __gt__(self,other): + return self.jd > other.jd + + def __le__(self,other): + return self.jd <= other.jd + + def __ge__(self,other): + return self.jd >= other.jd + + + def jyear(self): + """ Return the julian year using the already-converted + julian date + @return: Decimal Julian year + @rtype: float + """ + ans = jd2jyear(self.jd) + return ans + diff --git a/plio/date/julian2ls.py b/plio/date/julian2ls.py new file mode 100644 index 0000000000000000000000000000000000000000..963f9a3f05924331d80f07aebd006625fb76e69a --- /dev/null +++ b/plio/date/julian2ls.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python +import sys +import math +import numpy as np + + +def zero360(angles, rad=False): + """ + angles: ndarray: a scalar or vector of float angles + rad: boolean: flag whether angles are in radians + """ + + if rad: + half = math.pi + else: + half = 180.0 + + circle = 2.0 * half + + ii = np.floor(angles / circle) + bb = angles - circle * ii + #Negative values to positive + bb[bb < 0] += circle + return bb + + +def julian2ls(date, marsyear=None, reverse=False): + """ + date: Scalar or NumPy ndarray + marsyear: + reverse: Boolean + Reverse conversion from L_{s} to julian + + Original IDL from Hugh Keiffer + + References: + ----------- + [1] M. Allison and M. McEwen.' A post-Pathfinder evaluation of areocentric + solar coordinates with improved timing recipes for Mars seasonal/diurnal + climate studies'. Plan. Space Sci., 2000, v=48, pages = {215-235}, + + [2] http://www.giss.nasa.gov/tools/mars24/help/algorithm.html + + """ + + if not isinstance(date, np.ndarray): + date = np.asarray([date], dtype=np.float64) + + rec= np.array([-0.0043336633,0.0043630568,0.0039601861,0.029025116,-0.00042300026]) + dj2000 = 2451545.0 # JD of epoch J2000 + dj4m = 51544.50 # days offset from dj4 to djm + smja = 1.523679 # semi-major axis in AU, Table 2. last digit from Ref2 + meanmotion = 0.52402075 # mean motion, degrees/day. Constant in Table 2 + + a5 = np.array([71,57,39,37,21,20,18]) * 0.0001 + tau = np.array([2.2353,2.7543,1.1177,15.7866,2.1354,2.4694,32.8493]) + phi = np.array([49.409,168.173,191.837,21.736,15.704,95.528,49.095]) + + + if reverse: + #Convert from lsubs to julian + lsubs = date + delta_lsubs = lsubs - 250.99864 + rdelta_lsubs= np.radians(delta_lsubs) + dj4 = 51507.5 + 1.90826 * delta_lsubs -\ + 20.42 * np.sin(rdelta_lsubs) +\ + 0.72 * np.sin( 2. * rdelta_lsubs) # Eq. 14 + + brak = 686.9726 + 0.0043 * np.cos(rdelta_lsubs) -\ + 0.0003 * np.cos(2. * rdelta_lsubs) # in brackets in Eq 14 + if marsyear is None: + return False + ny=marsyear + 42 # orbits of mars since 1874.0 + dj4 += brak * (ny - 66) # last part of Eq. 14 + djtt= dj4 - dj4m # days from J2000 TT + tcen = djtt / 36525.0 # julian centuries from J2000 + tcor= 64.184 + tcen * (59.0 + tcen * (-51.2 + tcen * (-67.1 - tcen * 16.4))) # TT-UTC + djm=djtt - tcor / 86400.0 # convert correction to days and apply + out=djm # out is UTC + else: + #Convert from julian to lsubs + if date[0] < 1.e6: #Is this the intended functionality - what if date[0] < and date [1] > + djm = date + else: + djm = date - dj2000 + tcen = djm / 36525.0 # julian centuries from J2000 + tcor = 64.184 + tcen * (59.0 + tcen * (-51.2 + tcen * (-67.1 - tcen * 16.4))) # TT-UTC + djtt = djm + tcor / 86400.0 # convert correction to days and apply. get TT + dj4 = djtt + dj4m # A+M's MJD + nelements = date.size + pbs = np.zeros(nelements) + for i in range(nelements): + q = 0.985626 * djtt[i] + rq = np.radians(q / tau + phi) + pbsi = a5 * np.cos(rq) + pbs[i] = pbsi.sum() + meananomoly = 19.3870 + meanmotion * djtt # Mean anomoly M, Table 2 and Eq. 16 + rmeananomoly = np.radians(meananomoly) # M in radians + #Eqn of center, in brackets in eq 20 and all but the first term in eq 19. exclude pbs + eoc = (10.691 + 3.e-7 * djtt) * np.sin(rmeananomoly)\ + + 0.623 * np.sin(2. * rmeananomoly) + 0.050 * np.sin(3. * rmeananomoly)\ + + 0.005 * np.sin(4. * rmeananomoly) + 0.0005 * np.sin(5. * rmeananomoly) + + if reverse: + try: + j = reverse.size + except: + j=0 + rev = [0] + if j >= 5: + #User has supplied coefficients + rec = rev + if j < 5 and rev[0] < 5: + fd = 0.0 + else: + fd=rec[0] + rec[1]* np.sin(rdelta_lsubs)\ + + rec[2] * np.sin(2. * rdelta_lsubs)\ + + rec[3] * np.sin(3. * rdelta_lsubs)\ + + rec[4] * np.sin(4. * rdelta_lsubs) + #TODO: This looks like reversing is going to cause an fd undefined error if j >= 5 + out = out - ( pbs / meananomoly + fd) + else: + afms = 270.3863 + 0.52403840 * djtt # Eq. 17 + lsubs = afms + eoc + pbs # Eq. 19 LS in degrees + lsubs = zero360(lsubs) + out= lsubs + marstyear = 686.9728 # mean mars tropical year in terrestrial days + #marsysol=668.5991 #mars siderial year in sols + ny = np.floor((dj4 - 5668.690)/ marstyear ) # full Mars years from 1874 + myn = ny - 42 # climate MY + + if reverse: + return out + else: + return out, myn + +if __name__ == '__main__': + print(julian2ls(178, marsyear=40, reverse=True)) + diff --git a/plio/date/julian2season.py b/plio/date/julian2season.py new file mode 100644 index 0000000000000000000000000000000000000000..609ce8674324a187701756117ed20c8764af6421 --- /dev/null +++ b/plio/date/julian2season.py @@ -0,0 +1,28 @@ +import krc.config as config + +def j2season(_date): + """ + Ls date to a KRC season. + + Parameters + ----------- + _date : float + The input date to be converted + + Returns + ------- + startseason : int + The integer index to the start season + stopseasons : int + The integer index to the stop season + """ + date = _date + if date < config.RECORD_START_DATE: + remainder = (config.RECORD_START_DATE - date) / config.YEAR + date = date + int(remainder + 1.0) * config.YEAR + dateoffset = (date - config.RECORD_START_DATE) % config.YEAR + recordoffset = dateoffset / config.MARTIAN_DAY + startseason = int(recordoffset) + stopseason = startseason + 1 + return recordoffset, startseason, stopseason + diff --git a/plio/date/tests/__init__.py b/plio/date/tests/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/plio/date/tests/test_julian2ls.py b/plio/date/tests/test_julian2ls.py new file mode 100644 index 0000000000000000000000000000000000000000..c520df9b09fb8c7e8c577db69f39b253acf266d7 --- /dev/null +++ b/plio/date/tests/test_julian2ls.py @@ -0,0 +1,58 @@ +import numpy as np +import unittest + +from .. import julian2ls + + +class TestSmallJulian(unittest.TestCase): + def runTest(self): + smalldates = np.array([10834, 65432, 11111, 104]) + expected = np.array([177.77841982441714, 9.1150791399559239, + 343.68773181354481, 335.45117593903956]) + expectedmyn = np.array([40, 120, 40, 24]) + result, myn = julian2ls.julian2ls(smalldates) + np.testing.assert_allclose(expected, result, rtol=1e-4) + np.testing.assert_allclose(expectedmyn, myn) + + +class TestLargeJulian(unittest.TestCase): + + def runTest(self): + largedates = np.array([1100354, 1260543, 2020511, 1986287]) + expected = np.array([138.01380245992914, 289.56588061130606, + 103.53916196653154, 165.34239548092592]) + expectedmyn = np.array([-1943, -1710, -603, -653]) + result, myn = julian2ls.julian2ls(largedates) + np.testing.assert_allclose(expected, result, rtol=1e-4) + np.testing.assert_allclose(expectedmyn, myn) + + +class TestReverse(unittest.TestCase): + + def runTest(self): + expected = np.array([10834.400119543538]) + date = julian2ls.julian2ls(178, marsyear=40, reverse=True) + np.testing.assert_allclose(date, expected, rtol=1e-4) + +class TestLanders(unittest.TestCase): + """ + Testing against: http://www-mars.lmd.jussieu.fr/mars/time/martian_time.html + """ + def test_year_start(self): + result, myn = julian2ls.julian2ls(2435198.5) + self.assertEqual(myn, 0) + np.testing.assert_array_almost_equal(result, 354.743, 3) + + def test_curiosity(self): + result, myn = julian2ls.julian2ls(2456145.7207986116) + print(result) + np.testing.assert_array_almost_equal(result, 150.702, 3) + + def test_pathfinder(self): + result, myn = julian2ls.julian2ls(2450634.2061921293) + self.assertEqual(myn, 23) + np.testing.assert_array_almost_equal(result, 142.725, 3) + +if __name__ == '__main__': + + unittest.main()