Ruby script to calculate planetary distances.

jpl_distance_earth.rb

#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読込後に座標（位置・速度）を計算し、さらに
#  地球から太陽・惑星・月までの距離を計算する。
#  : 自作の RubyGems ライブラリ eph_jpl を使用する
#    （JPL DE430 のバイナリデータが準備済みであること）
#
# date          name            version
# 2016.07.13    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# 【引数】 ユリウス日
#          （省略可。省略時は現在日時のユリウス日を地球時TT(≒太陽系力学時TDB)とみなす）
#---------------------------------------------------------------------------------
#++
require 'eph_jpl'

class JplDistanceEarth
  FILE_BIN = "/path/to/JPLEPH"
  USAGE    = "【使用方法】 ./jpl_distance_earth.rb [ユリウス日]"
  CENTER   = 3
  KM       = false  # true: km, false: AU
  BODIES = {
     1 => "  水星 (Mercury)",  2 => "  金星 (Venus)  ",  3 => "  地球 (Earth)  ",
     4 => "  火星 (Mars)   ",  5 => "  木星 (Jupiter)",  6 => "  土星 (Saturn) ",
     7 => "天王星 (Uranus) ",  8 => "海王星 (Neptune)",  9 => "冥王星 (Pluto)  ",
    10 => "    月 (Moon)   ", 11 => "  太陽 (Sun)    "
  }

  def initialize
    @jd   = get_args
    @unit = KM ? "KM" : "AU"
  end

  def exec
    puts "JD(TT) = #{@jd}\n#{BODIES[3]}"
    BODIES.each do |target, target_name|
      next if target == 3
      eph = EphJpl.new(FILE_BIN, target, CENTER, @jd, KM)
      ps = eph.calc[0,3]
      d = calc_distance(ps)
      d = KM ? sprintf("%19.8f", d) : sprintf("%11.8f", d)
      puts "- #{target_name} : #{d} #{@unit}"
    end
  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: jd
  #=========================================================================
  def get_args
    jd = ARGV.shift
    return jd ? jd.to_f : gc2jd({
      year:  Time.now.strftime("%Y").to_i,
      month: Time.now.strftime("%m").to_i,
      day:   Time.now.strftime("%d").to_i,
      hour:  Time.now.strftime("%H").to_i,
      min:   Time.now.strftime("%M").to_i,
      sec:   Time.now.strftime("%S").to_i
    })
  rescue => e
    puts USAGE
    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` とする。
  #
  # @param:  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 \
         + (year / 400.0).floor \
         - (year / 100.0).floor \
         + (30.59 * (month - 2)).floor \
         + day \
         + 1721088
      # 時間(小数)部分計算
      t = (sec / 3600.0 + min / 60.0 + hour) / 24.0
      return jd + t
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 距離計算
  #
  # @param:  ps       (x-座標, y-座標, z-座標)
  # @return: distance (距離)
  #=========================================================================
  def calc_distance(ps)
    return Math.sqrt(ps.inject(0) { |s, a| s + a * a })
  rescue => e
    raise
  end
end

JplDistanceEarth.new.exec if __FILE__ == $0