Added a date module for julian and ls conversion

plio/date/astrodate.py

+""" 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
+

plio/date/julian2ls.py

+#!/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))
+

plio/date/julian2season.py

+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
+

plio/date/tests/test_julian2ls.py

+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()