jdcal.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440
  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 (https://opensource.org/licenses/bsd-license.php)
  32. """
  33. from __future__ import division
  34. from __future__ import print_function
  35. import math
  36. __version__ = "1.4.1"
  37. MJD_0 = 2400000.5
  38. MJD_JD2000 = 51544.5
  39. def ipart(x):
  40. """Return integer part of given number."""
  41. return math.modf(x)[1]
  42. def is_leap(year):
  43. """Leap year or not in the Gregorian calendar."""
  44. x = math.fmod(year, 4)
  45. y = math.fmod(year, 100)
  46. z = math.fmod(year, 400)
  47. # Divisible by 4 and,
  48. # either not divisible by 100 or divisible by 400.
  49. return not x and (y or not z)
  50. def gcal2jd(year, month, day):
  51. """Gregorian calendar date to Julian date.
  52. The input and output are for the proleptic Gregorian calendar,
  53. i.e., no consideration of historical usage of the calendar is
  54. made.
  55. Parameters
  56. ----------
  57. year : int
  58. Year as an integer.
  59. month : int
  60. Month as an integer.
  61. day : int
  62. Day as an integer.
  63. Returns
  64. -------
  65. jd1, jd2: 2-element tuple of floats
  66. When added together, the numbers give the Julian date for the
  67. given Gregorian calendar date. The first number is always
  68. MJD_0 i.e., 2451545.5. So the second is the MJD.
  69. Examples
  70. --------
  71. >>> gcal2jd(2000,1,1)
  72. (2400000.5, 51544.0)
  73. >>> 2400000.5 + 51544.0 + 0.5
  74. 2451545.0
  75. >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678.0, 2000,
  76. ....: 2012, 2245]
  77. >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
  78. >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
  79. >>> x = [gcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
  80. >>> for i in x: print i
  81. (2400000.5, -2395215.0)
  82. (2400000.5, -1451021.0)
  83. (2400000.5, -1062364.0)
  84. (2400000.5, -723762.0)
  85. (2400000.5, -679162.0)
  86. (2400000.5, -678774.0)
  87. (2400000.5, -678368.0)
  88. (2400000.5, -633797.0)
  89. (2400000.5, -65812.0)
  90. (2400000.5, 51827.0)
  91. (2400000.5, 56242.0)
  92. (2400000.5, 141393.0)
  93. Negative months and days are valid. For example, 2000/-2/-4 =>
  94. 1999/+12-2/-4 => 1999/10/-4 => 1999/9/30-4 => 1999/9/26.
  95. >>> gcal2jd(2000, -2, -4)
  96. (2400000.5, 51447.0)
  97. >>> gcal2jd(1999, 9, 26)
  98. (2400000.5, 51447.0)
  99. >>> gcal2jd(2000, 2, -1)
  100. (2400000.5, 51573.0)
  101. >>> gcal2jd(2000, 1, 30)
  102. (2400000.5, 51573.0)
  103. >>> gcal2jd(2000, 3, -1)
  104. (2400000.5, 51602.0)
  105. >>> gcal2jd(2000, 2, 28)
  106. (2400000.5, 51602.0)
  107. Month 0 becomes previous month.
  108. >>> gcal2jd(2000, 0, 1)
  109. (2400000.5, 51513.0)
  110. >>> gcal2jd(1999, 12, 1)
  111. (2400000.5, 51513.0)
  112. Day number 0 becomes last day of previous month.
  113. >>> gcal2jd(2000, 3, 0)
  114. (2400000.5, 51603.0)
  115. >>> gcal2jd(2000, 2, 29)
  116. (2400000.5, 51603.0)
  117. If `day` is greater than the number of days in `month`, then it
  118. gets carried over to the next month.
  119. >>> gcal2jd(2000,2,30)
  120. (2400000.5, 51604.0)
  121. >>> gcal2jd(2000,3,1)
  122. (2400000.5, 51604.0)
  123. >>> gcal2jd(2001,2,30)
  124. (2400000.5, 51970.0)
  125. >>> gcal2jd(2001,3,2)
  126. (2400000.5, 51970.0)
  127. Notes
  128. -----
  129. The returned Julian date is for mid-night of the given date. To
  130. find the Julian date for any time of the day, simply add time as a
  131. fraction of a day. For example Julian date for mid-day can be
  132. obtained by adding 0.5 to either the first part or the second
  133. part. The latter is preferable, since it will give the MJD for the
  134. date and time.
  135. BC dates should be given as -(BC - 1) where BC is the year. For
  136. example 1 BC == 0, 2 BC == -1, and so on.
  137. Negative numbers can be used for `month` and `day`. For example
  138. 2000, -1, 1 is the same as 1999, 11, 1.
  139. The Julian dates are proleptic Julian dates, i.e., values are
  140. returned without considering if Gregorian dates are valid for the
  141. given date.
  142. The input values are truncated to integers.
  143. """
  144. year = int(year)
  145. month = int(month)
  146. day = int(day)
  147. a = ipart((month - 14) / 12.0)
  148. jd = ipart((1461 * (year + 4800 + a)) / 4.0)
  149. jd += ipart((367 * (month - 2 - 12 * a)) / 12.0)
  150. x = ipart((year + 4900 + a) / 100.0)
  151. jd -= ipart((3 * x) / 4.0)
  152. jd += day - 2432075.5 # was 32075; add 2400000.5
  153. jd -= 0.5 # 0 hours; above JD is for midday, switch to midnight.
  154. return MJD_0, jd
  155. def jd2gcal(jd1, jd2):
  156. """Julian date to Gregorian calendar date and time of day.
  157. The input and output are for the proleptic Gregorian calendar,
  158. i.e., no consideration of historical usage of the calendar is
  159. made.
  160. Parameters
  161. ----------
  162. jd1, jd2: int
  163. Sum of the two numbers is taken as the given Julian date. For
  164. example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
  165. can be the MJD of the date and time. But any combination will
  166. work.
  167. Returns
  168. -------
  169. y, m, d, f : int, int, int, float
  170. Four element tuple containing year, month, day and the
  171. fractional part of the day in the Gregorian calendar. The first
  172. three are integers, and the last part is a float.
  173. Examples
  174. --------
  175. >>> jd2gcal(*gcal2jd(2000,1,1))
  176. (2000, 1, 1, 0.0)
  177. >>> jd2gcal(*gcal2jd(1950,1,1))
  178. (1950, 1, 1, 0.0)
  179. Out of range months and days are carried over to the next/previous
  180. year or next/previous month. See gcal2jd for more examples.
  181. >>> jd2gcal(*gcal2jd(1999,10,12))
  182. (1999, 10, 12, 0.0)
  183. >>> jd2gcal(*gcal2jd(2000,2,30))
  184. (2000, 3, 1, 0.0)
  185. >>> jd2gcal(*gcal2jd(-1999,10,12))
  186. (-1999, 10, 12, 0.0)
  187. >>> jd2gcal(*gcal2jd(2000, -2, -4))
  188. (1999, 9, 26, 0.0)
  189. >>> gcal2jd(2000,1,1)
  190. (2400000.5, 51544.0)
  191. >>> jd2gcal(2400000.5, 51544.0)
  192. (2000, 1, 1, 0.0)
  193. >>> jd2gcal(2400000.5, 51544.5)
  194. (2000, 1, 1, 0.5)
  195. >>> jd2gcal(2400000.5, 51544.245)
  196. (2000, 1, 1, 0.24500000000261934)
  197. >>> jd2gcal(2400000.5, 51544.1)
  198. (2000, 1, 1, 0.099999999998544808)
  199. >>> jd2gcal(2400000.5, 51544.75)
  200. (2000, 1, 1, 0.75)
  201. Notes
  202. -----
  203. The last element of the tuple is the same as
  204. (hh + mm / 60.0 + ss / 3600.0) / 24.0
  205. where hh, mm, and ss are the hour, minute and second of the day.
  206. See Also
  207. --------
  208. gcal2jd
  209. """
  210. from math import modf
  211. jd1_f, jd1_i = modf(jd1)
  212. jd2_f, jd2_i = modf(jd2)
  213. jd_i = jd1_i + jd2_i
  214. f = jd1_f + jd2_f
  215. # Set JD to noon of the current date. Fractional part is the
  216. # fraction from midnight of the current date.
  217. if -0.5 < f < 0.5:
  218. f += 0.5
  219. elif f >= 0.5:
  220. jd_i += 1
  221. f -= 0.5
  222. elif f <= -0.5:
  223. jd_i -= 1
  224. f += 1.5
  225. l = jd_i + 68569
  226. n = ipart((4 * l) / 146097.0)
  227. l -= ipart(((146097 * n) + 3) / 4.0)
  228. i = ipart((4000 * (l + 1)) / 1461001)
  229. l -= ipart((1461 * i) / 4.0) - 31
  230. j = ipart((80 * l) / 2447.0)
  231. day = l - ipart((2447 * j) / 80.0)
  232. l = ipart(j / 11.0)
  233. month = j + 2 - (12 * l)
  234. year = 100 * (n - 49) + i + l
  235. return int(year), int(month), int(day), f
  236. def jcal2jd(year, month, day):
  237. """Julian calendar date to Julian date.
  238. The input and output are for the proleptic Julian calendar,
  239. i.e., no consideration of historical usage of the calendar is
  240. made.
  241. Parameters
  242. ----------
  243. year : int
  244. Year as an integer.
  245. month : int
  246. Month as an integer.
  247. day : int
  248. Day as an integer.
  249. Returns
  250. -------
  251. jd1, jd2: 2-element tuple of floats
  252. When added together, the numbers give the Julian date for the
  253. given Julian calendar date. The first number is always
  254. MJD_0 i.e., 2451545.5. So the second is the MJD.
  255. Examples
  256. --------
  257. >>> jcal2jd(2000, 1, 1)
  258. (2400000.5, 51557.0)
  259. >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678, 2000,
  260. ...: 2012, 2245]
  261. >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12]
  262. >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
  263. >>> x = [jcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
  264. >>> for i in x: print i
  265. (2400000.5, -2395252.0)
  266. (2400000.5, -1451039.0)
  267. (2400000.5, -1062374.0)
  268. (2400000.5, -723765.0)
  269. (2400000.5, -679164.0)
  270. (2400000.5, -678776.0)
  271. (2400000.5, -678370.0)
  272. (2400000.5, -633798.0)
  273. (2400000.5, -65772.0)
  274. (2400000.5, 51871.0)
  275. (2400000.5, 56285.0)
  276. Notes
  277. -----
  278. Unlike `gcal2jd`, negative months and days can result in incorrect
  279. Julian dates.
  280. """
  281. year = int(year)
  282. month = int(month)
  283. day = int(day)
  284. jd = 367 * year
  285. x = ipart((month - 9) / 7.0)
  286. jd -= ipart((7 * (year + 5001 + x)) / 4.0)
  287. jd += ipart((275 * month) / 9.0)
  288. jd += day
  289. jd += 1729777 - 2400000.5 # Return 240000.5 as first part of JD.
  290. jd -= 0.5 # Convert midday to midnight.
  291. return MJD_0, jd
  292. def jd2jcal(jd1, jd2):
  293. """Julian calendar date for the given Julian date.
  294. The input and output are for the proleptic Julian calendar,
  295. i.e., no consideration of historical usage of the calendar is
  296. made.
  297. Parameters
  298. ----------
  299. jd1, jd2: int
  300. Sum of the two numbers is taken as the given Julian date. For
  301. example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
  302. can be the MJD of the date and time. But any combination will
  303. work.
  304. Returns
  305. -------
  306. y, m, d, f : int, int, int, float
  307. Four element tuple containing year, month, day and the
  308. fractional part of the day in the Julian calendar. The first
  309. three are integers, and the last part is a float.
  310. Examples
  311. --------
  312. >>> jd2jcal(*jcal2jd(2000, 1, 1))
  313. (2000, 1, 1, 0.0)
  314. >>> jd2jcal(*jcal2jd(-4000, 10, 11))
  315. (-4000, 10, 11, 0.0)
  316. >>> jcal2jd(2000, 1, 1)
  317. (2400000.5, 51557.0)
  318. >>> jd2jcal(2400000.5, 51557.0)
  319. (2000, 1, 1, 0.0)
  320. >>> jd2jcal(2400000.5, 51557.5)
  321. (2000, 1, 1, 0.5)
  322. >>> jd2jcal(2400000.5, 51557.245)
  323. (2000, 1, 1, 0.24500000000261934)
  324. >>> jd2jcal(2400000.5, 51557.1)
  325. (2000, 1, 1, 0.099999999998544808)
  326. >>> jd2jcal(2400000.5, 51557.75)
  327. (2000, 1, 1, 0.75)
  328. """
  329. from math import modf
  330. jd1_f, jd1_i = modf(jd1)
  331. jd2_f, jd2_i = modf(jd2)
  332. jd_i = jd1_i + jd2_i
  333. f = jd1_f + jd2_f
  334. # Set JD to noon of the current date. Fractional part is the
  335. # fraction from midnight of the current date.
  336. if -0.5 < f < 0.5:
  337. f += 0.5
  338. elif f >= 0.5:
  339. jd_i += 1
  340. f -= 0.5
  341. elif f <= -0.5:
  342. jd_i -= 1
  343. f += 1.5
  344. j = jd_i + 1402.0
  345. k = ipart((j - 1) / 1461.0)
  346. l = j - (1461.0 * k)
  347. n = ipart((l - 1) / 365.0) - ipart(l / 1461.0)
  348. i = l - (365.0 * n) + 30.0
  349. j = ipart((80.0 * i) / 2447.0)
  350. day = i - ipart((2447.0 * j) / 80.0)
  351. i = ipart(j / 11.0)
  352. month = j + 2 - (12.0 * i)
  353. year = (4 * k) + n + i - 4716.0
  354. return int(year), int(month), int(day), f