十進 BASIC で MJD
〜 修正ユリウス日を計算 〜
2023-07-07 作成 福島
TOP > asob > decbasic-mjd
[ TIPS | TOYS | OTAKU | LINK | MOVIE | CGI | AvTitle | ConfuTerm | HIST | AnSt | Asob | Shell ]

日付の加減算に関する問題

去年 (1 年前) の日付は「今年の西暦 - 1」を計算すれば求められます。
来年 (1 年後) の日付は「今年の西暦 + 1」ですね。
例: 2023-7-7 の 1 年前は 2022-7-7。

同様に、昨日の日付は「今日 - 1」で求められます。
例: 2023-7-7 の 1 日前は 2023-7-6。

では、10,950 日前の日付はどうでしょうか。(10,950 = 365 × 30)
本稿記述時の日付が 2023-7-7 なので、30 年前の日付は 1993-7-7。
···というわけには行きません。

1993 年 ～ 2023 年の間には「うるう年」が 7 回もあります。(1996, 2000, 2004, 2008, 2012, 2016, 2020)
このため、10,950 日前の日付は 1993-7-14 となります。

うるう年には、以下のような複雑なルールがあります。
  1. 4 年に一度はうるう年にする。
  2. 100 年に一度はうるう年にしない。
  3. 400 年に一度はうるう年にする。
西暦 2000 年は、上記 2 でうるう年にしないはずが、同 3 によってうるう年になります。
年月日の加減算に関する問題
時刻の計算をしてみます。

2 時間 35 分の映画が 18 時 30 分から上映されたら帰宅できるでしょうか。(例: シン・エヴァンゲリオン劇場版:||)
(途中で退席するとかはナシで)
答えは 21:05 なので、終電には余裕がありそうということが分かります。
未成年者の場合は、保護者の帯同が必要になるでしょう。

時計では 18 時 30 分 + 2 時間 35 分 = 21 時 5 分という計算をします。

同様に日付を計算して 2023-7-7 の 10 年 8 カ月 30 日後の日付はどうなるでしょう。
答えは 2033-3-4 ですが、上記うるう年のルールがあるため、即座に計算できる人はかなり少ないと思います。
(加算する 10 年 8 カ月 30 日間の中にもうるう年があることに注意)
それでも日付を計算したい
日付が暦(れき)のままだと計算が複雑になります。
一旦、基準日をかなり昔に設定して、そこから何日経過したか日数を揃えることにより、計算を単純にします。
その便利な基準日の定義が「修正ユリウス日」(MJD: Modified Julian Day) です。
「ユリウス通日」がオリジナルですが、そのままだと使いにくいので、これを使いやすくしたのが「修正ユリウス日」です。
「ユリウス～」という名称ですが、シーザーとは関係ありません。考案者は古典・年代学者スカリゲルで、その父 (天文・哲学者) の名前です。

すべての日付を「修正ユリウス日」という単純な整数に変換し、加減算が済んだのちに通常日付 (グレゴリオ暦) に戻す、という使い方をします。

十進 BASIC で実装

アルゴリズムは Wikipedia からのパクリです。
十進 BASIC の場合、INT() の挙動が他の多くのプログラミング言語と異なり、床関数 (負数の扱いが異なる) となります。
Python や C 言語等では int() が小数部の切り捨てになります。VB に移植するときは Math.Floor() を使ってください。

十進 BASIC の関数は、一度に複数の値を返すことができないため、
ここでは 3 つの値 y,m,d (それぞれ、年,月,日を表す) をひとつの数値に詰め込んで返しています。
呼び元では、返戻値を受領してから 3 つの値に分解しています。

y,m,d は整数です。時刻を考慮していません。(Wikipedia でも整数)

check-mjd.BAS
DECLARE EXTERNAL FUNCTION greg2mjd
DECLARE EXTERNAL FUNCTION jd2mjd
DECLARE EXTERNAL FUNCTION mjd2jd
DECLARE EXTERNAL FUNCTION mjd2greg
DECLARE EXTERNAL FUNCTION mjd2week
DECLARE EXTERNAL FUNCTION mjd2zodiac
DECLARE EXTERNAL FUNCTION mjd2stems

LET y = 2012
LET m=1
LET d=1
PRINT "Gregorian y,m,d"; y; m; d

LET mjd = greg2mjd(y,m,d)
LET w = mjd2week(mjd)
PRINT "mjd";mjd;mid$("WedThuFriSatSunMonTue",1+w*3,3)

LET ymd = mjd2greg(mjd)
! 返戻値を 3 つに分解する。
LET y = INT(ymd / 10000)
LET m = MOD(INT(ymd/100),100)
LET d = MOD(ymd,100)
PRINT "Gregorian y,m,d"; y; m; d

PRINT "---"

LET y = 1582
LET m = 2
LET d = 1
PRINT "Julian y,m,d"; y; m; d

LET mjd = jd2mjd(y,m,d)
PRINT "mjd"; mjd

LET ymd = mjd2jd(mjd)
! 返戻値を 3 つに分解する。
LET y = INT(ymd / 10000)
LET m = MOD(INT(ymd/100),100)
LET d = MOD(ymd,100)
PRINT "Julian y,m,d"; y; m; d

END



MERGE "mjd.BAS"

mjd.BAS (上記 check-mjd.BAS と同じ場所に保存する)
! mjd.BAS version 1.0 written by fuku@rouge.gr.jp
!
! ユリウス日 / 修正ユリウス日 / グレゴリオ暦 の相互変換を行う。
!
! https://ja.wikipedia.org/wiki/ユリウス通日#修正ユリウス日（MJD）
! より。

! Modified Julian Day (MJD)
!
!Mutual conversion:
!    Julian Date / Modified Julian Day / Gregorian Calendar



! Convert Gregorian-calendar to MJD.
EXTERNAL FUNCTION greg2mjd(y,m,d)
   IF m <= 2 THEN
      LET m = m + 12
      LET y = y - 1
   END IF

   LET greg2mjd = INT(365.25*y) + INT(y/400) - INT(y/ 100) &
& + INT(30.59*(m-2)) + d - 678912
END FUNCTION



! Convert Julian-date to MJD.
EXTERNAL FUNCTION jd2mjd(y,m,d)
   IF m <= 2 THEN
      LET m = m + 12
      LET y = y - 1
   END IF

   LET jd2mjd = INT(365.25*y) + INT(30.59*(m-2)) + d - 678914
END FUNCTION



! Convert MJD to Julian-date.
EXTERNAL FUNCTION mjd2jd(days)
   LET n = days + 678883
   LET a = 4 * n + 3
   LET b = 5 * INT(MOD(a,1461)/4) + 2

   LET y = INT(a/1461)
   LET m = INT(b/153) + 3
   LET d = INT(MOD(b,153)/5) + 1

   IF m > 12 THEN
      LET m = m - 12
      LET y = y + 1
   END IF
LET mjd2jd = y * 10000 + m * 100 + d

END FUNCTION



! Convert MJD to Gregorian-date.
EXTERNAL FUNCTION mjd2greg(days)
   LET n = days + 678881
   LET a = 4 * n + 3 + 4 * INT(3 / 4 * INT(4*(n+1)/146097+1))
   LET b = 5 * INT(MOD(a,1461)/4) + 2

   LET y = INT(a/1461)
   LET m = INT(b/153) + 3
   LET d = INT(MOD(b,153)/5) + 1

   IF m > 12 THEN
      LET m = m - 12
      LET y = y + 1
   END IF

   LET mjd2greg = y * 10000 + m * 100 + d
END FUNCTION



! Get day of week from MJD.
EXTERNAL FUNCTION mjd2week(days)
   LET mjd2week = MOD(days, 7)      ! 0:Wed, 1:Thu, .., 6:Tue
END FUNCTION



! Get Oriental-Zodiac from MJD.
EXTERNAL FUNCTION mjd2zodiac(days)
   LET mjd2zodiac = MOD(days, 12)   ! 0:Tora, 1:U, 2:Tatsu, .., 11:Ushi
END FUNCTION



! Get Heavenly-Stems from MJD.
EXTERNAL FUNCTION mjd2stems(days)
   LET mjd2stems = MOD(days, 10)    ! 0:Kinoe, 1:Kinoto, .., 9:Mizunoto
END FUNCTION