mk-mode BLOG

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

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

Ruby - JPL 天文暦バイナリデータの読み込み!

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

こんばんは。

NASA の機関である JPL(Jet Propulsion Laboratory) が惑星探査用に編纂・発行している月・惑星の暦の最新版 DE430 のバイナリ形式のデータを Ruby で読み込んでみました。

0. 前提条件

  • Ruby 2.3.0-p0 での作業を想定。
  • 使用するバイナリ形式データは、テキスト形式データ “ascp1950.430”, “ascp2050.430” を「JPL 天文暦データのバイナリ化!」の方法でバイナリ化したもの。(ファイル名は “JPLEPH” に変更)
  • あらかじめ用意されているバイナリデータを使用する場合は KSIZE 値が組み込まれていないので、別途テキスト形式のヘッダファイルを読み込む必要がある。
  • バイナリ形式データの仕様については「JPL 天文暦バイナリデータの仕様!」を参照。

1. Ruby スクリプトの作成

(バイナリデータの読み込みには IO.binread メソッドを使用)

jpl_read_de430.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
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
#! /usr/local/bin/ruby
# coding: utf-8
#---------------------------------------------------------------------------------
#= JPLEPH(JPL の DE430 バイナリデータ)読み込み
#  : データの読み込みテストのため、対象日時・天体番号は定数で設定している
#       1: 水星 (Mercury)
#       2: 金星 (Venus)
#       3: 地球 - 月の重心 (Earth-Moon barycenter)
#       4: 火星 (Mars)
#       5: 木星 (Jupiter)
#       6: 土星 (Saturn)
#       7: 天王星 (Uranus)
#       8: 海王星 (Neptune)
#       9: 冥王星 (Pluto)
#      10: 月(地心) (Moon (geocentric))
#      11: 太陽 (Sun)
#      12: 地球の章動(1980年IAUモデル) (Earth Nutations in longitude and obliquity (IAU 1980 model))
#      13: 月の秤動 (Lunar mantle libration)
#
# date          name            version
# 2016.03.11    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2016 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
#++
class JplReadDe430
  FILE_BIN = "JPLEPH"
  KSIZE    = 2036
  RECL     = 4
  JD       = 2457459.5  # <= 対象日時(2457459.5 = 2016-03-12 00:00:00)
  ASTR     = 1          # <= 天体番号

  def initialize
    @pos = 0      # レコード位置
  end

  def exec
    begin
      # ヘッダ(1レコード目)取得
      get_ttl     # TTL(タイトル)
      get_cnam    # CNAM(定数名)
      get_ss      # SS(ユリウス日(開始、終了)、分割日数)
      get_ncon    # NCON(定数の数)
      get_au      # AU(天文単位)
      get_emrat   # EMRAT(地球と月の質量比)
      get_ipt     # IPT(オフセット、係数の数、サブ区間数)(水星〜月の章動)
      get_numde   # NUMDE(DEバージョン番号)取得
      get_ipt_13  # IPT(オフセット、係数の数、サブ区間数)(月の秤動)
      p @ttl
      p @cnams
      p @sss
      p @ncon
      p @au
      p @emrat
      p @numde
      p @ipts

      # ヘッダ(2レコード目)取得
      get_cval    # CVAL(定数値)
      p @cvals

      # レコードインデックス取得
      @idx = ((JD - @sss[0]) / @sss[2]).truncate
      p @idx

      # 係数取得(対象のインデックス分を取得)
      get_coeff   # COEFF(係数)
      p @jds
      p @coeffs
    rescue => e
      $stderr.puts "[#{e.class}] #{e.message}"
      e.backtrace.each { |tr| $stderr.puts "\t#{tr}"}
      exit 1
    end
  end

  private

  # TTL 取得
  # - 84 byte * 3
  # - ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
  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 取得
  # - 6 byte * 400
  # - ASCII文字列(スペースを詰める/後続するnull文字やスペースを削除)
  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 取得
  # - 8 byte * 3
  # - 倍精度浮動小数点数(機種依存)
  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
  # - 4 byte * 1
  # - unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  def get_ncon
    recl = 4

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

  # AU 取得
  # - 8 byte * 1
  # - 倍精度浮動小数点数(機種依存)
  def get_au
    recl = 8

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

  # EMRAT 取得
  # - 8 byte * 1
  # - 倍精度浮動小数点数(機種依存)
  def get_emrat
    recl = 8

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

  # NUMDE
  # - 4 byte * 1
  # - unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  def get_numde
    recl = 4

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

  # IPT
  # - 4 byte * 12 * 3
  # - unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  def get_ipt
    recl = 4

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

  # IPT_13(月の秤動)
  # - 4 byte * 1 * 3
  # - unsigned int (符号なし整数, エンディアンと int のサイズに依存)
  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 取得
  # - 8 byte * @ncon
  # - 倍精度浮動小数点数(機種依存)
  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

  # CEFF 取得
  # - 8 byte * ?
  # - 倍精度浮動小数点数(機種依存)
  # - 地球・章動のみ要素数が 2 で、その他の要素数は 3
  # - @jds に当インデックスの開始・終了ユリウス日を格納
  # - @coeffs にその他の係数を格納
  def get_coeff
    @pos = KSIZE * RECL * (2 + @idx)
    recl = 8
    offset, cnt_coeff, cnt_sub = @ipts[ASTR - 1]
    n = ASTR == 12 ? 2 : 3

    begin
      ary = (0..(KSIZE / 2) - 1).map do |i|
        File.binread(FILE_BIN, recl, @pos + recl * i).unpack("d*")[0]
      end  # 該当インデックス分を全て取得
      @jds = ary[0, 2]
      ary = ary[offset - 1, cnt_coeff * n * cnt_sub]  # 該当惑星分のみ抽出
      @coeffs = ary.each_slice(cnt_coeff * n).map do |row|
        row.each_slice(cnt_coeff).to_a
      end
    rescue => e
      raise
    end
  end
end

JplReadDe430.new.exec

2. Ruby スクリプトの実行

バイナリ形式データ “JPLEPH” を Ruby スクリプトと同じディレクトリに配置して以下を実行する。

1
2
3
4
5
6
7
8
9
10
11
12
13
$ ./jpl_read_de430.rb
"JPL Planetary Ephemeris DE430/LE430\n
 Start Epoch: JED=  2287184.5 1549-DEC-21 00:00:00\n
 Final Epoch: JED=  2688976.5 2650-JAN-25 00:00:00"
["DENUM", "LENUM", "TDATEF", "TDATEB",

===< 中略 >===

[19870922.88638703, 3401048.340544464, -820016.1365685505,
 -21884.71452138689, 4458.980962630561, 113.7942416160033,
 -28.94032658038136, -0.7360655861367853, 0.2262947806520317,
 0.0053757421530248, -0.001954966141811797, -4.175805388132797e-05,
 1.792660121594814e-05, 3.35376232236008e-07]]]]

取得した値が出力される。(上記は視認性向上のため改行している)

必要であれば、出力を調整してご確認ください。


バイナリデータから目的の惑星・日付の係数が取得できるようになったので、チェビシェフの多項式を利用して座標(ICRS)を計算できることになります。

その後、赤道座標・黄道座標を計算できれば、こよみを正確に計算するために必要な太陽・月の黄経・黄緯等も算出できるはずです。(挑戦してみるつもりです)

以上。

Comments