3

Somewhat new to python but trying to use it for satellite orbit analysis. I'm really not getting the datetime object and methods. I want to be able to pass a datetime object to a function that accepts arguments in essentially the same format that datetime does (year, mon, day, hour, min, sec). The code below works but there has to be a better way. Please enlighten me. Thanks!

import jday
import datetime

jdate = jday.JD(2015,12,1,22,8,0) # example

now1 = datetime.datetime.now().strftime("%Y,%m,%d,%H,%M,%S")
now2 = now1.split(",")
now3 = [int(i) for i in now2]
jdatenow = jday.JD(*now3)
print jdatenow

jday module is ported over from Matlab of David Vallado's Astrodynamics source code.

import math as m

def JD(yr, mon, day, hr, min, sec):
        jd = 367.0 * yr - m.floor( 
        (7 * (yr + m.floor( (mon + 9) / 12.0) ) ) * 0.25 ) + m.floor(
        275 * mon / 9.0 ) + day + 1721013.5 + (
        (sec/60.0 + min ) / 60.0 + hr ) / 24.0
        return jd
jfs
  • 399,953
  • 195
  • 994
  • 1,670
Jesse Reich
  • 87
  • 1
  • 1
  • 11
  • would an inline version be better for you? like `now=[int(i) for i in datetime.datetime.now().strftime("%Y,%m,%d,%H,%M,%S").split(",")]` – Antonio Ragagnin Dec 04 '15 at 17:20
  • unrelated: the formula is wrong by a whole day in 1900 and 2100, see [details in my answer](http://stackoverflow.com/a/34092799/4279) – jfs Dec 09 '15 at 13:09

4 Answers4

2

Given you have ported over the JD code, and therefore have control over it as module jday, perhaps a decorator is what you are looking for. This has the obvious benefit of not breaking the original signature of the function (for existing client code), but adds the convenience of a date param as you requested.

Have also created a jday2 module which is the same as the original jday module but whose JD() function instead accepts a date object directly. This is the simplest solution if you don't need any backwards compatibility.

Please see working example code below:

jday.py

import math as m
import functools

def date_support_wrapper(f):
    """ Wraps JD and provides a way to pass a date param

    :param f: the original function
    :return: the wrapper around the original function
    """

    @functools.wraps(f)
    def wrap(*args, **kwargs):
        if 'date' in kwargs:
            d = kwargs['date']
            return f(yr=d.year, mon=d.month, day=d.day, hr=d.hour, min=d.minute, sec=d.second)

        return f(*args, **kwargs)

    return wrap


@date_support_wrapper
def JD(yr, mon, day, hr, min, sec):
    jd = 367.0 * yr - m.floor(
        (7 * (yr + m.floor((mon + 9) / 12.0))) * 0.25) + m.floor(
        275 * mon / 9.0) + day + 1721013.5 + (
                                                 (sec / 60.0 + min) / 60.0 + hr) / 24.0
    return jd

jday2.py

import math as m

def JD(dt):
    """ Same calculation as JD but accepts a date object argument

    :param dt: the date object
    :return: the JD result
    """
    yr, mon, day, hr, min, sec = dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second
    jd = 367.0 * yr - m.floor(
        (7 * (yr + m.floor((mon + 9) / 12.0))) * 0.25) + m.floor(
        275 * mon / 9.0) + day + 1721013.5 + (
                                                 (sec / 60.0 + min) / 60.0 + hr) / 24.0
    return jd

And example client code:

client.py

import datetime
import jday
import jday2

# The date we are interested in
a = dict(year=2015, month=12, day=1, hour=22, minute=8, second=0)
dt = datetime.datetime(**a)  # 2015-12-01 22:08:00

# The original signature of the function
jdate1 = jday.JD(a['year'], a["month"], a["day"], a["hour"], a["minute"], a["second"])
# 2457358.422222222

# The new signature that accepts a normal date object
# Note that we use keyword "date" argument
jdate2 = jday.JD(date=dt)
# 2457358.422222222

# The new signature that accepts a normal date object
jdate3 = jday2.JD(dt)
# 2457358.422222222
arcseldon
  • 35,523
  • 17
  • 121
  • 125
  • I'll have to read more about the the decorator concept but I think this is what I was hoping to do. Change the jday function to accept the date object. Thanks. – Jesse Reich Dec 03 '15 at 02:25
  • @JesseReich - please can you mark this as the accepted answer :) Decorators are a great way of using "aspect oriented programming" (AOP) to add / modify existing behaviour of a given function without actually altering that function itself. It just depends what you need - if you have no reason to support both function argument types, then simply modify the existing JD signature to accept a date object instead, and parse out the arguments you want from that date object. – arcseldon Dec 03 '15 at 05:49
  • @JesseReich - have just updated the answer (eventually..) to include a separate jday2.py module which may be all you need if you consider decorators too much for your needs. – arcseldon Dec 03 '15 at 06:14
1

Unfortunately you can't pass a datetime object directly to a function that expects integers. Is there any particular reason you can't call the function with the fields datetime provides?

now = datetime.now()
jdatenow = jday.JD(now.year, now.month, now.day, now.hour, now.minute, now.second)

Or if that's too bothersome, wrap it up:

from datetime import datetime

def my_JD(dt):
    return jday.JD(dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second)
jdatenow = my_JD(datetime.now())

This may be easier to understand than

jdatenow = jday.JD(*list(datetime.now().timetuple())[:-3])
Aaron B
  • 859
  • 11
  • 17
1

You could rewrite it without math.floor():

def jdate(year, month, day, hour, minute, second):
    day_fraction = ((second + 60 * minute) + 3600 * hour) / 86400.
    return (367 * year - (7 * (year + (month + 9) // 12)) // 4 +
            275 * month // 9 + day + (1721013.5 + day_fraction))

You could simplify it; if you use datetime arithmetic:

#!/usr/bin/env python3
from datetime import datetime, timedelta

DAY = timedelta(1)
JULIAN_EPOCH = datetime(2000, 1, 1, 12)  # noon (the epoch name is unrelated)
J2000_JD = timedelta(2451545)  # julian epoch in julian dates

def JD(dt):
    """Julian Date: JD(UTC)."""
    return (dt - JULIAN_EPOCH + J2000_JD) / DAY

To pass a datetime object to jdate(), you could use .timetuple() method:

import math

for time_tuple in [(1961, 1, 1), (1968, 2, 1), (1972, 1, 1), (1996, 1, 1)]:
    dt = datetime(*time_tuple)
    a, b = jdate(*dt.timetuple()[: 6]), JD(dt)
    print("{} UTC -> {} JD(UTC)".format(dt, b))
    assert math.isclose(a, b), (a, b)

Also, you could use dt.year, dt.month, dt.day, etc attributes if necessary.

Output

1961-01-01 00:00:00 UTC -> 2437300.5 JD(UTC)
1968-02-01 00:00:00 UTC -> 2439887.5 JD(UTC)
1972-01-01 00:00:00 UTC -> 2441317.5 JD(UTC)
1996-01-01 00:00:00 UTC -> 2450083.5 JD(UTC)

That is correct according to IERS web-site where a recommended "Julian Date" definition is provided.


The formulas produce different results for dates before March 1900 and after February 2100:

import jdcal  # pip install jdcal
import astropy.time  # pip install astropy

print("                UTC |    matlab |  datetime |   astropy | jdcal")
for year in [1900, 2000, 2100]:
    for time_tuple in [(year, 2, 28, 12), (year, 3, 1, 12)]:
        dt = datetime(*time_tuple)
        matlabJD = jdate(*dt.timetuple()[:6])
        datetimeJD = JD(dt)
        jdcalJD = sum(jdcal.gcal2jd(*dt.timetuple()[:3])) + .5
        astropyJD = astropy.time.Time(dt)
        print("{dt} | {matlabJD} | {datetimeJD} | {astropyJD.jd} | {jdcalJD}"
              .format(**vars()))

Output

                UTC |    matlab |  datetime |   astropy | jdcal
1900-02-28 12:00:00 | 2415078.0 | 2415079.0 | 2415079.0 | 2415079.0
1900-03-01 12:00:00 | 2415080.0 | 2415080.0 | 2415080.0 | 2415080.0
2000-02-28 12:00:00 | 2451603.0 | 2451603.0 | 2451603.0 | 2451603.0
2000-03-01 12:00:00 | 2451605.0 | 2451605.0 | 2451605.0 | 2451605.0
2100-02-28 12:00:00 | 2488128.0 | 2488128.0 | 2488128.0 | 2488128.0
2100-03-01 12:00:00 | 2488130.0 | 2488129.0 | 2488129.0 | 2488129.0

jdate() formula in your question thinks that 1900, 2100 are leap years. The datetime implementation, astropy and jdcal libraries produce the same results here.

Note: Julian day is an integer. JD() computes a Julian date that includes the fraction of the day, see definitions in the links.

As mentioned in the linked discussion, you should use an already made libraries and send patches if necessary instead of reinventing the wheel, to avoid simple errors due to leap years, floating points issues, wrong time scales, slight difference in the definitions for the term Julian day, etc.

Community
  • 1
  • 1
jfs
  • 399,953
  • 195
  • 994
  • 1,670
  • In satellite orbital analysis a fraction of a second could be significant, and there are several time scales that might be chosen depending on the details of the calculation. Output from whatever solution is ultimately chosen should properly label the time scale. A simple label of "UTC" would arouse suspicion in those receiving the output since various software libraries that attempt to support Julian dates have treated leap seconds in different ways. – Gerard Ashton Dec 04 '15 at 21:06
  • @GerardAshton: UTC means UTC in my answer. Click the link with IERS in the answer. – jfs Dec 04 '15 at 21:10
  • AJ.F. Sebastian: [An IAU source] http://www.iausofa.org/sofa_ts_c.pdf indicates UTC was initiated in 1960 so output labeling date/times in 1900 as UTC is suspicious on its face. – Gerard Ashton Dec 04 '15 at 21:24
  • @GerardAshton: the purpose of 1900 is to compare different implementations: the difference is a ***whole day*** that is much much larger any difference between time scales (such as UT1, TAI). – jfs Dec 04 '15 at 21:42
  • @GerardAshton: if you need a subsecond precision then you should probably use a uniform time scale (such as GPS) as recommended in the IERS link. UTC is non-uniform (due to leap seconds). – jfs Dec 04 '15 at 21:55
0

The shortest way I could figure out:

now = list(datetime.datetime.now().timetuple())[:-3]

Alex
  • 963
  • 9
  • 11