Ruby script to calculate JPL ephemeris with a selfmade gem library.

ephemeris_jpl.rb

#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読み込み、座標（位置・速度）を計算する
#  : 自作の RubyGems ライブラリ eph_jpl を使用する
#
# date          name            version
# 2016.06.19    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 【引数】
#   [第１] 対象天体番号（必須）
#           1: 水星            (Mercury)
#           2: 金星            (Venus)
#           3: 地球            (Earth)
#           4: 火星            (Mars)
#           5: 木星            (Jupiter)
#           6: 土星            (Saturn)
#           7: 天王星          (Uranus)
#           8: 海王星          (Neptune)
#           9: 冥王星          (Pluto)
#          10: 月              (Moon)
#          11: 太陽            (Sun)
#          12: 太陽系重心      (Solar system Barycenter)
#          13: 地球 - 月の重心 (Earth-Moon Barycenter)
#          14: 地球の章動      (Earth Nutations)
#          15: 月の秤動        (Lunar mantle Librations)
#   [第２] 基準天体番号（必須。 0, 1 - 13）
#          （ 0 は、対象天体番号が 14, 15 のときのみ）
#   [第３] ユリウス日（省略可。省略時は現在日時のユリウス日を地球時とみなす）
#---------------------------------------------------------------------------------
# 【注意事項】
#   * 求める座標は「赤道直角座標(ICRS)」
#   * 時刻系は「太陽系力学時(TDB)」である。（≒地球時(TT)）
#   * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、
#     天体番号が 14 の場合は、黄経における章動 Δψ, 黄道傾斜における章動 Δε の
#     角位置・角速度（4要素）、
#     天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。
#   * 対象天体番号 = 基準天体番号 は、無意味なので処理しない。
#---------------------------------------------------------------------------------
#++
require 'eph_jpl'

class EphemerisJpl
  FILE_BIN = "/home/masaru/src/ephemeris_jpl/JPLEPH"
  USAGE = <<-"EOS"
【使用方法】 ./ephemeris_jpl.rb 対象天体番号 基準天体番号 [ユリウス日]
【天体番号】（対象天体番号: 1 - 15, 基準天体番号: 0 - 13）
   1: 水星,  2: 金星,  3: 地球,  4: 火星,  5: 木星,
   6: 土星,  7: 天王星,  8: 海王星,  9: 冥王星, 10: 月,
  11: 太陽,  12: 太陽系重心,  13: 地球 - 月の重心,
  14: 地球の章動,  15: 月の秤動,
   0: 対象天体番号 14, 15 の場合に指定
  ※対象天体番号≠基準天体番号
EOS
  KM = false  # true: km, false: AU

  def initialize
    target, center, jd = get_args
    @eph = EphJpl.new(FILE_BIN, target, center, jd, KM)
  end

  def exec
    display(@eph.calc)
  rescue => e
    msg = "[#{e.class}] #{e.message}\n"
    msg << e.backtrace.each { |tr| "\t#{tr}"}.join("\n")
    $stderr.puts msg
    exit 1
  end

  private

  #=========================================================================
  # コマンドライン引数取得
  #
  # @param:  <none>
  # @return: [target, center, jd]
  #=========================================================================
  def get_args
    # 対象、基準天体番号
    unless ARGV[0] && ARGV[1]
      puts USAGE
      exit 0
    end
    target = ARGV[0].to_i
    center = ARGV[1].to_i
    if (target < 1 || target > 15) ||
       (center < 0 || center > 13) ||
        target == center ||
       (target == 14 && center != 0) ||
       (target == 15 && center != 0) ||
       (target != 14 && target != 15 && center == 0)
      puts USAGE
      exit 0
    end

    # ユリウス日
    jd = ARGV[2]
    if jd
      jd = jd.to_f
    else
      now = Time.now
      jd = gc2jd({
        year: now.year, month: now.month, day: now.day,
        hour: now.hour, min: now.min, sec: now.sec
      })
    end
    return [target, center, jd]
  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 を超えない最大の整数
  # * 「ユリウス日」でなく「準ユリウス日」を求めるなら、
  #   `+ 1721088` を `- 678912` とする。
  #
  # @params: hash = {
  #            year(int), month(int), day(int), hour(int), min(int), sec(int)
  #          }
  # @return: jd (ユリウス日)
  #=========================================================================
  def gc2jd(hash)
    year  = hash[:year]
    month = hash[:month]
    day   = hash[:day]
    hour  = hash[:hour]
    min   = hash[:min]
    sec   = hash[:sec]

    begin
      # 1月,2月は前年の13月,14月とする
      if month < 3
        year  -= 1
        month += 12
      end
      # 日付(整数)部分計算
      jd  = (365.25 * year).floor
      jd += (year / 400.0).floor
      jd -= (year / 100.0).floor
      jd += (30.59 * (month - 2)).floor
      jd += day
      jd += 1721088.5
      # 時間(小数)部分計算
      t  = sec / 3600.0
      t += min / 60.0
      t += hour
      t /= 24.0
      return jd + t
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 結果出力
  #
  # @param:  pvs (= Position-Velocity array)
  # @return: <none>
  #=========================================================================
  def display(pvs)
    strs =  sprintf("  Target: %2d", @eph.target)
    strs << " (#{@eph.target_name})\n"
    strs << sprintf("  Center: %2d", @eph.center)
    strs << " (#{@eph.center_name})" unless @eph.center == 0
    strs << "\n"
    strs << sprintf("      JD: %16.8f day\n", @eph.jd)
    strs << sprintf("     1AU: %11.1f km\n", @eph.bin[:au])
    case @eph.target
    when 14
      strs << sprintf("  Position(Δψ) = %30.20f rad\n",     pvs[0])
      strs << sprintf("  Position(Δε) = %30.20f rad\n",     pvs[1])
      strs << sprintf("  Velocity(Δψ) = %30.20f rad/day\n", pvs[2])
      strs << sprintf("  Velocity(Δε) = %30.20f rad/day\n", pvs[3])
    when 15
      strs << sprintf("  Position(φ) = %30.20f rad\n",     pvs[0])
      strs << sprintf("  Position(θ) = %30.20f rad\n",     pvs[1])
      strs << sprintf("  Position(ψ) = %30.20f rad\n",     pvs[2])
      strs << sprintf("  Velocity(φ) = %30.20f rad/day\n", pvs[3])
      strs << sprintf("  Velocity(θ) = %30.20f rad/day\n", pvs[4])
      strs << sprintf("  Velocity(ψ) = %30.20f rad/day\n", pvs[5])
    else
      unit_p, unit_v = @eph.unit.split(", ")
      strs << sprintf("  Position(x) = %30.20f #{unit_p}\n", pvs[0])
      strs << sprintf("  Position(y) = %30.20f #{unit_p}\n", pvs[1])
      strs << sprintf("  Position(z) = %30.20f #{unit_p}\n", pvs[2])
      strs << sprintf("  Velocity(x) = %30.20f #{unit_v}\n", pvs[3])
      strs << sprintf("  Velocity(y) = %30.20f #{unit_v}\n", pvs[4])
      strs << sprintf("  Velocity(z) = %30.20f #{unit_v}\n", pvs[5])
    end
    puts strs
  rescue => e
    raise
  end
end

EphemerisJpl.new.exec if __FILE__ == $0