123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440 |
- # -*- coding:utf-8 -*-
- """Functions for converting between Julian dates and calendar dates.
- A function for converting Gregorian calendar dates to Julian dates, and
- another function for converting Julian calendar dates to Julian dates
- are defined. Two functions for the reverse calculations are also
- defined.
- Different regions of the world switched to Gregorian calendar from
- Julian calendar on different dates. Having separate functions for Julian
- and Gregorian calendars allow maximum flexibility in choosing the
- relevant calendar.
- All the above functions are "proleptic". This means that they work for
- dates on which the concerned calendar is not valid. For example,
- Gregorian calendar was not used prior to around October 1582.
- Julian dates are stored in two floating point numbers (double). Julian
- dates, and Modified Julian dates, are large numbers. If only one number
- is used, then the precision of the time stored is limited. Using two
- numbers, time can be split in a manner that will allow maximum
- precision. For example, the first number could be the Julian date for
- the beginning of a day and the second number could be the fractional
- day. Calculations that need the latter part can now work with maximum
- precision.
- A function to test if a given Gregorian calendar year is a leap year is
- defined.
- Zero point of Modified Julian Date (MJD) and the MJD of 2000/1/1
- 12:00:00 are also given.
- This module is based on the TPM C library, by Jeffery W. Percival. The
- idea for splitting Julian date into two floating point numbers was
- inspired by the IAU SOFA C library.
- :author: Prasanth Nair
- :contact: prasanthhn@gmail.com
- :license: BSD (https://opensource.org/licenses/bsd-license.php)
- """
- from __future__ import division
- from __future__ import print_function
- import math
- __version__ = "1.4.1"
- MJD_0 = 2400000.5
- MJD_JD2000 = 51544.5
- def ipart(x):
- """Return integer part of given number."""
- return math.modf(x)[1]
- def is_leap(year):
- """Leap year or not in the Gregorian calendar."""
- x = math.fmod(year, 4)
- y = math.fmod(year, 100)
- z = math.fmod(year, 400)
- # Divisible by 4 and,
- # either not divisible by 100 or divisible by 400.
- return not x and (y or not z)
- def gcal2jd(year, month, day):
- """Gregorian calendar date to Julian date.
- The input and output are for the proleptic Gregorian calendar,
- i.e., no consideration of historical usage of the calendar is
- made.
- Parameters
- ----------
- year : int
- Year as an integer.
- month : int
- Month as an integer.
- day : int
- Day as an integer.
- Returns
- -------
- jd1, jd2: 2-element tuple of floats
- When added together, the numbers give the Julian date for the
- given Gregorian calendar date. The first number is always
- MJD_0 i.e., 2451545.5. So the second is the MJD.
- Examples
- --------
- >>> gcal2jd(2000,1,1)
- (2400000.5, 51544.0)
- >>> 2400000.5 + 51544.0 + 0.5
- 2451545.0
- >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678.0, 2000,
- ....: 2012, 2245]
- >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
- >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
- >>> x = [gcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
- >>> for i in x: print i
- (2400000.5, -2395215.0)
- (2400000.5, -1451021.0)
- (2400000.5, -1062364.0)
- (2400000.5, -723762.0)
- (2400000.5, -679162.0)
- (2400000.5, -678774.0)
- (2400000.5, -678368.0)
- (2400000.5, -633797.0)
- (2400000.5, -65812.0)
- (2400000.5, 51827.0)
- (2400000.5, 56242.0)
- (2400000.5, 141393.0)
- Negative months and days are valid. For example, 2000/-2/-4 =>
- 1999/+12-2/-4 => 1999/10/-4 => 1999/9/30-4 => 1999/9/26.
- >>> gcal2jd(2000, -2, -4)
- (2400000.5, 51447.0)
- >>> gcal2jd(1999, 9, 26)
- (2400000.5, 51447.0)
- >>> gcal2jd(2000, 2, -1)
- (2400000.5, 51573.0)
- >>> gcal2jd(2000, 1, 30)
- (2400000.5, 51573.0)
- >>> gcal2jd(2000, 3, -1)
- (2400000.5, 51602.0)
- >>> gcal2jd(2000, 2, 28)
- (2400000.5, 51602.0)
- Month 0 becomes previous month.
- >>> gcal2jd(2000, 0, 1)
- (2400000.5, 51513.0)
- >>> gcal2jd(1999, 12, 1)
- (2400000.5, 51513.0)
- Day number 0 becomes last day of previous month.
- >>> gcal2jd(2000, 3, 0)
- (2400000.5, 51603.0)
- >>> gcal2jd(2000, 2, 29)
- (2400000.5, 51603.0)
- If `day` is greater than the number of days in `month`, then it
- gets carried over to the next month.
- >>> gcal2jd(2000,2,30)
- (2400000.5, 51604.0)
- >>> gcal2jd(2000,3,1)
- (2400000.5, 51604.0)
- >>> gcal2jd(2001,2,30)
- (2400000.5, 51970.0)
- >>> gcal2jd(2001,3,2)
- (2400000.5, 51970.0)
- Notes
- -----
- The returned Julian date is for mid-night of the given date. To
- find the Julian date for any time of the day, simply add time as a
- fraction of a day. For example Julian date for mid-day can be
- obtained by adding 0.5 to either the first part or the second
- part. The latter is preferable, since it will give the MJD for the
- date and time.
- BC dates should be given as -(BC - 1) where BC is the year. For
- example 1 BC == 0, 2 BC == -1, and so on.
- Negative numbers can be used for `month` and `day`. For example
- 2000, -1, 1 is the same as 1999, 11, 1.
- The Julian dates are proleptic Julian dates, i.e., values are
- returned without considering if Gregorian dates are valid for the
- given date.
- The input values are truncated to integers.
- """
- year = int(year)
- month = int(month)
- day = int(day)
- a = ipart((month - 14) / 12.0)
- jd = ipart((1461 * (year + 4800 + a)) / 4.0)
- jd += ipart((367 * (month - 2 - 12 * a)) / 12.0)
- x = ipart((year + 4900 + a) / 100.0)
- jd -= ipart((3 * x) / 4.0)
- jd += day - 2432075.5 # was 32075; add 2400000.5
- jd -= 0.5 # 0 hours; above JD is for midday, switch to midnight.
- return MJD_0, jd
- def jd2gcal(jd1, jd2):
- """Julian date to Gregorian calendar date and time of day.
- The input and output are for the proleptic Gregorian calendar,
- i.e., no consideration of historical usage of the calendar is
- made.
- Parameters
- ----------
- jd1, jd2: int
- Sum of the two numbers is taken as the given Julian date. For
- example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
- can be the MJD of the date and time. But any combination will
- work.
- Returns
- -------
- y, m, d, f : int, int, int, float
- Four element tuple containing year, month, day and the
- fractional part of the day in the Gregorian calendar. The first
- three are integers, and the last part is a float.
- Examples
- --------
- >>> jd2gcal(*gcal2jd(2000,1,1))
- (2000, 1, 1, 0.0)
- >>> jd2gcal(*gcal2jd(1950,1,1))
- (1950, 1, 1, 0.0)
- Out of range months and days are carried over to the next/previous
- year or next/previous month. See gcal2jd for more examples.
- >>> jd2gcal(*gcal2jd(1999,10,12))
- (1999, 10, 12, 0.0)
- >>> jd2gcal(*gcal2jd(2000,2,30))
- (2000, 3, 1, 0.0)
- >>> jd2gcal(*gcal2jd(-1999,10,12))
- (-1999, 10, 12, 0.0)
- >>> jd2gcal(*gcal2jd(2000, -2, -4))
- (1999, 9, 26, 0.0)
- >>> gcal2jd(2000,1,1)
- (2400000.5, 51544.0)
- >>> jd2gcal(2400000.5, 51544.0)
- (2000, 1, 1, 0.0)
- >>> jd2gcal(2400000.5, 51544.5)
- (2000, 1, 1, 0.5)
- >>> jd2gcal(2400000.5, 51544.245)
- (2000, 1, 1, 0.24500000000261934)
- >>> jd2gcal(2400000.5, 51544.1)
- (2000, 1, 1, 0.099999999998544808)
- >>> jd2gcal(2400000.5, 51544.75)
- (2000, 1, 1, 0.75)
- Notes
- -----
- The last element of the tuple is the same as
- (hh + mm / 60.0 + ss / 3600.0) / 24.0
- where hh, mm, and ss are the hour, minute and second of the day.
- See Also
- --------
- gcal2jd
- """
- from math import modf
- jd1_f, jd1_i = modf(jd1)
- jd2_f, jd2_i = modf(jd2)
- jd_i = jd1_i + jd2_i
- f = jd1_f + jd2_f
- # Set JD to noon of the current date. Fractional part is the
- # fraction from midnight of the current date.
- if -0.5 < f < 0.5:
- f += 0.5
- elif f >= 0.5:
- jd_i += 1
- f -= 0.5
- elif f <= -0.5:
- jd_i -= 1
- f += 1.5
- l = jd_i + 68569
- n = ipart((4 * l) / 146097.0)
- l -= ipart(((146097 * n) + 3) / 4.0)
- i = ipart((4000 * (l + 1)) / 1461001)
- l -= ipart((1461 * i) / 4.0) - 31
- j = ipart((80 * l) / 2447.0)
- day = l - ipart((2447 * j) / 80.0)
- l = ipart(j / 11.0)
- month = j + 2 - (12 * l)
- year = 100 * (n - 49) + i + l
- return int(year), int(month), int(day), f
- def jcal2jd(year, month, day):
- """Julian calendar date to Julian date.
- The input and output are for the proleptic Julian calendar,
- i.e., no consideration of historical usage of the calendar is
- made.
- Parameters
- ----------
- year : int
- Year as an integer.
- month : int
- Month as an integer.
- day : int
- Day as an integer.
- Returns
- -------
- jd1, jd2: 2-element tuple of floats
- When added together, the numbers give the Julian date for the
- given Julian calendar date. The first number is always
- MJD_0 i.e., 2451545.5. So the second is the MJD.
- Examples
- --------
- >>> jcal2jd(2000, 1, 1)
- (2400000.5, 51557.0)
- >>> year = [-4699, -2114, -1050, -123, -1, 0, 1, 123, 1678, 2000,
- ...: 2012, 2245]
- >>> month = [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12]
- >>> day = [1, 12, 23, 14, 25, 16, 27, 8, 9, 10, 11, 31]
- >>> x = [jcal2jd(y, m, d) for y, m, d in zip(year, month, day)]
- >>> for i in x: print i
- (2400000.5, -2395252.0)
- (2400000.5, -1451039.0)
- (2400000.5, -1062374.0)
- (2400000.5, -723765.0)
- (2400000.5, -679164.0)
- (2400000.5, -678776.0)
- (2400000.5, -678370.0)
- (2400000.5, -633798.0)
- (2400000.5, -65772.0)
- (2400000.5, 51871.0)
- (2400000.5, 56285.0)
- Notes
- -----
- Unlike `gcal2jd`, negative months and days can result in incorrect
- Julian dates.
- """
- year = int(year)
- month = int(month)
- day = int(day)
- jd = 367 * year
- x = ipart((month - 9) / 7.0)
- jd -= ipart((7 * (year + 5001 + x)) / 4.0)
- jd += ipart((275 * month) / 9.0)
- jd += day
- jd += 1729777 - 2400000.5 # Return 240000.5 as first part of JD.
- jd -= 0.5 # Convert midday to midnight.
- return MJD_0, jd
- def jd2jcal(jd1, jd2):
- """Julian calendar date for the given Julian date.
- The input and output are for the proleptic Julian calendar,
- i.e., no consideration of historical usage of the calendar is
- made.
- Parameters
- ----------
- jd1, jd2: int
- Sum of the two numbers is taken as the given Julian date. For
- example `jd1` can be the zero point of MJD (MJD_0) and `jd2`
- can be the MJD of the date and time. But any combination will
- work.
- Returns
- -------
- y, m, d, f : int, int, int, float
- Four element tuple containing year, month, day and the
- fractional part of the day in the Julian calendar. The first
- three are integers, and the last part is a float.
- Examples
- --------
- >>> jd2jcal(*jcal2jd(2000, 1, 1))
- (2000, 1, 1, 0.0)
- >>> jd2jcal(*jcal2jd(-4000, 10, 11))
- (-4000, 10, 11, 0.0)
- >>> jcal2jd(2000, 1, 1)
- (2400000.5, 51557.0)
- >>> jd2jcal(2400000.5, 51557.0)
- (2000, 1, 1, 0.0)
- >>> jd2jcal(2400000.5, 51557.5)
- (2000, 1, 1, 0.5)
- >>> jd2jcal(2400000.5, 51557.245)
- (2000, 1, 1, 0.24500000000261934)
- >>> jd2jcal(2400000.5, 51557.1)
- (2000, 1, 1, 0.099999999998544808)
- >>> jd2jcal(2400000.5, 51557.75)
- (2000, 1, 1, 0.75)
- """
- from math import modf
- jd1_f, jd1_i = modf(jd1)
- jd2_f, jd2_i = modf(jd2)
- jd_i = jd1_i + jd2_i
- f = jd1_f + jd2_f
- # Set JD to noon of the current date. Fractional part is the
- # fraction from midnight of the current date.
- if -0.5 < f < 0.5:
- f += 0.5
- elif f >= 0.5:
- jd_i += 1
- f -= 0.5
- elif f <= -0.5:
- jd_i -= 1
- f += 1.5
- j = jd_i + 1402.0
- k = ipart((j - 1) / 1461.0)
- l = j - (1461.0 * k)
- n = ipart((l - 1) / 365.0) - ipart(l / 1461.0)
- i = l - (365.0 * n) + 30.0
- j = ipart((80.0 * i) / 2447.0)
- day = i - ipart((2447.0 * j) / 80.0)
- i = ipart(j / 11.0)
- month = j + 2 - (12.0 * i)
- year = (4 * k) + n + i - 4716.0
- return int(year), int(month), int(day), f
|