When I was investigating the calculation method of Julian day, I was wondering which formula was correct depending on the source, so I organized it for the time being. I also made a comparison, and although I found that there were differences in the results of each calculation method, it is unclear which one is correct or whether it is possible to say that one method is correct in the first place. As it is.

Julian day calculation method

[Wikipedia / Julian Day Number (JDN)](https://ja.wikipedia.org/wiki/%E3%83%A6%E3%83%AA%E3%82%A6%E3%82% According to B9% E9% 80% 9A% E6% 97% A5 # Julian_Day_Number_ (JDN))

The conversion formula was devised by Fliegel and Van Flandern [13], Hatcher [14], Meeus [15]. However, conversion formulas that have been arranged are often used [16].

It seems that the formula has several shapes. The formulas in wikipedia and the formulas in other documents are listed below.

The form of each formula is quite different, but is it really equivalent?

wikipedia method

[Wikipedia / Sidereal time # Sidereal time calculation method](https://ja.wikipedia.org/wiki/%E6%81%92%E6%98%9F%E6%99%82#%E6%81%92 % E6% 98% 9F% E6% 99% 82% E3% 81% AE% E8% A8% 88% E7% AE% 97% E6% B3% 95)

\begin{aligned}
JD = & \left\lfloor 365.25 Y \right\rfloor + \left \lfloor \frac{Y}{400} \right\rfloor - \left\lfloor \frac{Y}{100} \right\rfloor + \left\lfloor 30.59 \left( M-2 \right) \right\rfloor + D + 1721088.5 \\
& + \frac{h}{24} + \frac{m}{1440} + \frac{s}{86400} \\
\\
& \lfloor x \rfloor is the floor function
\end{aligned}

Let the year in the current Gregorian calendar of UT be Y, the month be M, the day be D, the hour be h, the minute be m, and the second be s. However, January and February are substituted as March and 14 of the previous year (Y value is decremented by -1) (Example: In the case of February 5, 2013, Y = 2012, M = 14, D = Five).

Fliegel et al.'S method

Fliegel, H. F. and Van Flandern, T. C., "A Machine Algorithm for Processing Calendar Dates," Communications of the ACM 11, p. 657, 1968.

\begin{aligned}
JD(I,J,K) = & K - 32075 + 1461 * (I + 4800 + (J - 14) / 12) / 4 \\
& + 367 * (J - 2 - (J - 14) / 12 * 12) / 12 \\
& - 3 * ((I + 4900 + (J - 14) / 12 ) / 100 ) / 4
\end{aligned}

calendar date (I = year; J = month, a number from 1 to 12; K = day of month) to a Julian Date (JD)

COMPUTATION AND MEASUREMENT PROBLEM 11.

In FORTRAN integer arithmetic, multiplication and division are performed left to right in the order of occurrence, and the absolute value of each result is truncated to the next lower integer value after each operation, so that both 2/12 and -2/12 become 0.

Multiplication and division are performed from left to right, each time the absolute value of the calculation is rounded down to the next smaller integer value.

wikipedia / Julian day#Converting Gregorian calendar date to Julian Day Number

The algorithm[61] is valid for all (possibly proleptic) Gregorian calendar dates after November 23, −4713. Divisions are integer divisions, fractional parts are ignored.

If you change the expression so that it is easier to compare with others,

\begin{aligned}
JD = & D - 32075 + \left\lfloor \frac{ 1461 \left( Y + 4800 + \left\lfloor \frac{M - 14}{12} \right\rfloor \right)}{4} \right\rfloor \\
& + \left\lfloor \frac{ 367 \left( M - 2 - 12 \left\lfloor \frac{M - 14}{12} \right\rfloor \right) }{12} \right\rfloor \\
& - \left\lfloor \frac{ 3 \left\lfloor \frac{ Y + 4900 + \left\lfloor \frac{M - 14}{12} \right\rfloor }{100} \right\rfloor }{ 4 } \right\rfloor
\end{aligned}

Hatcher's method

Hatcher, D. A., Simple formulae for Julian day numbers and calendar dates, Quarterly Journal of the Royal Astronomical Society, v. 25, p. 53-55, 1984

\begin{aligned}
Y' & = Y + y - ((n + m - 1 - M) / n) \text{INT} \\
M' & = (M - m) \text{MOD } n \\
D' & = D - 1 \\
J & = ((pY' + q) / r) \text{INT} + ((sM' + t) / u) \text{INT} + D' -j
\end{aligned}

Y, M, D and Y', M', D' are the year, month and day of month

(x) INT is the integral part of x, and x is assumed to be positive; (x) MOD y is the positive remainder on divideing x by y.

(x) INT is the integer part of x, and x is assumed to be positive. (x) MOD is the positive remainder of x divided by y.

For conversion from the Gregorian calendar

y	j	m	n	r	p	q	v	u	s	t	w
4716	1401+g	3	12	4	1461	-	3	5	153	2	2

g = (((Y' + 184) / 100) \text{INT} \times 3/4) \text{INT} - 38

Therefore, if you substitute it in the above formula,

\begin{aligned}
Y' = & Y + 4716 - \left( \frac{14 - M}{12} \right) \text{INT} \\
M' = & (M - 3) \text{MOD } 12 \\
D' = & D - 1 \\
J = & \left( \frac{1461 Y'}{4} \right) \text{INT} + \left( \frac{153M' + 2}{5} \right) \text{INT} + D' \\
& - \left( 1401 + \left( \frac{Y' + 184}{100} \right) \text{INT} \times \frac{3}{4} \right) \text{INT} - 38 \\
= & \left( 365.25Y' \right) \text{INT} + \left( 30.6M' + 0.4 \right) \text{INT} + D' \\
& - \left( 1401 + \left( \frac{Y' + 184}{100} \right) \text{INT} \times \frac{3}{4} \right) \text{INT} - 38
\end{aligned}

Meeus method

Meeus, J., Astronomical Algorithms, 1998

\begin{aligned}
A & = \text{INT} \left( \frac{Y}{100} \right) \\
B & = 2 - A + \text{INT} \left( \frac{A}{4} \right) \\
JD & = \text{INT} (365.25 (Y + 4716)) + \text{INT} (30.6001(M + 1)) + D + B - 1524.5
\end{aligned}

Let Y be the year, M the month number (1 for January, 2 for February, etc., to 12 for December), and D the day of month (with decimals, if any) of the given calendar date.

INT(x) the greatest integer less than or equal to x.

for instance, INT(-7.83) = -8

If you change the expression so that it is easier to compare with others,

\begin{aligned}
JD & = \left\lfloor 365.25 (Y + 4716) \right\rfloor + \left\lfloor 30.6001(M + 1) \right\rfloor + D + 2 - \left\lfloor \frac{Y}{100} \right\rfloor + \left\lfloor \frac{ \left\lfloor \frac{Y}{100} \right\rfloor }{4} \right\rfloor - 1524.5 \\
& = \left\lfloor 365.25 Y \right\rfloor + \left\lfloor 30.6001(M + 1) \right\rfloor + D - \left\lfloor \frac{Y}{100} \right\rfloor + \left\lfloor \frac{Y}{400} \right\rfloor + ‭1,720,996.5‬ \\
& = 365 Y + \left\lfloor \frac{Y}{4} \right\rfloor - \left\lfloor \frac{Y}{100} \right\rfloor + \left\lfloor \frac{Y}{400} \right\rfloor + \left\lfloor 30.6001(M + 1) \right\rfloor + D + ‭1,720,996.5‬
\end{aligned}

Vallado et al.'S method

CelesTrak / Revisiting Spacetrack Report #3 - AIAA 2006-6753-Sourcecode(C++) Vallado, David A., Paul Crawford, Richard Hujsak, and T.S. Kelso, "Revisiting Spacetrack Report #3," presented at the AIAA/AAS Astrodynamics Specialist Conference, Keystone, CO, 2006 August 21–24.

\begin{aligned}
jd = & 367.0 * year - \\
& floor((7 * (year + floor((mon + 9) / 12.0))) * 0.25) + \\
& floor(275 * mon / 9.0) + day + 1721013.5 + \\
& ((sec / 60.0 + minute) / 60.0 + hr) / 24.0 \\
\end{aligned}

If you change the expression so that it is easier to compare with others,

\begin{aligned}
JD = & 367 Y - \left\lfloor \frac{7}{4} \left(Y + \left\lfloor \frac{M + 9}{12} \right\rfloor \right) \right\rfloor + \left\lfloor \frac{275 M}{9} \right\rfloor + D + 1721013.5 \\
& + \frac{h}{24} + \frac{m}{1440} + \frac{s}{86400}
\end{aligned}

Howard D. Curtis's method

ScienceDirect / Julian Day Number Howard D. Curtis, in Orbital Mechanics for Engineering Students (Fourth Edition), 2020

\begin{aligned}
J_0 = & 367 y - \text{INT} \frac{7 y + \text{INT} \frac{m + 9}{12}}{4} + \text{INT} \frac{275 m}{9} + d + 1,721,013.5 \\
& 1901 \le y \le 2099 \\
& 1 \le m \le 12 \\
& 1 \le d \le 31
\end{aligned}

INT(x) means retaining only the integer portion of x, without rounding (or, in other words, round toward zero). For example, INT(− 3.9) = −3 and INT(3.9) = 3.

ScienceDirect / Julian Day Number Howard D. Curtis, in Orbital Mechanics for Engineering Students (Third Edition), 2014

\begin{aligned}
J_0 = & 367 y - \text{INT} \left\{ \frac{7 \left[ y + \text{INT} \left( \frac{m + 9}{12} \right) \right] }{4} \right\} + \text{INT} \left( \frac{275 m}{9} \right) + d + 1,721,013.5 \\
& 1901 \le y \le 2099 \\
& 1 \le m \le 12 \\
& 1 \le d \le 31
\end{aligned}

INT (x) means to retain only the integer portion of x, without rounding (or, in other words, round toward zero), that is, INT (−3.9) = −3 and INT (3.9) = 3

In the above Fourth Edition and Third Edition, the range of the numerator 7 in the second term is different.

boost library method

boost/date_time/gregorian_calendar.ipp

  //! Convert a ymd_type into a day number
  /*! The day number is an absolute number of days since the start of count
   */
  template<typename ymd_type_, typename date_int_type_>
  BOOST_DATE_TIME_INLINE
  date_int_type_
  gregorian_calendar_base<ymd_type_,date_int_type_>::day_number(const ymd_type& ymd)
  {
    unsigned short a = static_cast<unsigned short>((14-ymd.month)/12);
    unsigned short y = static_cast<unsigned short>(ymd.year + 4800 - a);
    unsigned short m = static_cast<unsigned short>(ymd.month + 12*a - 3);
    unsigned long  d = ymd.day + ((153*m + 2)/5) + 365*y + (y/4) - (y/100) + (y/400) - 32045;
    return static_cast<date_int_type>(d);
  }

  //! Convert a year-month-day into the julian day number
  /*! Since this implementation uses julian day internally, this is the same as the day_number.
   */
  template<typename ymd_type_, typename date_int_type_>
  BOOST_DATE_TIME_INLINE
  date_int_type_
  gregorian_calendar_base<ymd_type_,date_int_type_>::julian_day_number(const ymd_type& ymd)
  {
    return day_number(ymd);
  }

If you change the expression so that it is easier to compare with others,

\begin{aligned}
a & = \text{INT} \left( \frac{14 - M}{12} \right) \\
y & = Y + 4800 - a \\
m & = M + 12 a - 3 \\
JD & = \text{INT} \left( D + \frac{153m + 2}{5} + 365 y + \frac{y}{4} - \frac{y}{100} + \frac{y}{400} - 32045 \right)
\end{aligned}

PHP method

PHP> Manual> Function Reference> Date and Time Related> Calendar> Calendar Functions gregoriantojd

gregoriantojd — Convert Gregorius days to Julius integration days

Although it is a PHP function, it seems that the entity is implemented in C as follows.

php-src/ext/calendar/calendar.stub.php

`php:calendar.stub.ph`


function gregoriantojd(int $month, int $day, int $year): int {}

php-src/ext/calendar/calendar.c

`calendar.c`


/* {{{ proto int gregoriantojd(int month, int day, int year)
   Converts a gregorian calendar date to julian day count */
PHP_FUNCTION(gregoriantojd)
{
	zend_long year, month, day;

	if (zend_parse_parameters(ZEND_NUM_ARGS(), "lll", &month, &day, &year) == FAILURE) {
		RETURN_THROWS();
	}

	RETURN_LONG(GregorianToSdn(year, month, day));
}
/* }}} */

php-src/ext/calendar/gregor.c

`gregor.c`


zend_long GregorianToSdn(
						   int inputYear,
						   int inputMonth,
						   int inputDay)
{
	zend_long year;
	int month;

	/* check for invalid dates */
	if (inputYear == 0 || inputYear < -4714 ||
		inputMonth <= 0 || inputMonth > 12 ||
		inputDay <= 0 || inputDay > 31) {
		return (0);
	}
	/* check for dates before SDN 1 (Nov 25, 4714 B.C.) */
	if (inputYear == -4714) {
		if (inputMonth < 11) {
			return (0);
		}
		if (inputMonth == 11 && inputDay < 25) {
			return (0);
		}
	}
	/* Make year always a positive number. */
	if (inputYear < 0) {
		year = inputYear + 4801;
	} else {
		year = inputYear + 4800;
	}

	/* Adjust the start of the year. */
	if (inputMonth > 2) {
		month = inputMonth - 3;
	} else {
		month = inputMonth + 9;
		year--;
	}

	return (((year / 100) * DAYS_PER_400_YEARS) / 4
			+ ((year % 100) * DAYS_PER_4_YEARS) / 4
			+ (month * DAYS_PER_5_MONTHS + 2) / 5
			+ inputDay
			- GREGOR_SDN_OFFSET);
}

If you change the expression so that it is easier to compare with others,

\begin{aligned}
Y' & = \left\{
    \begin{aligned}
        Y + 4801 \quad (Y < 0) \\
        Y + 4800 \quad (Y \ge 0)
    \end{aligned}
\right. \\
M' & = \left\{
    \begin{aligned}
        & M - 3 & (M > 2) \\
        & M + 9, \quad Y' = Y' - 1 & (M \le 2)
    \end{aligned}
\right. \\
JD & = \text{INT} \left( \frac{ \frac{Y'}{100} \times 146097 }{4} + \frac{(Y' \% 100) \times 1461}{4} + \frac{153 M' + 2}{5} + D - 32045 \right)
\end{aligned}

pyorbital method

pyorbital/pyorbital/__init__.py

`init.py`


def dt2np(utc_time):
    try:
        return np.datetime64(utc_time)
    except ValueError:
        return utc_time.astype('datetime64[ns]')

pyorbital/pyorbital/astronomy.py

`astronomy.py`


def jdays2000(utc_time):
    """Get the days since year 2000.
    """
    return _days(dt2np(utc_time) - np.datetime64('2000-01-01T12:00'))


def jdays(utc_time):
    """Get the julian day of *utc_time*.
    """
    return jdays2000(utc_time) + 2451545


def _days(dt):
    """Get the days (floating point) from *d_t*.
    """
    return dt / np.timedelta64(1, 'D')

In other words

JD =  UTC(Y,M,D,h,m,s) - UTC(2000, 1, 1, 12, 0, 0) + 2451545

National Astronomical Observatory of Japan Web Page Method

National Astronomical Observatory of Japan / Julian Day

Dates follow the rules of the Gregorian calendar after October 15, 1582, and the Julian calendar before that. The year before the first year of the era is 0 year. For this reason, there is a year-to-year difference from the method of setting this to 1 BC.

Julian days can be calculated by filling out a web form. Output to a CSV file for comparison with the results of other calculation methods. Since it is difficult to input manually, send a request to the Web form with python, get the result and output it to CSV.

Python code

Below is a CSV graph generated by the above code. On this website, -4712/01/01 12:00 is the minimum value (Julian day = 0.0 days).

--- Around 4712 (-4712/01/01 ~ -4711/12/01)

--Around 0 years (0001/01/01 ~ -0001/12/01)

--October 1582, when the Gregorian calendar began (1582/10/01 ~ 1582/10/31)

The day after 1582/10/04, the last day of the Julian calendar, is the first day of the Gregorian calendar, 1582/10/15, so the web page is returning results in that way. In addition, the values from 10/05 to 10/14 do not exist. Perhaps because of that, 10/15 is after 10/24. Anyway, the Julian day on 10/04 is 2299160 and the Julian day on 10/15 is 2299161, and the linear result that the value does not fly is maintained.

Comparison of calculation results of each method

Python code for comparison

Horizontal axis: Gregorian calendar (I feel that it is inappropriate to call the Gregorian calendar before 1582/10/15 ...) Vertical axis: Julian day

Regarding NAOJ (calculation result on the website of the National Astronomical Observatory of Japan), since it is worried about the load and time to calculate from 4713 to 2000 in 1-month increments, only the following period is calculated, and it is shown in the graph. Is also plotted only in the following period.

--- Around 4712 (-4712/01/01 ~ -4711/12/01) --Around 0 years (0001/01/01 ~ -0001/12/01) --October 1582, when the Gregorian calendar began (1582/10/01 ~ 1582/10/31)

-4713-2000

-4713/01/01 ~ 2000/12/01 is calculated in 1-month increments (1st of each month).

From a bird's-eye view, the values are almost the same for all calculation methods, and the results are almost linear.

-Enlarged view around 4713

Calculated from -4713/01/01 to -4711/12/01 in 1-month increments (1st of each month).

Regarding NAOJ (calculation result on the National Astronomical Observatory of Japan website), -4712/01/01 12:00 is the minimum value (Julian day = 0.0 days) on this website, so it started in -4712.

Enlarged view around 0 years

Calculated from -0001/01/01 to 0001/12/01 in 1-month increments (1st of each month).

Only php is discontinuous from December 0001 to January 0000, and it is thought that php treats "1 BC is followed by 1 BC". From the Julian day value, 0 BC is the same as 1 BC. On the contrary, it seems that the method other than php is treated as "1 BC → 0 BC → 1 BC".

Enlarged view around year 0, but excluding year 0

Calculated from -0001/01/01 to -0001/12/01 and 0001/01/01 to 0001/12/01 in 1-month increments (1st day of each month).

As mentioned above, if you try to draw a graph excluding 0 years, only php will be continuous, and other methods will be discontinuous.

Enlarged view of October 1582, when the Gregorian calendar began

Calculated from 1582/10/01 to 1582/10/31 in 1-day increments.

Regarding NAOJ (calculation result on the website of the National Astronomical Observatory of Japan), the day after 1582/10/04, which is the last day of the Julian calendar, will be 1582/10/15, which is the first day of the Gregorian calendar. There is no value for / 05-10 / 14 because it is not included in the calculation of Julian day.

Difference between maximum and minimum values

-4713/01/01 ~ 2000/12/01 is calculated in 1-month increments (1st of each month).

Since the calculated Julian day values vary between calculation methods, I plotted the difference between the maximum and minimum values to see how much the variation was. --- Difference in 4713/01/01 to -0001/12/01: 417.5 to 381.5

0000/01/01 ～ 2000/12/01：　17.5 ～ 1.5 --1582/10/01: 13.5 (As mentioned above, 1582/10 calculates NAOJ, and NAOJ has the maximum value and the difference is large)

The closer to the present, the smaller the difference.

Compared conclusions

Although it was found that there are differences in the results of each calculation method and the tendency, it remains unclear which one is correct or whether it is possible to say that one method is correct in the first place.

Time in Julian day

The full date of the Julian calendar after 12:00 UT.

\begin{aligned}
JD & = JDN + \frac{hour - 12}{24} + \frac{minute}{1440} + \frac{second}{86400} \\
JD & : \text{Julian Date} \\
JDN & : \text{Julian Day Number}
\end{aligned}

wikipedia / Julian day#Finding Julian date given Julian day number and time of day

Other

Modified Julian Day: MJD

[Wikipedia / Julian Day #Modified Julian Day (MJD)](https://ja.wikipedia.org/wiki/%E3%83%A6%E3%83%AA%E3%82%A6%E3%82% B9% E9% 80% 9A% E6% 97% A5 #% E4% BF% AE% E6% AD% A3% E3% 83% A6% E3% 83% AA% E3% 82% A6% E3% 82% B9 % E6% 97% A5% EF% BC% 88MJD% EF% BC% 89) wikipedia / Julian day#Variants Meeus, J., Astronomical Algorithms, 1998

\begin{aligned}
MJD & = JD - 2400000.5 \\
MJD & :Modified Julian Day\\
JD & :Julian day
\end{aligned}

Truncated Julian Day: TJD

Which is correct, with or without the floor function?

TJD = JD - 2440000.5

wikipedia / Julian day#Variants

TJD = \left\lfloor JD - 2440000.5 \right\rfloor

J2000

wikipedia / Epoch (astronomy)#Julian Dates and J2000

date that is an interval of x Julian years of 365.25 days away from the epoch J2000 = JD 2451545.0 (TT), still corresponding (in spite of the use of the prefix "J" or word "Julian") to the Gregorian calendar date of January 1, 2000, at 12h TT (about 64 seconds before noon UTC on the same calendar day).[10](See also Julian year (astronomy).) Like the Besselian epoch, an arbitrary Julian epoch is therefore related to the Julian date by

J = 2000 + ( \text{Julian date} -2451545.0 ) \div 365.25

Other reference sites

Julian Day Playing on "Julian Day" The Julian Period

[Wikipedia / Gregorian calendar](https://ja.wikipedia.org/wiki/%E3%82%B0%E3%83%AC%E3%82%B4%E3%83%AA%E3%82%AA% E6% 9A% A6)

[Python] Acquisition of Web data by GET / POST request (Requests module) [Python] Scraping tables using Beautiful Soup Beautiful Soup in 10 minutes

ticklabels empty when not interactive #6103

[PYTHON] Julian day calculation method

Julian day calculation method

wikipedia method

Fliegel et al.'S method

Hatcher's method

Meeus method

Vallado et al.'S method

Howard D. Curtis's method

boost library method

PHP method

php:calendar.stub.ph

calendar.c

gregor.c

pyorbital method

__init__.py

astronomy.py

National Astronomical Observatory of Japan Web Page Method

Comparison of calculation results of each method

-4713-2000

-Enlarged view around 4713

Enlarged view around 0 years

Enlarged view around year 0, but excluding year 0

Enlarged view of October 1582, when the Gregorian calendar began

Difference between maximum and minimum values

Compared conclusions

Time in Julian day

Other

Modified Julian Day: MJD

Other reference sites

`php:calendar.stub.ph`

`calendar.c`

`gregor.c`

`init.py`

`astronomy.py`