Ruby script to calculate GMST with IAU1982 theory.

calc_gmst_iau_1982.rb

#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= グリニジ平均恒星時 GMST(= Greenwich Mean Sidereal Time)の計算
#  : IAU1982 版
#
#
# Date          Author          Version
# 2018.06.15    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2018 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 引数 : 日時(UT1（世界時1）)
#          書式：YYYYMMDD or YYYYMMDDHHMMSS
#          無指定なら現在(システム日時)を UT1 とみなす。
#---------------------------------------------------------------------------------
#++
require 'date'

class CalcGmstIau82
  PI2 = Math::PI * 2    # => 6.283185307179586
  D2R = Math::PI / 180  # => 0.017453292519943295

  def initialize
    @ut1 = get_arg
  end

  #===========================================================================
  # Execution
  #===========================================================================
  def exec
    begin
      puts @ut1.strftime("%Y-%m-%d %H:%M:%S UT1")
      jd_ut1 = gc2jd(@ut1)
      puts "JD(UT1): #{jd_ut1}"
      gmst = calc_gmst(jd_ut1)
      gmst_d = gmst * 180 / Math::PI
      gmst_h = deg2hms(gmst_d)
      puts "GMST = #{gmst} rad.\n" +
           "     = #{gmst_d} deg.\n" +
           "     = #{gmst_h}"
    rescue => e
      msg = "[#{e.class}] #{e.message}\n"
      msg << e.backtrace.map { |tr| "\t#{tr}"}.join("\n")
      $stderr.puts msg
      exit 1
    end
  end

  private

  #===========================================================================
  # コマンドライン引数の取得
  # * コマンドライン引数無指定なら現在日時
  #
  # @return: Time オプジェクト
  #===========================================================================
  def get_arg
    unless arg = ARGV.shift
      t = Time.now
      return Time.new(t.year, t.month, t.day, t.hour, t.min, t.sec)
    end
    exit 0 unless arg =~ /^\d{8}$|^\d{14}$/
    year, month, day = arg[ 0, 4].to_i, arg[ 4, 2].to_i, arg[ 6, 2].to_i
    exit 0 unless Date.valid_date?(year, month, day)
    hour, min,   sec = arg[ 8, 2].to_i, arg[10, 2].to_i, arg[12, 2].to_i
    exit 0 if hour > 23 || min > 59 || sec > 59
    return Time.new(year, month, day, hour, min, sec)
  rescue => e
    raise
  end

  #===========================================================================
  # 年月日(グレゴリオ暦) -> ユリウス日(JD) 変換
  # : フリーゲルの公式を使用
  #     JD = int(365.25 * year)
  #        + int(year / 400)
  #        - int(year / 100)
  #        + int(30.59 * (month - 2))
  #        + day
  #        + 1721088
  #   ※上記の int(x) は厳密には、x を超えない最大の整数
  #     (ちなみに、準JDを求めるなら `+ 1721088` が `- 678912` となる)
  #
  # @param:  ut1 (UT1)
  # @return: jd (Julian Day)
  #===========================================================================
  def gc2jd(ut1)
    year, month, day = ut1.year, ut1.month, ut1.day
    hour, min, sec   = ut1.hour, ut1.min, ut1.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  = (sec / 3600.0 \
       + min / 60.0 \
       + hour) / 24.0
    return jd + t
  rescue => e
    raise
  end

  #===========================================================================
  # GMST（グリニッジ平均恒星時）計算
  # : IAU1982理論(by David Vallado)によるもの
  #     GMST = 18h 41m 50.54841s
  #          + 8640184.812866s T + 0.093104s T^2 - 0.0000062s T^3 
  #     (但し、 T = 2000年1月1日12時(UT1)からのユリウス世紀単位)
  #
  # @param:  jd_ut1 (UT1 に対するユリウス日
  # @return: gmst (グリニッジ平均恒星時(単位:radian))
  #===========================================================================
  def calc_gmst(jd_ut1)
    t_ut1= (jd_ut1 - 2451545.0) / 36525.0
    gmst =  67310.54841 \
         + (876600.0 * 3600.0 + 8640184.812866 \
         + (0.093104 \
         -  6.2e-6 * t_ut1) * t_ut1) * t_ut1
    gmst = (gmst * D2R / 240.0) % PI2
    gmst += PI2 if gmst < 0.0
    return gmst
  rescue => e
    raise
  end

  #===========================================================================
  # 99.999° -> 99h99m99s 変換
  #
  # @param:  deg  (Degree)
  # @return: "99 h 99 m 99.999 s"
  #===========================================================================
  def deg2hms(deg)
    sign = ""

    begin
      h = (deg / 15.0).truncate
      _m = (deg - h * 15.0) * 4.0
      m = _m.truncate
      s = (_m - m) * 60.0
      if s < 0
        s *= -1
        sign = "-"
      end
      return "%s%2d h %02d m %06.3f s" % [sign, h, m, s]
    rescue => e
      raise
    end
  end
end

CalcGmstIau82.new.exec if __FILE__ == $0