Ruby - JPL 天文暦データから ICRS 座標を計算！

#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読み込み、座標（位置・速度）を計算する
#
# date          name            version
# 2016.03.12    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)」
#   * 天体番号は、係数データの番号（並び順）と若干異なるので注意する。
#     （特に、天体番号 3, 10, と 12 以降）
#     係数データの並び順：
#       水星(1)、金星(2)、地球 - 月の重心(3)、火星(4)、木星(5)、土星(6)、
#       天王星(7)、海王星(8)、冥王星(9)、月（地心）(10)、太陽(11)、
#       地球の章動(12)、月の秤動(13)
#   * 時刻系は「太陽系力学時(TDB)」である。（≒地球時(TT)）
#   * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、
#     天体番号が 14 の場合は、黄経における章動 Δψ, 黄道傾斜における章動 Δε の
#     角位置・角速度（4要素）、
#     天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。
#   * 対象天体番号 = 基準天体番号 は、無意味なので処理しない。
#   * 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。
#   * その他、JPL 提供の FORTRAN プログラム "testeph.f" を参考にした。
#---------------------------------------------------------------------------------
#++
class JplCalcDe430
  USAGE = <<-"EOS"
【使用方法】 ./jpl_calc_de430.rb 対象天体番号 基準天体番号 [ユリウス日]
【天体番号】（対象天体番号: 1 - 15, 基準天体番号: 0 - 13）
   1: 水星,  2: 金星,  3: 地球,  4: 火星,  5: 木星,
   6: 土星,  7: 天王星,  8: 海王星,  9: 冥王星, 10: 月,
  11: 太陽,  12: 太陽系重心,  13: 地球 - 月の重心,
  14: 地球の章動,  15: 月の秤動,
   0: 対象天体番号 14, 15 の場合に指定
  ※対象天体番号≠基準天体番号
EOS
  FILE_BIN = "JPLEPH"
  KSIZE = 2036
  RECL  = 4
  ASTRS = [
    "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus",
    "Neptune", "Pluto", "Moon", "Sun", "Solar system Barycenter",
    "Earth-Moon barycenter", "Earth Nutations", "Lunar mantle Librations"
  ]
  KIND = 2      # 計算区分（0: 計算しない、1: 位置のみ計算、2: 位置・速度を計算）
                # （但し、現時点では 1 の場合の処理は実装していない）
  BARY = true   # 基準フラグ（true: 太陽系重心が基準 false: 太陽が基準）
  KM   = false  # 単位フラグ（true: km, km/sec, false: AU, AU/day）

  def initialize
    get_args                       # コマンドライン引数取得
    @pos      = 0                  # レコード位置
    @list     = Array.new(12, 0)   # 計算対象フラグ一覧配列
    @ipts     = Array.new          # 定数一覧配列
    @coeffs   = Array.new          # 係数一覧配列
    @idx_subs = Array.new(12)      # サブ区間のインデックス
    @tcs      = Array.new(12)      # サブ区間内の時間（チェビシェフ多項式用に正規化）
    @pvs      = Array.new(11).map { |a| Array.new(6, 0.0) }  # 位置・角度データ配列
    @pvs_2    = Array.new(13)      # 位置・角度データ配列（対象天体と基準天体の差算出用）
    @rrds     = Array.new(6, 0.0)  # 算出データ（対象天体と基準天体の差）配列
  end

  def exec
    begin
      # === ヘッダ情報取得
      get_ttl     # TTL（タイトル）
      get_cnam    # CNAM（定数名）
      get_ss      # SS（ユリウス日（開始、終了）、分割日数）
      check_jd    # 引数のユリウス日をチェック
      get_ncon    # NCON（定数の数）
      get_au      # AU（天文単位）
      get_emrat   # EMRAT（地球と月の質量比）
      get_ipt     # IPT（オフセット、係数の数、サブ区間数）（水星〜地球の章動）
      get_numde   # NUMDE（DEバージョン番号）取得
      get_ipt_13  # IPT（オフセット、係数の数、サブ区間数）（月の秤動）
      get_cval    # CVAL（定数値）
      get_jdepoc  # JDEPOC（元期）

      # === 各種計算
      get_list    # 計算対象フラグ一覧（係数データの並びに対応）取得
      get_coeff   # 係数取得（対象、基準天体のインデックス分を取得）
      #p @ttl
      #p @cnams
      #p @sss
      #p @ncon
      #p @au
      #p @emrat
      #p @ver_de
      #p @ipts
      #p @cvals
      #p @jdepoc
      #p @list
      #p @coeffs
      #exit

      # 補間（11:太陽）
      @pv_sun = interpolate(11)
      # 補間（1:水星〜10:月）
      0.upto(9) do |i|
        next if @list[i] == 0
        @pvs[i] = interpolate(i + 1)
        next if i > 8
        next if BARY
        @pvs[i] = @pvs[i].map.with_index do |pv, j|
          pv - @pv_sun[j]
        end
      end
      # 補間（14:地球の章動）
      if @list[10] > 0 && @ipts[11][1] > 0
        @p_nut = interpolate(14)
      end
      # 補間（15:月の秤動）
      if @list[11] > 0 && @ipts[12][1] > 0
        @pvs[10] = interpolate(15)
      end

      # 対象天体と基準天体の差
      case
      when @astrs[0] == 14
        @rrds = @p_nut + [0.0, 0.0] if @ipts[11][1] > 0
      when @astrs[0] == 15
        @rrds = @pvs[10] if @ipts[12][1] > 0
      else
        0.upto(9) { |i| @pvs_2[i] = @pvs[i] }
        @pvs_2[10] = @pv_sun           if @astrs.include?(11)
        @pvs_2[11] = Array.new(6, 0.0) if @astrs.include?(12)
        @pvs_2[12] = @pvs[2]           if @astrs.include?(13)
        if @astrs[0] * @astrs[1] == 30 || @astrs[0] + @astrs[1] == 13
          @pvs_2[2] = Array.new(6, 0.0)
        else
          @pvs_2[2] = @pvs[2].map.with_index do |pv, i|
            pv - @pvs[9][i] / (1.0 + @emrat)
          end unless @list[2] == 0
          @pvs_2[9] = @pvs_2[2].map.with_index do |pv, i|
            pv + @pvs[9][i]
          end unless @list[9] == 0
        end
        0.upto(5) do |i|
          @rrds[i] = @pvs_2[@astrs[0] - 1][i] - @pvs_2[@astrs[1] - 1][i]
        end
      end

      # === 結果出力
      display
    rescue => e
      $stderr.puts "[#{e.class}] #{e.message}"
      e.backtrace.each { |tr| $stderr.puts "\t#{tr}"}
      exit 1
    end
  end

  private

  #=========================================================================
  # コマンドライン引数取得
  #
  #   * 第１引数の値をインスタンス変数 @astrs[0] に設定する
  #   * 第２引数の値をインスタンス変数 @astrs[1] に設定する
  #   * 第３引数の値をインスタンス変数 @jd に設定する
  #     （但し、第３引数が存在しない場合は現在日時を設定する）
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  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
    @astrs = [target, center]

    # ユリウス日
    @jd = ARGV[2]
    @jd = @jd ? @jd.to_f : gc_to_jd({
      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
    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 gc_to_jd(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

      # 時間(小数)部分計算
      t  = sec / 3600.0
      t += min / 60.0
      t += hour
      t /= 24.0

      return jd + t
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 引数のユリウス日をチェック
  #
  #   * @jd の値をヘッダ情報のユリウス日（開始、終了）と比較・チェックする
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def check_jd
    if @jd < @sss[0] || @jd >= @sss[1]
      puts "Please input JD s.t. #{@sss[0]} <= JD < #{@sss[1]}."
      exit 0
    end
  rescue => e
    raise
  end

  #=========================================================================
  # TTL 取得
  #
  #   * タイトルをインスタンス変数 @ttl に設定する
  #   * 84 byte * 3
  #   * ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_ttl
    recl = 84

    begin
      @ttl = (0..2).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("A*")[0]
      end.join("\n")
      @pos += recl * 3
    rescue => e
      raise
    end
  end

  #=========================================================================
  # CNAM 取得
  #
  #   * 定数名をインスタンス変数（配列） @cnams に設定する
  #   * 6 byte * 400
  #   * ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_cnam
    recl = 6

    begin
      @cnams = (0..399).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("A*")[0]
      end
      @pos += recl * 400
    rescue => e
      raise
    end
  end

  #=========================================================================
  # SS 取得
  #
  #   * ユリウス日（開始、終了）、分割日数をインスタンス変数（配列） @sss に
  #     設定する
  #   * 8 byte * 3
  #   * 倍精度浮動小数点数(機種依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_ss
    recl = 8

    begin
      @sss = (0..2).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
      end
      @pos += recl * 3
    rescue => e
      raise
    end
  end

  #=========================================================================
  # NCON 取得
  #
  #   * 定数の数をインスタンス変数 @ncon に設定する
  #   * 4 byte * 1
  #   * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_ncon
    recl = 4

    begin
      @ncon = File.binread(FILE_BIN, recl, @pos).unpack("I*")[0]
      @pos += recl
    rescue => e
      raise
    end
  end

  #=========================================================================
  # AU 取得
  #
  #   * AU（天文単位）の値をインスタンス変数 @au に設定する
  #   * 8 byte * 1
  #   * 倍精度浮動小数点数(機種依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_au
    recl = 8

    begin
      @au = File.binread(FILE_BIN, recl, @pos).unpack("d*")[0]
      @pos += recl
    rescue => e
      raise
    end
  end

  #=========================================================================
  # EMRAT 取得
  #
  #   * EMRAT（地球と月の質量比）の値をインスタンス変数 @emrat に設定する
  #   * 8 byte * 1
  #   * 倍精度浮動小数点数(機種依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_emrat
    recl = 8

    begin
      @emrat = File.binread(FILE_BIN, recl, @pos).unpack("d*")[0]
      @pos += recl
    rescue => e
      raise
    end
  end

  #=========================================================================
  # NUMDE 取得
  #
  #   * DEバージョン番号(NUMDE)をインスタンス変数 @ver_de に設定する
  #   * 4 byte * 1
  #   * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_numde
    recl = 4

    begin
      @ver_de = File.binread(FILE_BIN, recl, @pos).unpack("I*")[0]
      @pos += recl
    rescue => e
      raise
    end
  end

  #=========================================================================
  # IPT 取得
  #
  #   * IPT（オフセット、係数の数、サブ区間数）（水星〜地球の章動）の値をイ
  #     ンスタンス変数（配列） @ipts に設定する
  #   * 4 byte * 12 * 3
  #   * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_ipt
    recl = 4

    begin
      @ipts = (0..11).map do |i|
        ary = (0..2).map do |j|
          File.binread(FILE_BIN, recl, @pos + recl * j).unpack("I*")[0]
        end
        @pos += recl * 3
        ary
      end
    rescue => e
      raise
    end
  end

  #=========================================================================
  # IPT_13（月の秤動）
  #
  #   * IPT（オフセット、係数の数、サブ区間数）（月の秤動）の値をインスタン
  #     ス変数（配列） @ipts に設定する
  #   * 4 byte * 1 * 3
  #   * unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  #   * 取得したバイト数だけレコード位置 @pos を進める
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_ipt_13
    recl = 4

    begin
      @ipts << (0..2).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("I*")[0]
      end
      @pos += recl * 3
    rescue => e
      raise
    end
  end

  #=========================================================================
  # CVAL 取得
  #
  #   * レコード位置 @pos を2レコード目の先頭へ進める
  #   * 定数の値をインスタンス変数（配列） @cvals に設定する
  #   * 8 byte * @ncon
  #   * 倍精度浮動小数点数(機種依存)
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_cval
    @pos = KSIZE * RECL
    recl = 8

    begin
      @cvals = (0..@ncon - 1).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
      end
    rescue => e
      raise
    end
  end

  #=========================================================================
  # JDEPOC 取得
  #
  #   * @cvals（定数値配列）の中から「元期」をインスタンス変数 @jdepoc に設
  #     定する
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_jdepoc
    @jdepoc = @cvals[4]
  rescue => e
    raise
  end

  #=========================================================================
  # 計算対象フラグ一覧取得
  #
  #   * チェビシェフ多項式による計算が必要な天体の一覧をインスタンス変数（配
  #     列） @list に設定する
  #   * 配列の並び順（係数データの並び順から「太陽」を除外した13個）
  #     @list = [
  #       水星, 金星, 地球 - 月の重心, 火星, 木星, 土星, 天王星, 海王星,
  #       冥王星, 月（地心）, 地球の章動, 月の秤動
  #     ]
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_list
    if @astrs[0] == 14
      @list[10] = KIND if @ipts[11][1] > 0
      return
    end
    if @astrs[0] == 15
      @list[11] = KIND if @ipts[12][1] > 0
      return
    end
    @astrs.each do |k|
      @list[k - 1] = KIND if k <= 10
      @list[2]     = KIND if k == 10
      @list[9]     = KIND if k ==  3
      @list[2]     = KIND if k == 13
    end
  rescue => e
    raise
  end

  #=========================================================================
  # COEFF 取得
  #
  #   * レコード位置 @pos を3レコード目の先頭へ進める
  #   * 対象区間のユリウス日（開始、終了）をインスタンス変数（配列） @jds に
  #     設定する
  #   * 係数の値をインスタンス変数（配列） @coeffs に設定する
  #     （x, y, z やサブ区間毎に分割して格納）
  #   * 8 byte * ?
  #   * 倍精度浮動小数点数(機種依存)
  #   * @coeffs に対象のインデックスの全ての係数を格納
  #   * 最初の2要素は当該データの開始・終了ユリウス日
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def get_coeff
    idx = ((@jd - @sss[0]) / @sss[2]).floor  # レコードインデックス
    @pos = KSIZE * RECL * (2 + idx)
    recl = 8

    begin
      items = (0..(KSIZE / 2) - 1).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
      end
      @jds = [items.shift, items.shift]
      @ipts.each_with_index do |ipt, i|
        n = i == 11 ? 2 : 3  # 要素数
        ary_1 = Array.new
        ipt[2].times do |j|
          ary_0 = Array.new
          n.times do |k|
            ary_0 << items.shift(ipt[1])
          end
          ary_1 << ary_0
        end
        @coeffs << ary_1
      end
    rescue => e
      raise
    end
  end

  #=========================================================================
  # 補間
  #
  #   * 使用するチェビシェフ多項式の係数は、
  #   * 天体番号が 1 〜 13 の場合は、 x, y, z の位置・速度（6要素）、
  #     天体番号が 14 の場合は、 Δψ, Δε の角位置・角速度（4要素）、
  #     天体番号が 15 の場合は、 φ, θ, ψ の角位置・角速度（6要素）。
  #   * 天体番号が 12 の場合は、 x, y, z の位置・速度の値は全て 0.0 とする。
  #
  #   @params: astr（天体番号）
  #   @return: pvs = [
  #     x 位置, y 位置, z 位置,
  #     x 速度, y 速度, z 速度
  #   ]
  #     - 14（地球の章動）の場合は、
  #         pvs = [
  #           Δψ の角位置, Δε の角位置,
  #           Δψ の角速度, Δε の角速度
  #         ]
  #     - 15（月の秤動）の場合は、
  #         pvs = [
  #           φ の角位置, θ の角位置, ψ の角位置,
  #           φ の角速度, θ の角速度, ψ の角速度
  #         ]
  #=========================================================================
  def interpolate(astr)
    pvs = Array.new

    begin
      tc, idx_sub = norm_time(astr)
      n_item = astr == 14 ? 2 : 3  # 要素数
      i_ipt  = astr > 13 ? astr - 3 : astr - 1
      i_coef = astr > 13 ? astr - 3 : astr - 1

      # 位置
      ary_p = [1, tc]
      2.upto(@ipts[i_ipt][1] - 1) do |i|
        ary_p << 2 * tc * ary_p[-1] - ary_p[-2]
      end  # 各項
      0.upto(n_item - 1) do |i|
        val = 0
        0.upto(@ipts[i_ipt][1] - 1) do |j|
          val += @coeffs[i_coef][idx_sub][i][j] * ary_p[j]
        end
        val /= @au if !KM && astr < 14
        pvs << val
      end  # 値

      # 速度
      ary_v = [0, 1, 2 * 2 * tc]
      3.upto(@ipts[i_ipt][1] - 1) do |i|
        ary_v << 2 * tc * ary_v[-1] + 2 * ary_p[i - 1] - ary_v[-2]
      end  # 各項
      0.upto(n_item - 1) do |i|
        val = 0
        0.upto(@ipts[i_ipt][1] - 1) do |j|
          val += @coeffs[i_coef][idx_sub][i][j] * ary_v[j] * 2 * @ipts[i_ipt][2] / @sss[2].to_f
        end
        val /= KM ? 86400.0 : @au if astr < 14
        pvs << val
      end  # 値

      return pvs
    rescue => e
      raise
    end
  end

  #=========================================================================
  # チェビシェフ多項式用に時刻を正規化、サブ区間のインデックス算出
  #
  #   @params: astr（天体番号）
  #   @return: [チェビシェフ時間, サブ区間のインデックス]
  #=========================================================================
  def norm_time(astr)
    idx = astr > 13 ? astr - 2 : astr
    jd_start = @jds[0]
    tc = (@jd - jd_start) / @sss[2].to_f
    temp = tc * @ipts[idx - 1][2]
    idx = (temp - tc.floor).floor         # サブ区間のインデックス
    tc = (temp % 1.0 + tc.floor) * 2 - 1  # チェビシェフ時間
    return [tc, idx]
  rescue => e
    raise
  end

  #=========================================================================
  # 結果出力
  #
  #   @params: <none>
  #   @return: <none>
  #=========================================================================
  def display
    strs =  "  Target: %2d" % @astrs[0]
    strs << " (#{ASTRS[@astrs[0] - 1]})\n"
    strs << "  Center: %2d" % @astrs[1]
    strs << " (#{ASTRS[@astrs[1] - 1]})" unless @astrs[1] == 0
    strs << "\n"
    strs << "      JD: %16.8f day\n" % @jd
    strs << "     1AU: %11.1f km\n" % @au
    puts strs

    strs = ""
    case @astrs[0]
    when 14
      strs << "  Position(Δψ) = %32.20f rad\n"     % @rrds[0]
      strs << "  Position(Δε) = %32.20f rad\n"     % @rrds[1]
      strs << "  Velocity(Δψ) = %32.20f rad/day\n" % @rrds[2]
      strs << "  Velocity(Δε) = %32.20f rad/day\n" % @rrds[3]
    when 15
      strs << "  Position(φ) = %32.20f rad\n"     % @rrds[0]
      strs << "  Position(θ) = %32.20f rad\n"     % @rrds[1]
      strs << "  Position(ψ) = %32.20f rad\n"     % @rrds[2]
      strs << "  Velocity(φ) = %32.20f rad/day\n" % @rrds[3]
      strs << "  Velocity(θ) = %32.20f rad/day\n" % @rrds[4]
      strs << "  Velocity(ψ) = %32.20f rad/day\n" % @rrds[5]
    else
      unit_0, unit_1 = KM ? ["km", "sec"] : ["AU", "day"]
      strs << "  Position(x) = %32.20f #{unit_0}\n"           % @rrds[0]
      strs << "  Position(y) = %32.20f #{unit_0}\n"           % @rrds[1]
      strs << "  Position(z) = %32.20f #{unit_0}\n"           % @rrds[2]
      strs << "  Velocity(x) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[3]
      strs << "  Velocity(y) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[4]
      strs << "  Velocity(z) = %32.20f #{unit_0}/#{unit_1}\n" % @rrds[5]
    end
    puts strs
  rescue => e
    raise
  end
end

JplCalcDe430.new.exec if __FILE__ == $0

$ ./jpl_calc_de430.rb 11 3 2457465.5
  Target: 11 (Sun)
  Center:  3 (Earth)
      JD: 2457465.50000000 day
     1AU: 149597870.7 km
  Position(x) =         0.99443659220700997281 AU
  Position(y) =        -0.03816291768957833647 AU
  Position(z) =        -0.01655177670960058384 AU
  Velocity(x) =         0.00099333756561238804 AU/day
  Velocity(y) =         0.01582779844821733881 AU/day
  Velocity(z) =         0.00686186627679565356 AU/day

de#  -- date -- -- jed -- t# c# x# -- coordinate ---

430  2016.12.01 2457723.5  8  5  2       -8.02225801549924900000

$ ./jpl_calc_de430.rb 8 5 2457723.5
  Target:  8 (Neptune)
  Center:  5 (Jupiter)
      JD: 2457723.50000000 day
     1AU: 149597870.7 km
  Position(x) =        33.70070119505614059108 AU
  Position(y) =        -8.02225801549925066070 AU
  Position(z) =        -4.10459769835983934172 AU
  Velocity(x) =        -0.00000061831278598047 AU/Day
  Velocity(y) =         0.00929863072586037051 AU/Day
  Velocity(z) =         0.00393151443263963212 AU/Day