jdcal.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525
  1. # -*- coding:utf-8 -*-
  2. """Functions for converting between Julian dates and calendar dates.
  3. A function for converting Gregorian calendar dates to Julian dates, and
  4. another function for converting Julian calendar dates to Julian dates
  5. are defined. Two functions for the reverse calculations are also
  6. defined.
  7. Different regions of the world switched to Gregorian calendar from
  8. Julian calendar on different dates. Having separate functions for Julian
  9. and Gregorian calendars allow maximum flexibility in choosing the
  10. relevant calendar.
  11. All the above functions are "proleptic". This means that they work for
  12. dates on which the concerned calendar is not valid. For example,
  13. Gregorian calendar was not used prior to around October 1582.
  14. Julian dates are stored in two floating point numbers (double). Julian
  15. dates, and Modified Julian dates, are large numbers. If only one number
  16. is used, then the precision of the time stored is limited. Using two
  17. numbers, time can be split in a manner that will allow maximum
  18. precision. For example, the first number could be the Julian date for
  19. the beginning of a day and the second number could be the fractional
  20. day. Calculations that need the latter part can now work with maximum
  21. precision.
  22. A function to test if a given Gregorian calendar year is a leap year is
  23. defined.
  24. Zero point of Modified Julian Date (MJD) and the MJD of 2000/1/1
  25. 12:00:00 are also given.
  26. This module is based on the TPM C library, by Jeffery W. Percival. The
  27. idea for splitting Julian date into two floating point numbers was
  28. inspired by the IAU SOFA C library.
  29. :author: Prasanth Nair
  30. :contact: prasanthhn@gmail.com
  31. :license: BSD (http://www.opensource.org/licenses/bsd-license.php)
  32. """
  33. from __future__ import division
  34. from __future__ import print_function
  35. import math
  36. __version__ = "1.0"
  37. MJD_0 = 2400000.5
  38. MJD_JD2000 = 51544.5
  39. def fpart(x):
  40. """Return fractional part of given number."""
  41. return math.modf(x)[0]
  42. def ipart(x):
  43. """Return integer part of given number."""
  44. return math.modf(x)[1]
  45. def is_leap(year):
  46. """Leap year or not in the Gregorian calendar."""
  47. x = math.fmod(year, 4)
  48. y = math.fmod(year, 100)
  49. z = math.fmod(year, 400)
  50. # Divisible by 4 and,
  51. # either not divisible by 100 or divisible by 400.
  52. return not x and (y or not z)
  53. def gcal2jd(year, month, day):
  54. """Gregorian calendar date to Julian date.
  55. The input and output are for the proleptic Gregorian calendar,
  56. i.e., no consideration of historical usage of the calendar is
  57. made.
  58. Parameters
  59. ----------
  60. year : int
  61. Year as an integer.
  62. month : int
  63. Month as an integer.
  64. day : int
  65. Day as an integer.
  66. Returns
  67. -------
  68. jd1, jd2: 2-element tuple of floats
  69. When added together, the numbers give the Julian date for the
  70. given Gregorian calendar date. The first number is always
  71. MJD_0 i.e., 2451545.5. So the second is the MJD.
  72. Examples
  73. --------
  74. >>> gcal2jd(2000,1,1)
  75. (2400000.5, 51544.0)
  76. >>> 2400000.5 + 51544.0 + 0.5
  77. 2451545.0
  78. >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678.0, 2000,
  79. ....: 2012, 2245]
  80. >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
  81. >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
  82. >>> x = [gcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
  83. >>> for i in x: print i
  84. (2400000.5, -2395215.0)
  85. (2400000.5, -1451021.0)
  86. (2400000.5, -1062364.0)
  87. (2400000.5, -723762.0)
  88. (2400000.5, -679162.0)
  89. (2400000.5, -678774.0)
  90. (2400000.5, -678368.0)
  91. (2400000.5, -633797.0)
  92. (2400000.5, -65812.0)
  93. (2400000.5, 51827.0)
  94. (2400000.5, 56242.0)
  95. (2400000.5, 141393.0)
  96. Negative months and days are valid. For example, 2000/-2/-4 =>
  97. 1999/+12-2/-4 => 1999/10/-4 => 1999/9/30-4 => 1999/9/26.
  98. >>> gcal2jd(2000, -2, -4)
  99. (2400000.5, 51447.0)
  100. >>> gcal2jd(1999, 9, 26)
  101. (2400000.5, 51447.0)
  102. >>> gcal2jd(2000, 2, -1)
  103. (2400000.5, 51573.0)
  104. >>> gcal2jd(2000, 1, 30)
  105. (2400000.5, 51573.0)
  106. >>> gcal2jd(2000, 3, -1)
  107. (2400000.5, 51602.0)
  108. >>> gcal2jd(2000, 2, 28)
  109. (2400000.5, 51602.0)
  110. Month 0 becomes previous month.
  111. >>> gcal2jd(2000, 0, 1)
  112. (2400000.5, 51513.0)
  113. >>> gcal2jd(1999, 12, 1)
  114. (2400000.5, 51513.0)
  115. Day number 0 becomes last day of previous month.
  116. >>> gcal2jd(2000, 3, 0)
  117. (2400000.5, 51603.0)
  118. >>> gcal2jd(2000, 2, 29)
  119. (2400000.5, 51603.0)
  120. If `day` is greater than the number of days in `month`, then it
  121. gets carried over to the next month.
  122. >>> gcal2jd(2000,2,30)
  123. (2400000.5, 51604.0)
  124. >>> gcal2jd(2000,3,1)
  125. (2400000.5, 51604.0)
  126. >>> gcal2jd(2001,2,30)
  127. (2400000.5, 51970.0)
  128. >>> gcal2jd(2001,3,2)
  129. (2400000.5, 51970.0)
  130. Notes
  131. -----
  132. The returned Julian date is for mid-night of the given date. To
  133. find the Julian date for any time of the day, simply add time as a
  134. fraction of a day. For example Julian date for mid-day can be
  135. obtained by adding 0.5 to either the first part or the second
  136. part. The latter is preferable, since it will give the MJD for the
  137. date and time.
  138. BC dates should be given as -(BC - 1) where BC is the year. For
  139. example 1 BC == 0, 2 BC == -1, and so on.
  140. Negative numbers can be used for `month` and `day`. For example
  141. 2000, -1, 1 is the same as 1999, 11, 1.
  142. The Julian dates are proleptic Julian dates, i.e., values are
  143. returned without considering if Gregorian dates are valid for the
  144. given date.
  145. The input values are truncated to integers.
  146. """
  147. year = int(year)
  148. month = int(month)
  149. day = int(day)
  150. a = ipart((month - 14) / 12.0)
  151. jd = ipart((1461 * (year + 4800 + a)) / 4.0)
  152. jd += ipart((367 * (month - 2 - 12 * a)) / 12.0)
  153. x = ipart((year + 4900 + a) / 100.0)
  154. jd -= ipart((3 * x) / 4.0)
  155. jd += day - 2432075.5 # was 32075; add 2400000.5
  156. jd -= 0.5 # 0 hours; above JD is for midday, switch to midnight.
  157. return MJD_0, jd
  158. def jd2gcal(jd1, jd2):
  159. """Julian date to Gregorian calendar date and time of day.
  160. The input and output are for the proleptic Gregorian calendar,
  161. i.e., no consideration of historical usage of the calendar is
  162. made.
  163. Parameters
  164. ----------
  165. jd1, jd2: int
  166. Sum of the two numbers is taken as the given Julian date. For
  167. example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
  168. can be the MJD of the date and time. But any combination will
  169. work.
  170. Returns
  171. -------
  172. y, m, d, f : int, int, int, float
  173. Four element tuple containing year, month, day and the
  174. fractional part of the day in the Gregorian calendar. The first
  175. three are integers, and the last part is a float.
  176. Examples
  177. --------
  178. >>> jd2gcal(*gcal2jd(2000,1,1))
  179. (2000, 1, 1, 0.0)
  180. >>> jd2gcal(*gcal2jd(1950,1,1))
  181. (1950, 1, 1, 0.0)
  182. Out of range months and days are carried over to the next/previous
  183. year or next/previous month. See gcal2jd for more examples.
  184. >>> jd2gcal(*gcal2jd(1999,10,12))
  185. (1999, 10, 12, 0.0)
  186. >>> jd2gcal(*gcal2jd(2000,2,30))
  187. (2000, 3, 1, 0.0)
  188. >>> jd2gcal(*gcal2jd(-1999,10,12))
  189. (-1999, 10, 12, 0.0)
  190. >>> jd2gcal(*gcal2jd(2000, -2, -4))
  191. (1999, 9, 26, 0.0)
  192. >>> gcal2jd(2000,1,1)
  193. (2400000.5, 51544.0)
  194. >>> jd2gcal(2400000.5, 51544.0)
  195. (2000, 1, 1, 0.0)
  196. >>> jd2gcal(2400000.5, 51544.5)
  197. (2000, 1, 1, 0.5)
  198. >>> jd2gcal(2400000.5, 51544.245)
  199. (2000, 1, 1, 0.24500000000261934)
  200. >>> jd2gcal(2400000.5, 51544.1)
  201. (2000, 1, 1, 0.099999999998544808)
  202. >>> jd2gcal(2400000.5, 51544.75)
  203. (2000, 1, 1, 0.75)
  204. Notes
  205. -----
  206. The last element of the tuple is the same as
  207. (hh + mm / 60.0 + ss / 3600.0) / 24.0
  208. where hh, mm, and ss are the hour, minute and second of the day.
  209. See Also
  210. --------
  211. gcal2jd
  212. """
  213. from math import modf
  214. jd1_f, jd1_i = modf(jd1)
  215. jd2_f, jd2_i = modf(jd2)
  216. jd_i = jd1_i + jd2_i
  217. f = jd1_f + jd2_f
  218. # Set JD to noon of the current date. Fractional part is the
  219. # fraction from midnight of the current date.
  220. if -0.5 < f < 0.5:
  221. f += 0.5
  222. elif f >= 0.5:
  223. jd_i += 1
  224. f -= 0.5
  225. elif f <= -0.5:
  226. jd_i -= 1
  227. f += 1.5
  228. l = jd_i + 68569
  229. n = ipart((4 * l) / 146097.0)
  230. l -= ipart(((146097 * n) + 3) / 4.0)
  231. i = ipart((4000 * (l + 1)) / 1461001)
  232. l -= ipart((1461 * i) / 4.0) - 31
  233. j = ipart((80 * l) / 2447.0)
  234. day = l - ipart((2447 * j) / 80.0)
  235. l = ipart(j / 11.0)
  236. month = j + 2 - (12 * l)
  237. year = 100 * (n - 49) + i + l
  238. return int(year), int(month), int(day), f
  239. def jcal2jd(year, month, day):
  240. """Julian calendar date to Julian date.
  241. The input and output are for the proleptic Julian calendar,
  242. i.e., no consideration of historical usage of the calendar is
  243. made.
  244. Parameters
  245. ----------
  246. year : int
  247. Year as an integer.
  248. month : int
  249. Month as an integer.
  250. day : int
  251. Day as an integer.
  252. Returns
  253. -------
  254. jd1, jd2: 2-element tuple of floats
  255. When added together, the numbers give the Julian date for the
  256. given Julian calendar date. The first number is always
  257. MJD_0 i.e., 2451545.5. So the second is the MJD.
  258. Examples
  259. --------
  260. >>> jcal2jd(2000, 1, 1)
  261. (2400000.5, 51557.0)
  262. >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678, 2000,
  263. ...: 2012, 2245]
  264. >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12]
  265. >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
  266. >>> x = [jcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
  267. >>> for i in x: print i
  268. (2400000.5, -2395252.0)
  269. (2400000.5, -1451039.0)
  270. (2400000.5, -1062374.0)
  271. (2400000.5, -723765.0)
  272. (2400000.5, -679164.0)
  273. (2400000.5, -678776.0)
  274. (2400000.5, -678370.0)
  275. (2400000.5, -633798.0)
  276. (2400000.5, -65772.0)
  277. (2400000.5, 51871.0)
  278. (2400000.5, 56285.0)
  279. Notes
  280. -----
  281. Unlike `gcal2jd`, negative months and days can result in incorrect
  282. Julian dates.
  283. """
  284. year = int(year)
  285. month = int(month)
  286. day = int(day)
  287. jd = 367 * year
  288. x = ipart((month - 9) / 7.0)
  289. jd -= ipart((7 * (year + 5001 + x)) / 4.0)
  290. jd += ipart((275 * month) / 9.0)
  291. jd += day
  292. jd += 1729777 - 2400000.5 # Return 240000.5 as first part of JD.
  293. jd -= 0.5 # Convert midday to midnight.
  294. return MJD_0, jd
  295. def jd2jcal(jd1, jd2):
  296. """Julian calendar date for the given Julian date.
  297. The input and output are for the proleptic Julian calendar,
  298. i.e., no consideration of historical usage of the calendar is
  299. made.
  300. Parameters
  301. ----------
  302. jd1, jd2: int
  303. Sum of the two numbers is taken as the given Julian date. For
  304. example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
  305. can be the MJD of the date and time. But any combination will
  306. work.
  307. Returns
  308. -------
  309. y, m, d, f : int, int, int, float
  310. Four element tuple containing year, month, day and the
  311. fractional part of the day in the Julian calendar. The first
  312. three are integers, and the last part is a float.
  313. Examples
  314. --------
  315. >>> jd2jcal(*jcal2jd(2000, 1, 1))
  316. (2000, 1, 1, 0.0)
  317. >>> jd2jcal(*jcal2jd(-4000, 10, 11))
  318. (-4000, 10, 11, 0.0)
  319. >>> jcal2jd(2000, 1, 1)
  320. (2400000.5, 51557.0)
  321. >>> jd2jcal(2400000.5, 51557.0)
  322. (2000, 1, 1, 0.0)
  323. >>> jd2jcal(2400000.5, 51557.5)
  324. (2000, 1, 1, 0.5)
  325. >>> jd2jcal(2400000.5, 51557.245)
  326. (2000, 1, 1, 0.24500000000261934)
  327. >>> jd2jcal(2400000.5, 51557.1)
  328. (2000, 1, 1, 0.099999999998544808)
  329. >>> jd2jcal(2400000.5, 51557.75)
  330. (2000, 1, 1, 0.75)
  331. """
  332. from math import modf
  333. jd1_f, jd1_i = modf(jd1)
  334. jd2_f, jd2_i = modf(jd2)
  335. jd_i = jd1_i + jd2_i
  336. f = jd1_f + jd2_f
  337. # Set JD to noon of the current date. Fractional part is the
  338. # fraction from midnight of the current date.
  339. if -0.5 < f < 0.5:
  340. f += 0.5
  341. elif f >= 0.5:
  342. jd_i += 1
  343. f -= 0.5
  344. elif f <= -0.5:
  345. jd_i -= 1
  346. f += 1.5
  347. j = jd_i + 1402.0
  348. k = ipart((j - 1) / 1461.0)
  349. l = j - (1461.0 * k)
  350. n = ipart((l - 1) / 365.0) - ipart(l / 1461.0)
  351. i = l - (365.0 * n) + 30.0
  352. j = ipart((80.0 * i) / 2447.0)
  353. day = i - ipart((2447.0 * j) / 80.0)
  354. i = ipart(j / 11.0)
  355. month = j + 2 - (12.0 * i)
  356. year = (4 * k) + n + i - 4716.0
  357. return int(year), int(month), int(day), f
  358. # Some tests.
  359. def _test_gcal2jd_with_sla_cldj():
  360. """Compare gcal2jd with slalib.sla_cldj."""
  361. import random
  362. try:
  363. from pyslalib import slalib
  364. except ImportError:
  365. print("SLALIB (PySLALIB not available).")
  366. return 1
  367. n = 1000
  368. mday = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
  369. # sla_cldj needs year > -4699 i.e., 4700 BC.
  370. year = [random.randint(-4699, 2200) for i in range(n)]
  371. month = [random.randint(1, 12) for i in range(n)]
  372. day = [random.randint(1, 31) for i in range(n)]
  373. for i in range(n):
  374. x = 0
  375. if is_leap(year[i]) and month[i] == 2:
  376. x = 1
  377. if day[i] > mday[month[i]] + x:
  378. day[i] = mday[month[i]]
  379. jd_jdc = [gcal2jd(y, m, d)[1]
  380. for y, m, d in zip(year, month, day)]
  381. jd_sla = [slalib.sla_cldj(y, m, d)[0]
  382. for y, m, d in zip(year, month, day)]
  383. diff = [abs(i - j) for i, j in zip(jd_sla, jd_jdc)]
  384. assert max(diff) <= 1e-8
  385. assert min(diff) <= 1e-8
  386. def _test_jd2gcal():
  387. """Check jd2gcal as reverse of gcal2jd."""
  388. import random
  389. n = 1000
  390. mday = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
  391. year = [random.randint(-4699, 2200) for i in range(n)]
  392. month = [random.randint(1, 12) for i in range(n)]
  393. day = [random.randint(1, 31) for i in range(n)]
  394. for i in range(n):
  395. x = 0
  396. if is_leap(year[i]) and month[i] == 2:
  397. x = 1
  398. if day[i] > mday[month[i]] + x:
  399. day[i] = mday[month[i]]
  400. jd = [gcal2jd(y, m, d)[1]
  401. for y, m, d in zip(year, month, day)]
  402. x = [jd2gcal(MJD_0, i) for i in jd]
  403. for i in range(n):
  404. assert x[i][0] == year[i]
  405. assert x[i][1] == month[i]
  406. assert x[i][2] == day[i]
  407. assert x[i][3] <= 1e-15
  408. def _test_jd2jcal():
  409. """Check jd2jcal as reverse of jcal2jd."""
  410. import random
  411. n = 1000
  412. year = [random.randint(-4699, 2200) for i in range(n)]
  413. month = [random.randint(1, 12) for i in range(n)]
  414. day = [random.randint(1, 28) for i in range(n)]
  415. jd = [jcal2jd(y, m, d)[1]
  416. for y, m, d in zip(year, month, day)]
  417. x = [jd2gcal(MJD_0, i) for i in jd]
  418. for i in range(n):
  419. assert x[i][0] == year[i]
  420. assert x[i][1] == month[i]
  421. assert x[i][2] == day[i]
  422. assert x[i][3] <= 1e-15