ユリウス日の算出方法

ユリウス日の算出方法を調べていたところ、出典によって算出式がマチマチでどれが正しいのか困ってしまったため、とりあえず整理してみたメモです。
各算出方法で結果に差が出ることは分かったものの、どれが正しいのか、あるいは、そもそも唯一の正しい方法があるタイプの問題なのかどうか、このメモでは結局まだ不明です。

ユリウス日の算出方法

wikipedia / ユリウス通日#Julian Day Number (JDN) によると、

換算式は、Fliegel and Van Flandern[13]、Hatcher[14]、Meeus[15]によって考案されている。ただしこれらに整理を施した換算式が使われることも多い[16]。

とのことであり、式にはいくつかの形がある模様。
以下に、wikipediaでの式と、他文献等の式を列挙します。それぞれ式の形がかなり違っており、ほんとに等価なのか疑問です。

wikipediaの方法

wikipedia / 恒星時#恒星時の計算法

$$
\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は床関数
\end{aligned}
$$

UTの現在のグレゴリオ暦での年をY、月をM、日をD、時間をh、分をm、秒をsとする。ただし、1月と2月はそれぞれ前年（Yの値を-1する）の13月、14月として代入する（例: 2013年2月5日の場合、Y=2012, M=14, D=5）。

Fliegelらの方法

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.

乗算と除算は左から右に実行され、そのつど計算結果の絶対値は、その次に小さい整数値に切り下げられる。

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.

他と比較しやすいよう表現を変えると、

$$
\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の方法

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はxの整数部であり、xは正と想定される。(x) MODはxをyで除算した正の剰余。

グレゴリオ暦からの変換の場合は

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
$$

と記載されているため、上記の式に代入すると

$$
\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の方法

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

他と比較しやすいよう表現を変えると、

$$
\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らの方法

CelesTrak / Revisiting Spacetrack Report #3 – AIAA 2006-6753 – Source code (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}
$$

他と比較しやすいよう表現を変えると、

$$
\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の方法

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

上記のFourth EditionとThird Editionでは、第2項の分子の7がかかっている範囲が異なっている。

boostライブラリの方法

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);
  }

他と比較しやすいよう表現を変えると、

$$
\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の方法

PHP > マニュアル > 関数リファレンス > 日付および時刻関連 > カレンダー > カレンダー関数
gregoriantojd

gregoriantojd — グレゴリウス日をユリウス積算日に変換する

PHPの関数ではあるが、以下の通り実体はCで実装されている模様。

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

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

php-src/ext/calendar/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

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);
}

他と比較しやすいよう表現を変えると、

$$
\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の方法

pyorbital/pyorbital/__init__.py

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

pyorbital/pyorbital/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')

つまり

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

国立天文台のWebページの方法

国立天文台 / ユリウス日

1582年10月15日以後はグレゴリオ暦、それより前はユリウス暦の規則に従った日付となります。
紀元元年の前年を0年としています。このため、これを紀元前1年とする方法とは1年ずつ差異があります。

Webフォームに入力することでユリウス日を算出することができます。
他の算出方法での結果と比べるために、CSVファイルに出力しておきます。手入力だと大変なため、pythonでWebフォームへリクエストを送信して結果を取得しCSVに出力します。

Pythonコード

以下は上記のコードで生成したCSVのグラフ。なお、このWebサイトでは -4712/01/01 12:00 が最小値（ユリウス日＝0.0日）。

-4712年周辺（-4712/01/01 ～ -4711/12/01）

0年周辺（0001/01/01 ～ -0001/12/01）

グレゴリオ暦が始まった1582年10月（1582/10/01 ～ 1582/10/31）

ユリウス暦の最終日である1582/10/04の次の日はグレゴリオ暦の最初日である1582/10/15であるため、Webページもそのような形で結果を返してきています。
また、10/05～10/14という値が存在しない形になっており、その影響か、10/24の次が10/15になってしまっています。
ともかく、10/04のユリウス日は2299160、10/15のユリウス日は2299161となっており、値が飛ぶことのなく線形に値が維持されています。

各方法の算出結果の比較

比較用のPythonコード

横軸：グレゴリオ暦（1582/10/15よりも前をグレゴリオ暦と呼ぶのは不適切な気もしますが…）
縦軸：ユリウス日

なお、NAOJ（国立天文台のWebサイトでの算出結果）については、-4713年～2000年を1ヶ月刻みで算出するのは負荷と時間が心配なので以下の期間だけを算出しているため、グラフにも以下の期間でしかプロットされていません。

-4712年周辺（-4712/01/01 ～ -4711/12/01）
0年周辺（0001/01/01 ～ -0001/12/01）
グレゴリオ暦が始まった1582年10月（1582/10/01 ～ 1582/10/31）

-4713年～2000年

-4713/01/01 ～ 2000/12/01 を1ヶ月刻み（各月の1日）で算出。

俯瞰するとどの算出方法でもほぼ同じ値で、ほぼ線形な結果になっています。

-4713年周辺の拡大図

-4713/01/01 ～ -4711/12/01 を1ヶ月刻み（各月の1日）で算出。

NAOJ（国立天文台のWebサイトでの算出結果）については、このWebサイトでは -4712/01/01 12:00 が最小値（ユリウス日＝0.0日）のため、-4712年から始まっています。

0年周辺の拡大図

-0001/01/01 ～ 0001/12/01 を1ヶ月刻み（各月の1日）で算出。

phpのみ、-0001年12月から0000年1月で不連続となっており、phpでは「紀元前1年の次は紀元後1年」という扱いをしていると考えられます。ユリウス日の値から見れば、紀元0年は紀元前1年と同じということになります。
逆にphp以外の方法では、「紀元前1年→紀元0年→紀元後1年」という扱いをしているようです。

0年周辺の拡大図、ただし0年は除外して描画

-0001/01/01 ～ -0001/12/01、0001/01/01 ～ 0001/12/01 を1ヶ月刻み（各月の1日）で算出。

前述の通り、ためしに0年を除外してグラフを描画するとphpだけが連続になり、他の方法は不連続になります。

グレゴリオ暦が始まった1582年10月の拡大図

1582/10/01 ～ 1582/10/31 を1日刻みで算出。

NAOJ（国立天文台のWebサイトでの算出結果）については、10/05～10/14がユリウス日の算出対象外になっているようで、値がありません。
これは、ユリウス暦の最終日である1582/10/04の次の日はグレゴリオ暦の最初日である1582/10/15になるということに準拠しているためと思われます。

最大値と最小値の差

-4713/01/01 ～ 2000/12/01 を1ヶ月刻み（各月の1日）で算出。

算出されたユリウス日の値は算出方法間でバラツキがあるため、バラツキがどれくらいかを見るために最大値と最小値の差をプロットしてみました。

-4713/01/01 ～ -0001/12/01 における差：　417.5 ～ 381.5
0000/01/01 ～ 2000/12/01：　17.5 ～ 1.5
1582/10/01：　13.5（前述の通り1582/10はNAOJを算出しており、NAOJが最大値を出していて差が大きくなっている）

現在に近づくほど差は小さくなってきています。

比較した結論（未完）

各算出方法で結果に差が出ることは分かったものの、どれが正しいのか、あるいは、そもそも唯一の正しい方法があるタイプの問題なのかどうか、結局まだよく分からないという状況です。

ユリウス日における時刻

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

その他

Modified Julian Day: MJD（修正ユリウス日）

wikipedia / ユリウス通日#修正ユリウス日（MJD）
wikipedia / Julian day#Variants
Meeus, J., Astronomical Algorithms, 1998

$$
\begin{aligned}
MJD & = JD - 2400000.5 \\
MJD & : 修正ユリウス日 \\
JD & : ユリウス日
\end{aligned}
$$

Truncated Julian Day: TJD

床関数のあるのと無いの、どちらが正しいのか？

wikipedia / 恒星時#恒星時の計算法

$$
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
$$

他参考サイト

wikipedia / グレゴリオ暦

ticklabels empty when not interactive #6103

ユリウス日の算出方法

ユリウス日の算出方法

wikipediaの方法

Fliegelらの方法

Hatcherの方法

Meeusの方法

Valladoらの方法

Howard D. Curtisの方法

boostライブラリの方法

PHPの方法

pyorbitalの方法

国立天文台のWebページの方法

各方法の算出結果の比較

-4713年～2000年

-4713年周辺の拡大図

0年周辺の拡大図

0年周辺の拡大図、ただし0年は除外して描画

グレゴリオ暦が始まった1582年10月の拡大図

最大値と最小値の差

比較した結論（未完）

ユリウス日における時刻

その他

Modified Julian Day: MJD（修正ユリウス日）

Truncated Julian Day: TJD

J2000

他参考サイト

いいね:

コメントを残すコメントをキャンセル

ユリウス日の算出方法

ユリウス日の算出方法

wikipediaの方法

Fliegelらの方法

Hatcherの方法

Meeusの方法

Valladoらの方法

Howard D. Curtisの方法

boostライブラリの方法

PHPの方法

pyorbitalの方法

国立天文台のWebページの方法

各方法の算出結果の比較

-4713年～2000年

-4713年周辺の拡大図

0年周辺の拡大図

0年周辺の拡大図、ただし0年は除外して描画

グレゴリオ暦が始まった1582年10月の拡大図

最大値と最小値の差

比較した結論（未完）

ユリウス日における時刻

その他

Modified Julian Day: MJD（修正ユリウス日）

Truncated Julian Day: TJD

J2000

他参考サイト

共有:

いいね:

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル