mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Ruby - JPL 天文暦 gem の作成!

[ プログラミング, 暦・カレンダー ] [ Ruby, RubyGems, カレンダー ]

こんばんは。

以前、「Ruby - JPL 天文暦データから ICRS 座標を計算!」について紹介しました。

その際に使用した Ruby スクリプトを改変して gem ライブラリ化しました。
対象となる天体の番号・中心となる天体の番号・ユリウス日を指定すると、そのユリウス日の中心天体から見た対象天体の位置(直交座標)と速度を計算して返します。(座標系は ICRS(国際地球基準座標系))

以下では、今回作成した gem の簡単な利用方法をご紹介します。

0. 前提条件

  • Ruby 2.3.1-p112 での作業を想定。
  • 自作した gem ライブラリの名称は eph_jpl で、使用するバイナリデータは DE430.
  • JPL 天文暦とは何か?バイナリデータとは何か?等々は、過去記事「Ruby - JPL 天文暦データから ICRS 座標を計算!」等をご参照ください。

1. インストール

1
$ sudo gem install eph_jpl

2. Ruby スクリプトの作成例

ephemeris_jpl.rb
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
#! /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] 対象天体番号(必須)
#           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)
#   [第2] 基準天体番号(必須。 0, 1 - 13)
#          ( 0 は、対象天体番号が 14, 15 のときのみ)
#   [第3] ユリウス日(省略可。省略時は現在日時のユリウス日を地球時とみなす)
#---------------------------------------------------------------------------------
# 【注意事項】
#   * 求める座標は「赤道直角座標(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

exit 0 unless __FILE__ == $0
EphemerisJpl.new.exec

3. サンプルスクリプトの実行

[ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/] からバイナリファイル “linux_p1550p2650.430” を取得後に適当なディレクトリに配置し、ファイル名を上記のスクリプト内 FILE_BIN と合わせる。(上記のスクリプトでは JPLEPH を想定)
(このバイナリファイルはテキストファイルから生成することも可能。過去記事「JPL 天文暦データのバイナリ化!」を参照)

1
2
3
4
5
6
7
8
9
10
11
$ ./ephemeris_jpl.rb 11 3 2457558.5
  Target: 11 (Sun)
  Center:  3 (Earth)
      JD: 2457558.50000000 day
     1AU: 149597870.7 km
  Position(x) =         0.03679229733787987150 AU
  Position(y) =         0.93165937859412728539 AU
  Position(z) =         0.40388562663893279314 AU
  Velocity(x) =        -0.01690604098444144221 AU/day
  Velocity(y) =         0.00062873365728935307 AU/day
  Velocity(z) =         0.00027267631528302664 AU/day

計算結果の検証は、[ftp://ssd.jpl.nasa.gov/pub/eph/planets/Linux/de430/testpo.430] の内容と比較して行う。(gem ライブラリ側で検証済みなので、問題ないはず)

4. gem ライブラリ


以前作成した Ruby スクリプトを何かと応用(他の座標系への変換等)したかったので gem ライブラリ化した次第です。

以上。

Comments