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

Updated:


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

今回は、 Python で同様の実装をしてみました。(Python でのバイナリデータの読み込み方法に関する記録と考えてもよい)

0. 前提条件

  • LMDE 2 (Linux Mint Debian Edition 2; 64bit) での作業を想定。
  • Python 3.6.4 (エンコード:UTF-8)での作業を想定。
  • 使用するバイナリ形式データは、テキスト形式データ “ascp1950.430”, “ascp2050.430” を「JPL 天文暦データのバイナリ化!」の方法でバイナリ化したもの。(ファイル名は “JPLEPH” に変更)
  • バイナリ形式データの仕様については「JPL 天文暦バイナリデータの仕様!」を参照。

1. Python スクリプトの作成

File: jpl_read_de430.py

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
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
#! /usr/local/bin/python3.6
"""
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
  2018.03.18    mk-mode.com     1.00 新規作成

Copyright(C) 2018 mk-mode.com All Rights Reserved.
"""
import struct
import sys
import traceback


class JplReadDe430:
    FILE_BIN = "JPLEPH"
    KSIZE    = 2036
    RECL     = 4
    JD       = 2457459.5  # 対象日時(= 2016-03-12 00:00:00)
    ASTR     = 1          # 天体番号

    def __init__(self):
        self.pos = 0      # レコード位置

    def exec(self):
        """ Execution """
        try:
            # バイナリファイル読み込み
            with open(self.FILE_BIN, "rb") as f:
                # ヘッダ(1レコード目)取得
                self.__get_ttl(f)     # TTL
                self.__get_cnam(f)    # CNAM
                self.__get_ss(f)      # SS
                self.__get_ncon(f)    # NCON
                self.__get_au(f)      # AU
                self.__get_emrat(f)   # EMRAT
                self.__get_ipt(f)     # IPT
                self.__get_numde(f)   # NUMDE
                self.__get_ipt_13(f)  # IPT
                # ヘッダ(2レコード目)取得
                self.__get_cval(f)    # CVAL
                # レコードインデックス計算
                self.idx = int((self.JD - self.sss[0]) / self.sss[2])
                # 係数取得(対象のインデックス分を取得)
                self.__get_coeff(f)   # COEFF
            # 結果出力
            self.__display()
        except Exception as e:
            raise

    def __get_ttl(self, f):
        """ TTL(タイトル)取得
            - 84 byte * 3
            - ASCII文字列(後続のnull文字やスペースを削除)

        :param file_object f
        """
        len_rec = 84
        self.ttl = ""
        try:
            for i in range(3):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack(str(len_rec) + "s", f.read(len_rec))[0]
                if i != 0:
                    self.ttl += "\n"
                self.ttl += a.decode("utf-8").rstrip()
            self.pos += len_rec * 3
        except Exception as e:
            raise

    def __get_cnam(self, f):
        """ CNAM(定数名)取得
            - 6 byte * 400
            - ASCII文字列(後続のnull文字やスペースを削除)

        :param file_object f
        """
        len_rec = 6
        self.cnams = []
        try:
            for i in range(400):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack(str(len_rec) + "s", f.read(len_rec))[0]
                self.cnams.append(a.decode("utf-8").rstrip())
            self.pos += len_rec * 400
        except Exception as e:
            raise

    def __get_ss(self, f):
        """ SS(ユリウス日(開始、終了)、分割日数)取得
            - 8 byte * 3
            - 倍精度浮動小数点数(機種依存)

        :param file_object f
        """
        len_rec = 8
        self.sss = []
        try:
            for i in range(3):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack("d", f.read(len_rec))[0]
                self.sss.append(a)
            self.pos += len_rec * 3
        except Exception as e:
            raise

    def __get_ncon(self, f):
        """ NCON(定数の数)取得
            - 4 byte * 1
            - unsigned int (符号なし整数, エンディアンと int のサイズに依存)

        :param file_object f
        """
        len_rec = 4
        try:
            f.seek(self.pos)
            a = struct.unpack("I", f.read(len_rec))[0]
            self.ncon = a
            self.pos += len_rec
        except Exception as e:
            raise

    def __get_au(self, f):
        """ AU(天文単位)取得
            - 8 byte * 1
            - 倍精度浮動小数点数(機種依存)

        :param file_object f
        """
        len_rec = 8
        try:
            f.seek(self.pos)
            a = struct.unpack("d", f.read(len_rec))[0]
            self.au = a
            self.pos += len_rec
        except Exception as e:
            raise

    def __get_emrat(self, f):
        """ EMRAT(地球と月の質量比)取得
            - 8 byte * 1
            - 倍精度浮動小数点数(機種依存)

        :param file_object f
        """
        len_rec = 8
        try:
            f.seek(self.pos)
            a = struct.unpack("d", f.read(len_rec))[0]
            self.emrat = a
            self.pos += len_rec
        except Exception as e:
            raise

    def __get_numde(self, f):
        """ NUMDE(DEバージョン番号)取得
            - 4 byte * 1
            - unsigned int (符号なし整数, エンディアンと int のサイズに依存)

        :param file_object f
        """
        len_rec = 4
        try:
            f.seek(self.pos)
            a = struct.unpack("I", f.read(len_rec))[0]
            self.numde = a
            self.pos += len_rec
        except Exception as e:
            raise

    def __get_ipt(self, f):
        """ IPT(オフセット、係数の数、サブ区間数)(水星〜月の章動)取得
            - 4 byte * 12 * 3
            - unsigned int (符号なし整数, エンディアンと int のサイズに依存)

        :param file_object f
        """
        len_rec = 4
        self.ipts = []
        try:
            for i in range(12):
                l = []
                for j in range(3):
                    f.seek(self.pos + len_rec * j)
                    a = struct.unpack("I", f.read(len_rec))[0]
                    l.append(a)
                self.ipts.append(l)
                self.pos += len_rec * 3
        except Exception as e:
            raise

    def __get_ipt_13(self, f):
        """ IPT_13(オフセット、係数の数、サブ区間数)(月の秤動)取得
            - 4 byte * 1 * 3
            - unsigned int (符号なし整数, エンディアンと int のサイズに依存)

        :param file_object f
        """
        len_rec = 4
        try:
            l = []
            for i in range(3):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack("I", f.read(len_rec))[0]
                l.append(a)
            self.ipts.append(l)
            self.pos += len_rec * 3
        except Exception as e:
            raise

    def __get_cval(self, f):
        """ CVAL(定数値)取得
            - 8 byte * @ncon
            - 倍精度浮動小数点数(機種依存)

        :param file_object f
        """
        self.pos = self.KSIZE * self.RECL
        len_rec = 8
        self.cvals = []
        try:
            for i in range(self.ncon):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack("d", f.read(len_rec))[0]
                self.cvals.append(a)
        except Exception as e:
            raise

    def __get_coeff(self, f):
        """ CEFF(係数)取得
            - 8 byte * ?
            - 倍精度浮動小数点数(機種依存)
            - 地球・章動のみ要素数が 2 で、その他の要素数は 3
            - self.jds に当インデックスの開始・終了ユリウス日を格納
            - self.coeffs にその他の係数を格納

        :param file_object f
        """
        self.pos = self.KSIZE * self.RECL * (2 + self.idx)
        offset, cnt_coeff, cnt_sub = self.ipts[self.ASTR - 1]
        n = 2 if self.ASTR == 12 else 3
        len_rec = 8
        self.coeffs = []
        try:
            l = []
            for i in range(self.KSIZE // 2):
                f.seek(self.pos + len_rec * i)
                a = struct.unpack("d", f.read(len_rec))[0]
                l.append(a)
            self.jds = l[0:2]
            l = l[(offset - 1):(cnt_coeff * n * cnt_sub + offset - 1)]
            for a in [
                l[x:x + cnt_coeff * n]
                for x in range(0, len(l), cnt_coeff * n)
            ]:
                c = []
                for b in [
                    a[x:x + cnt_coeff]
                    for x in range(0, len(a), cnt_coeff)
                ]:
                    c.append(b)
                self.coeffs.append(c)
        except Exception as e:
            raise

    def __display(self):
        """ Display """
        try:
            print(self.ttl)
            print(self.cnams)
            print(self.sss)
            print(self.ncon)
            print(self.au)
            print(self.emrat)
            print(self.numde)
            print(self.ipts)
            print(self.cvals)
            print(self.idx)
            print(self.jds)
            print(self.coeffs)
        except Exception as e:
            raise


if __name__ == '__main__':
    try:
        obj = JplReadDe430()
        obj.exec()
    except Exception as e:
        traceback.print_exc()
        sys.exit(1)

2. Python スクリプトの実行

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

$ ./jpl_read_de430.py
JPL Planetary Ephemeris DE430/LE430
Start Epoch: JED=  2287184.5 1549-DEC-21 00:00:00
Final Epoch: JED=  2688976.5 2650-JAN-25 00:00:00
['DENUM', 'LENUM', 'TDATEF', 'TDATEB', 'JDEPOC', 
      :
===< 中略 >===
      :
[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]]]

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

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

3. 参考サイト

アンパックする際の書式文字列については以下のページを参照。


バイナリデータから目的の惑星・日付の係数が取得できるようになったので、天体の視位置計算に利用できることになります。

以上。





 

Sponsored Link

 

Comments