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.