Ruby - 「グレゴリオ暦 -> ユリウス日」変換の検証！

２種の計算式で「グレゴリオ暦 -> ユリウス日」の変換を行い、結果が同じになるかを検証してみました。

0. 前提条件

• Ruby 2.3.1-p112 での作業を想定。
• 使用する２種の計算式は「フリーゲルの公式」とその他の計算式。
• ここでの「ユリウス日」は JD(Julian Day) であり、JDN(Julian Day Number), CJD(Chronological Julian Day), MJD(Modified Julian Date) ではない。

1. 検証用 Ruby スクリプトの作成

以下のように作成してみた。（うるう年の２月から３月への変わり目、年末から年始への変わり目、うるう年でない年の２月から３月への変わり目をチェックするようにしている）

File: verify_gc2jd.rb

#! /usr/local/bin/ruby
# coding: utf-8
#=グレゴリオ暦 -> ユリウス日 変換の検証(UTC ベース)
#
# : ２種類の計算式で グレゴリオ暦 -> ユリウス日 の変換を行い、
#   相違がないかを検証する。
#
# date          name            version
# 2016.06.14    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------

class VerifyGc2Jd
  def run
    begin
      55.upto(59) { |sec| generate(2016,  2, 28, 23, 59, sec) }
       0.upto( 5) { |sec| generate(2016,  2, 29,  0,  0, sec) }
      puts "..."
      55.upto(59) { |sec| generate(2016,  2, 29, 23, 59, sec) }
       0.upto( 5) { |sec| generate(2016,  3,  1,  0,  0, sec) }
      puts "..."
      55.upto(59) { |sec| generate(2016, 12, 31, 23, 59, sec) }
       0.upto( 5) { |sec| generate(2017,  1,  1,  0,  0, sec) }
      puts "..."
      55.upto(59) { |sec| generate(2017,  2, 28, 23, 59, sec) }
       0.upto( 5) { |sec| generate(2017,  3,  1,  0,  0, sec) }
    rescue => e
      msg = "[#{e.class}] #{e.message}\n"
      msg << e.backtrace.each { |tr| "\t#{tr}" }.join("\n")
      $stderr.puts msg
    end
  end

  private

  def generate(*dt)
    str =  Time.new(*dt).strftime("%Y-%m-%d %H:%M:%SUTC, JD: ")
    jd_1, jd_2 = gc2jd1(*dt), gc2jd2(*dt)
    str << "#{jd_1}, #{jd_2}"
    str << " (NOT MATCH!)" unless jd_1 == jd_2
    puts str
  rescue => e
    raise
  end

  #=========================================================================
  # グレゴリオ暦 -> ユリウス日(JD)
  #
  # * フリーゲルの公式を使用する
  #     month = 1, 2 の場合、year = year - 1, month = month + 12 とする
  #     JD = int(365.25 * year)
  #        + int(year / 400)
  #        - int(year / 100)
  #        + int(30.59 * (month - 2))
  #        + day
  #        + 1721088.5
  #   ※上記の int(x) は厳密には、x を超えない最大の整数
  #   ※正子を整数とする JDN を求めるなら、 1721088.5 を 1721089 に置き換える
  #   ※修正ユリウス日 MJDを求めるなら + 1721088.5 を - 678912 に置き換える。
  #
  # @param:  year
  # @param:  month
  # @param:  day
  # @param:  hour
  # @param:  minute
  # @param:  second
  # @return: JD
  #=========================================================================
  def gc2jd1(year, month, day, hour, min, sec)
    # 1月,2月は前年の13月,14月とする
    if month < 3
      year  -= 1
      month += 12
    end
    # 日付(整数)部分計算
    jd  = (365.25 * year).truncate \
        + (year / 400.0).truncate \
        - (year / 100.0).truncate \
        + (30.59 * (month - 2)).truncate \
        + day \
        + 1721088.5
    # 時間(小数)部分計算
    t  = (hour + min / 60.0 + sec / 3600.0) / 24.0
    return sprintf("%18.10f", jd + t)
  rescue => e
    raise
  end

  #=========================================================================
  # グレゴリオ暦 -> ユリウス日(JD)
  #
  # * フリーゲルの公式ではない計算式を使用する
  #   [The Jurian Pieiod](http://www.tondering.dk/claus/cal/julperiod.php#formula)
  #     a = int((14 - month) / 12)
  #     y = year + 4800 - a
  #     m = month + 12a - 3
  #    JD = day + int((153m + 2) / 5) + 365y
  #       + int(y / 4) - int(y / 100) - int (y / 400) - 32045.5
  #   ※上記の int(x) は厳密には、x を超えない最大の整数
  #   ※正子を整数とする JDN を求めるなら、 32045.5 を 32045 に置き換える
  #   ※修正ユリウス日 MJD を求めるなら、 32045.5 を 2432046 に置き換える
  #
  # @param:  year
  # @param:  month
  # @param:  day
  # @param:  hour
  # @param:  minute
  # @param:  second
  # @return: JD
  #=========================================================================
  def gc2jd2(year, month, day, hour, min, sec)
    # 日付(整数)部分計算
    a = ((14 - month) / 12.0).truncate
    y = year + 4800 - a
    m = month + 12 * a - 3
    jd = day \
       + ((153 * m + 2) / 5.0).truncate \
       + 365 * y \
       + (y / 4.0).truncate \
       - (y / 100.0).truncate \
       + (y / 400.0).truncate \
       - 32045.5
    # 時間(小数)部分計算
    t  = (hour + min / 60.0 + sec / 3600.0) / 24.0
    return sprintf("%18.10f", jd + t)
  rescue => e
    raise
  end
end

VerifyGc2Jd.new.run if __FILE__ == $0

• Gist - Ruby script to verify of converting Gregorian Calendar to Julian Day.

2. 検証用 Ruby スクリプトの実行

$ ./verify_gc2jd.rb
2016-02-28 23:59:55UTC -> JD: 2457447.4999421295, 2457447.4999421295
2016-02-28 23:59:56UTC -> JD: 2457447.4999537035, 2457447.4999537035
2016-02-28 23:59:57UTC -> JD: 2457447.4999652780, 2457447.4999652780
2016-02-28 23:59:58UTC -> JD: 2457447.4999768520, 2457447.4999768520
2016-02-28 23:59:59UTC -> JD: 2457447.4999884260, 2457447.4999884260
2016-02-29 00:00:00UTC -> JD: 2457447.5000000000, 2457447.5000000000
2016-02-29 00:00:01UTC -> JD: 2457447.5000115740, 2457447.5000115740
2016-02-29 00:00:02UTC -> JD: 2457447.5000231480, 2457447.5000231480
2016-02-29 00:00:03UTC -> JD: 2457447.5000347220, 2457447.5000347220
2016-02-29 00:00:04UTC -> JD: 2457447.5000462965, 2457447.5000462965
2016-02-29 00:00:05UTC -> JD: 2457447.5000578705, 2457447.5000578705
...
2016-02-29 23:59:55UTC -> JD: 2457448.4999421295, 2457448.4999421295
2016-02-29 23:59:56UTC -> JD: 2457448.4999537035, 2457448.4999537035
2016-02-29 23:59:57UTC -> JD: 2457448.4999652780, 2457448.4999652780
2016-02-29 23:59:58UTC -> JD: 2457448.4999768520, 2457448.4999768520
2016-02-29 23:59:59UTC -> JD: 2457448.4999884260, 2457448.4999884260
2016-03-01 00:00:00UTC -> JD: 2457448.5000000000, 2457448.5000000000
2016-03-01 00:00:01UTC -> JD: 2457448.5000115740, 2457448.5000115740
2016-03-01 00:00:02UTC -> JD: 2457448.5000231480, 2457448.5000231480
2016-03-01 00:00:03UTC -> JD: 2457448.5000347220, 2457448.5000347220
2016-03-01 00:00:04UTC -> JD: 2457448.5000462965, 2457448.5000462965
2016-03-01 00:00:05UTC -> JD: 2457448.5000578705, 2457448.5000578705
...
2016-12-31 23:59:55UTC -> JD: 2457754.4999421295, 2457754.4999421295
2016-12-31 23:59:56UTC -> JD: 2457754.4999537035, 2457754.4999537035
2016-12-31 23:59:57UTC -> JD: 2457754.4999652780, 2457754.4999652780
2016-12-31 23:59:58UTC -> JD: 2457754.4999768520, 2457754.4999768520
2016-12-31 23:59:59UTC -> JD: 2457754.4999884260, 2457754.4999884260
2017-01-01 00:00:00UTC -> JD: 2457754.5000000000, 2457754.5000000000
2017-01-01 00:00:01UTC -> JD: 2457754.5000115740, 2457754.5000115740
2017-01-01 00:00:02UTC -> JD: 2457754.5000231480, 2457754.5000231480
2017-01-01 00:00:03UTC -> JD: 2457754.5000347220, 2457754.5000347220
2017-01-01 00:00:04UTC -> JD: 2457754.5000462965, 2457754.5000462965
2017-01-01 00:00:05UTC -> JD: 2457754.5000578705, 2457754.5000578705
...
2017-02-28 23:59:55UTC -> JD: 2457813.4999421295, 2457813.4999421295
2017-02-28 23:59:56UTC -> JD: 2457813.4999537035, 2457813.4999537035
2017-02-28 23:59:57UTC -> JD: 2457813.4999652780, 2457813.4999652780
2017-02-28 23:59:58UTC -> JD: 2457813.4999768520, 2457813.4999768520
2017-02-28 23:59:59UTC -> JD: 2457813.4999884260, 2457813.4999884260
2017-03-01 00:00:00UTC -> JD: 2457813.5000000000, 2457813.5000000000
2017-03-01 00:00:01UTC -> JD: 2457813.5000115740, 2457813.5000115740
2017-03-01 00:00:02UTC -> JD: 2457813.5000231480, 2457813.5000231480
2017-03-01 00:00:03UTC -> JD: 2457813.5000347220, 2457813.5000347220
2017-03-01 00:00:04UTC -> JD: 2457813.5000462965, 2457813.5000462965
2017-03-01 00:00:05UTC -> JD: 2457813.5000578705, 2457813.5000578705

結果として、一致した。

---

グレゴリオ暦からユリウス日を計算するのによく「フリーゲルの公式」を使用するが、別の計算式でも問題ないことが分かりました。

Web 上で検索すると他にも計算式がヒットしますが、どれも同じでしょうね。（単に計算式をどう変形させるかだけの違い）

以上。