グレゴリオ暦と修正ユリウス日

はじめに

 ある日からある日までの日数計算や特定の日の曜日計算は、まれに必要になることがあります。 専用の関数があれば簡単なのですが、HSPには日付の計算をする命令が用意されていません。 無いものは作るしか無いのですが、作ろうとすると実装は容易ではありません。 ここでは日付計算に必要な技術とHSP3での実装例を順を追って解説してみます。

 なお、ここで解説した実装例はHSP用のモジュールとして配布していますので、必要な方はモジュールをダウンロードしてご利用ください。  →修正ユリウス日換算モジュール

 ではまずは必要な用語の解説から。

グレゴリオ暦

 今日、我々が一般に使用している暦のこと。
 例えばこれを書いている今日は西暦「2017年 2月25日」です。これがグレゴリオ暦です。 1582年からユリウス暦の後に採用された暦だそうです。 ユリウス暦のことをあわせて書くと話がややこしくなるので特別な記述がない場合は、ここでは1582年以前もグレゴリオ暦表記することにします。

 グレゴリオ暦は、時刻の単位が組み合わさっていて、さらに時刻はから出来ています。 さらにさらにはごとに最大値がばらばら。 このバラバラの単位が問題で、計算しようとすると頭を抱えることになります。 これに加えて年は年で、4年毎にうるう年がありますし、それも100年ごとに適用されない年がある。でも400年ごとには適用されるという…超面倒くさい仕様!
 これでは日数計算するのは大変です。計算を簡単にするには全て同じ1個の単位で扱えることが必要です。 それを可能にするのが「ユリウス通日」です。

ユリウス通日(ユリウス日)(Julian Day、JD)

 紀元前4714年11月24日(ユリウス暦 紀元前4713年1月1日)の正午からの日数。 時刻値は小数値が用いられます。 これのお陰で、日時が1個の単位、1つの数値になるので日数計算がグンと簡単になりました。

 しかしこれも問題があります。基準となる0日が昔過ぎて、最近の日付を表そうとした時に桁数が膨大になってしまいます。 これは昔の記憶容量が小さいコンピューターでは大きな問題となったようで修正ユリウス日が開発されました。

修正ユリウス日(Modified Julian Date、MJD)

 修正ユリウス日は、ユリウス日から2400000.5を差し引いたもの。1858年11月17日が修正ユリウス日 0。 ユリウス日を使って最近の日付を計算しようとすると無駄に桁が大きくなります。6000年以上も前の日付が基準ですからね。 これが昔のコンピュータには不向きだったことから、1957年に修正ユリウス日が考案されました。

 ここではユリウス通日ではなく、修正ユリウス日について記述していきます。PC向きのようなので。

紀元前

 西暦よりも昔のこと。
 西暦1年の前年は紀元前1年です。 過去に行くほど年の数字が大きくなっていきます。 西暦0年や紀元前0年は存在しません。これはまた計算がめんどくさそうですね。
 なお、計算では便宜上、紀元前1年は0として計算に使用します。紀元前2年は-1ですね。

ということで用語解説はここまでです。

修正ユリウス日の例

 計算にする前に、計算結果が正しいか確認する方法を確保しておく必要があります。

 検索すると計算してくれるサイトがすぐに見つかります。 単に知りたいだけなら手軽な下記サイトをお勧めします。

高精度計算サイト keisan
http://keisan.casio.jp/exec/system/1177639469

国立天文台(NAOJ)
http://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/date2jd.cgi

 上記のサイトを使って計算した結果がこちらになります。 この数字を使って自分の計算結果があっているかを確認しながら作っていきます。

BC/AC修正ユリウス日備考
BC-4712 1 1-2399963BC 4713年
AC 200 3 1 -605833
AC 15821015 -100840
AC 2000 1 1 51544
AC 1 1 1 -678575
BC 01231 -678576BC 1年
AC 18581117 0
AC 1886 4 4 10000
AC 1913 821 20000
AC 1831 7 2 -10000
AC 1000 1 1 -313698

フリーゲルの公式

 グレゴリオ暦→修正ユリウス日 の変換を行う計算式というものがあります。 今回はこれを使用します。式はWikipediaに載っています。

グレゴリオ暦 year 年 month 月 day 日午前0時の修正ユリウス日は次式で表される。

MJD = [365.25y] + [y/400] ー [y/100] + [30.59(m-2)] + d - 678912

  • [x]は床関数。
  • 1月、2月は前年の13月、14月として計算する。
if (month=1月) or (month=2月) {
	y = year - 1
	m = month + 12
} else {
	y = year
	m = month
}
d = day

 という感じです。式だけ見ると一見簡単そうですよね。 とりあえず動くものなら実装は簡単で、[ ]はint関数(小数点以下切り捨て)で書けば紀元後の計算であれば問題なく動作します。 でも実際には床関数なので紀元前を計算しようとすると計算誤差が出始めます。簡単そうに見えたんですけどね…。

なお、y,m,dの変換は分かりやすいようにif文で書きましたが、HSP3での実装はこんな感じです。

m = mod(Month-3,12) + 3
d = Day
y = double( Year + floor((-3.0+Month)/12.0) )

床関数

 詳しくは検索して調べてください。検索したほうが分かりやすい解説がすぐに複数見つかります。 HSP的には「負の方向へ丸め」。hspmathのfloor関数がこの計算をやってくれます。

 これで紀元前の日付から求めた修正ユリウス日も正しい結果が得られるようになったはずです。 次は逆に修正ユリウス日から元の年月日を求めてみます。

修正ユリウス日から年月日へ

 修正ユリウス日から元の年月日に戻せないと使い物になりません。 計算式はWikipediaに載っているので見てみましょう。

ユリウス通日 グレゴリオ暦からの換算式 - Wikipedia
https://ja.wikipedia.org/wiki/ユリウス通日

 ひとつ下の「逆換算」の項目が修正ユリウス日から年月日へ変換する換算式です。 この式で出ているn,y,m,dは「グレゴリオ暦からの換算式」に書かれています。 フリーゲルの公式の項で説明に使用しているy,m,dとは異なる数値なので注意してください。(Wikipediaの表記が紛らわしい!)

y,m,dから年月日にするのはこのように実装しました。

Day   = d + 1
Month = (m+2)¥12+1
Year  = y - floor((-3.0+Month)/12.0)

 相互に変換できるようになりました。これでようやく実践に使えるようになりましたね! …と、思いますよね。思いましたよッ!くそーッ!まだ計算結果合わないよ!!

 記載された計算式をHSPに書き写しても正しい年月日は得られません。 乗除の扱いが違うようです。次は乗除について考えます。

乗除

 割り算の余りの計算のことです。modulo、modとか書いたりもします。
HSPで割り算の余りを計算すると次のようになります。
5¥3  → 2
-5¥3 → -2
5¥-3 → 2

 CやC++、Javaでもこのような動作になるそうです。 しかし、実際に期待されている計算結果はこうです。
5¥3  → 2
-5¥3 → 1
5¥-3 → -1

 負数が入ると結果が違いますね。RubyやPythonはこのような結果を返すそうです。 逆換算の式の途中にmodがあるからと言ってこれを¥(HSP3の乗除)で書いてしまうと、正しい計算結果が得られないわけです。 困ったものです。

このあたりの内容については下記サイトに説明がありますので読んでみてください。
おともだちティータイム
http://shunirr.hatenablog.jp/entry/20120409/1333993409

 あきらめるわけにも行かず、また代用できる関数も無いようなのでここは自作するしかありません。 そこでこのようにHSPで実装してみました。適当ですがとりあえず期待通り動くようなのでよしとしましょう。

#define ctype mod(%1,%2) (((%1)-(int((%1)/(%2))-1)*(%2)) ¥ (%2))

結果:
mod( 5,3) → 2
mod(-5,3) → 1
mod(5,-3) → -1

 これでHSPでも期待したような割り算の余りの値を得られるようになりました。 作ったこの機能を先程の式に組み込むと、修正ユリウス日から正しい年月日が取得できるようになったと思います。

まとまってないまとめ

 最初はWikipediaに計算式載ってるからそれ使えば簡単だろ。ぐらいに考えてたんですが考えが甘かったです。 しかしおかげで床関数や負数の乗除など新しい知識も得られましたし、計算結果の精度もムダに高いものが出来ました。 さて、次は作ったこのモジュールを使って…

 何も作る予定がない。

 そもそも最近作ってたのはもっと別のもので、日付なんか使わないツール作ってて、このモジュールは息抜きだったんでした。 元の作業に戻るか…。

関連記事

  1. 修正ユリウス日換算モジュール ダウンロード mod_mjd Ver.1.00 (2017/...
  2. 測地座標を直交座標に変換 目次 はじめに 地図データを入手 計算式を入手 用語解説 測...
  3. 測地座標をメルカトル図法に変換 目次 あらすじ メルカトル図法 換算式 ミラー図法 換算式 ...