mk-mode BLOG

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

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

Python - 章動の計算(IAU2000A 理論)!

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

こんばんは。

以前、天体位置の回転に使用する章動を Ruby で実装して計算してみました。(IAU2000A 理論を使用)

今回は、同じアルゴリズムを Python で実装してみました。

0. 前提条件、事前知識

  • Python 3.6.5 での作業を想定。
  • 章動の計算には、 IAU SOFA(International Astronomical Union, Standards of Fundamental Astronomy) の提供する C ソースコードに実装されているアルゴリズム “nut00a.c” を使用する。
  • IAU SOFA のソースコードには、 MHB2000(Mathews-Herring-Buffett, 2000) の理論や IERS2003(International Earth Rotation & Reference Systems service, 2003) の理論の使用が混在していることに留意。
  • ここでは「章動(しょうどう)」そのものが何かについては詳細には説明しないが、簡単に説明すると、章動には
    • 黄経における章動()
    • 黄道傾斜における章動()
      があり、それぞれが
    • 日月章動(luni-solar nutation)
    • 惑星章動(planetary nutation)
      で構成されている。
  • また、算出アルゴリズムについてもここでは詳細には説明しない。(と言うより、煩雑で自分には説明できない)
    参考サイトやソーススクリプトを参照のこと。
  • 歳差・章動の変換については、過去記事「歳差・章動の変換行列について!」も参照のこと。
  • 今回は IAU2000A 理論に基づいて計算しているが、簡略版の IAU2000B 理論 “nut00b.c” も存在する。
  • IAU2000A 理論に若干の補正を加えた IAU2006 理論も存在するが、今回は非考慮。

1. 係数データの準備

計算に使用する係数データを用意する。

IERS Conventions Center (2003)FTP archive ページから、日月章動用の係数データ Lunisolar components と惑星章動用の係数データ Planetary components を取得して、扱いやすいように整形する。

整形したものは以下にアップしている。

2. Ruby スクリプトの作成

IAU SOFA の提供する C ソースコードに実装されているアルゴリズム “nut00a.c” に則って作成している。

nutation_model.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
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
#! /usr/local/bin/python3.6
"""
章動の計算
: IAU2000A 章動理論(MHB2000, IERS2003)による
  黄経における章動(Δψ), 黄道傾斜における章動(Δε) の計算

* IAU SOFA(International Astronomical Union, Standards of Fundamental Astronomy)
  の提供する C ソースコード "nut00a.c" で実装されているアルゴリズムを使用する。
* 係数データファイルの項目について
  - 日月章動(luni-solar nutation, "dat_ls.txt")
    (左から) L L' F D Om PS PST PC EC ECT ES
  - 惑星章動(planetary nutation, "dat_pl.txt)
    (左から) L L' F D Om Lm Lv Le LM Lj Ls Lu Ln Pa PS PC ES EC
* 参考サイト
  - [SOFA Library Issue 2012-03-01 for ANSI C: Complete List](http://www.iausofa.org/2012_0301_C/sofa/)
  - [USNO Circular 179](http://aa.usno.navy.mil/publications/docs/Circular_179.php)
  - [IERS Conventions Center](http://62.161.69.131/iers/conv2003/conv2003_c5.html)

  Date          Author          Version
  2018.04.11    mk-mode.com     1.00 新規作成

Copyright(C) 2018 mk-mode.com All Rights Reserved.
---
引数 : 日時(TT(地球時))
         書式:YYYYMMDD or YYYYMMDDHHMMSS or YYYYMMDDHHMMSSUUUUUU
         無指定なら現在(システム日時)を地球時とみなす。
"""
from datetime import datetime
from datetime import date
import math
import re
import sys
import traceback

class NutationModel:
    DAT_LS = "nut_ls.txt"
    DAT_PL = "nut_pl.txt"
    PI     = 3.141592653589793238462643     # PI
    PI2    = 6.283185307179586476925287     # 2 * PI
    AS2R   = 4.848136811095359935899141e-6  # Arcseconds to radians
    TURNAS = 1296000.0                      # Arcseconds in a full circle
    U2R    = AS2R / 1e7                     # Units of 0.1 microarcsecond to radians
    R2D    = 57.29577951308232087679815     # Radians to degrees
    D2S    = 3600.0                         # Degrees to seconds

    def __init__(self):
        self.__get_arg()
        self.__get_data()

    def exec(self):
        """ 実行 """
        try:
            jd = self.__calc_jd(self.tt)
            t  = self.__calc_t(jd)
            dpsi_ls, deps_ls = self.__calc_lunisolar(t)
            dpsi_pl, deps_pl = self.__calc_planetary(t)
            self.dpsi = dpsi_ls + dpsi_pl
            self.deps = deps_ls + deps_pl
            self.dpsi_d = self.dpsi * self.R2D
            self.deps_d = self.deps * self.R2D
            self.dpsi_s = self.dpsi_d * self.D2S
            self.deps_s = self.deps_d * self.D2S
            self.__display()
        except Exception as e:
            raise

    def __get_arg(self):
        """ コマンドライン引数の取得
            * コマンドライン引数で指定した日時を self.tt に設定
            * コマンドライン引数が存在しなければ、現在時刻を self.tt に設定
        """
        try:
            if len(sys.argv) < 2:
                self.tt = datetime.now()
                return
            if re.search(r"^(\d{8}|\d{14}|\d{20})$", sys.argv[1]):
                dt = sys.argv[1].ljust(20, "0")
            else:
                sys.exit(0)
            try:
                self.tt = datetime.strptime(dt, "%Y%m%d%H%M%S%f")
            except ValueError as e:
                print("Invalid date!")
                sys.exit(0)
        except Exception as e:
            raise

    def __get_data(self):
        """ テキストファイル(DAT_LS, DAT_PL)からデータ取得
            * luni-solar の最初の5列、planetary の最初の14列は整数に、
              残りの列は浮動小数点*10000にする
            * 読み込みデータは self.dat_ls, self.dat_pl に格納
        """
        self.dat_ls = []
        self.dat_pl = []
        try:
            with open(self.DAT_LS, "r") as f:
                data = f.read()
                for l in re.split('\n', data)[1:]:
                    l = re.sub(r'^\s+', "", l)
                    items = re.split(r'\s+', l)
                    if len(items) < 2:
                        break
                    items = [int(x) for x in items[:5]] \
                          + [int(re.sub(r'\.', "", x)) for x in items[5:]]
                    self.dat_ls.append(items)
            with open(self.DAT_PL, "r") as f:
                data = f.read()
                for l in re.split('\n', data)[1:]:
                    l = re.sub(r'^\s+', "", l)
                    items = re.split(r'\s+', l)
                    if len(items) < 2:
                        break
                    items = [int(x) for x in items[:14]] \
                          + [int(re.sub(r'\.', "", x)) for x in items[14:]]
                    self.dat_pl.append(items)
        except Exception as e:
            raise

    def __calc_jd(self, tt):
        """ ユリウス日の計算
            * 地球時 self.tt のユリウス日を計算し、self.jd に設定

        :param  datetime tt: 地球時
        :return float      : ユリウス日
        """
        year, month,  day    = tt.year, tt.month,  tt.day
        hour, minute, second = tt.hour, tt.minute, tt.second
        try:
            if month < 3:
                year  -= 1
                month += 12
            d = int(365.25 * year) + year // 400  - year // 100 \
              + int(30.59 * (month - 2)) + day + 1721088.5
            t  = (second / 3600 + minute / 60 + hour) / 24
            return d + t
        except Exception as e:
            raise

    def __calc_t(self, jd):
        """ ユリウス世紀数の計算
            * ユリウス日 self.jd のユリウス世紀数を計算し、 self.t に設定

        :param  float jd: ユリウス日
        :return float   : ユリウス世紀数
        """
        try:
            return (jd - 2451545) / 36525
        except Exception as e:
            raise

    def __calc_lunisolar(self, t):
        """ 日月章動(luni-solar nutation)の計算

        :param  float t: ユリウス世紀数
        :return list: [黄経における章動(Δψ), 黄道傾斜における章動(Δε)]
        """
        dp, de = 0.0, 0.0
        try:
            l  = self.__l_iers2003(t)
            lp = self.__lp_mhb2000(t)
            f  = self.__f_iers2003(t)
            d  = self.__d_mhb2000(t)
            om = self.__om_iers2003(t)
            for x in reversed(self.dat_ls):
                arg = (x[0] * l + x[1] * lp + x[2] * f \
                     + x[3] * d + x[4] * om) % self.PI2
                sarg, carg = math.sin(arg), math.cos(arg)
                dp += (x[5] + x[6] * t) * sarg + x[ 7] * carg
                de += (x[8] + x[9] * t) * carg + x[10] * sarg
            return [dp * self.U2R, de * self.U2R]
        except Exception as e:
            raise

    def __calc_planetary(self, t):
        """ 惑星章動(planetary nutation)

        :param  float t: ユリウス世紀数
        :return list: [黄経における章動(Δψ), 黄道傾斜における章動(Δε)]
        """
        dp, de = 0.0, 0.0
        try:
            l  = self.__l_mhb2000(t)
            f  = self.__f_mhb2000(t)
            d  = self.__d_mhb2000_2(t)
            om = self.__om_mhb2000(t)
            pa = self.__pa_iers2003(t)
            me = self.__lme_iers2003(t)
            ve = self.__lve_iers2003(t)
            ea = self.__lea_iers2003(t)
            ma = self.__lma_iers2003(t)
            ju = self.__lju_iers2003(t)
            sa = self.__lsa_iers2003(t)
            ur = self.__lur_iers2003(t)
            ne = self.__lne_mhb2000(t)
            for x in reversed(self.dat_pl):
                arg = (x[ 0] * l  + x[ 2] * f  + x[ 3] * d  + x[ 4] * om \
                     + x[ 5] * me + x[ 6] * ve + x[ 7] * ea + x[ 8] * ma \
                     + x[ 9] * ju + x[10] * sa + x[11] * ur + x[12] * ne \
                     + x[13] * pa) % self.PI2
                sarg, carg = math.sin(arg), math.cos(arg)
                dp += x[14] * sarg + x[15] * carg
                de += x[16] * sarg + x[17] * carg
            return [dp * self.U2R, de * self.U2R]
        except Exception as e:
            raise

    def __l_iers2003(self, t):
        """ Mean anomaly of the Moon (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Mean anomaly of the Moon
        """
        try:
            return ((    485868.249036    \
                  + (1717915923.2178      \
                  + (        31.8792      \
                  + (         0.051635    \
                  + (        -0.00024470) \
                  * t) * t) * t) * t) % self.TURNAS) * self.AS2R
        except Exception as e:
            raise

    def __lp_mhb2000(self, t):
        """ Mean anomaly of the Sun (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean anomaly of the Sun
        """
        try:
            return ((  1287104.79305     \
                  + (129596581.0481      \
                  + (       -0.5532      \
                  + (        0.000136    \
                  + (       -0.00001149) \
                  * t) * t) * t) * t) % self.TURNAS) * self.AS2R
        except Exception as e:
            raise

    def __f_iers2003(self, t):
        """ Mean longitude of the Moon minus that of the ascending node
            (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Mean longitude of the Moon minus that of the
                         ascending node
        """
        try:
            return ((     335779.526232   \
                  + (1739527262.8478      \
                  + (       -12.7512      \
                  + (        -0.001037    \
                  + (         0.00000417) \
                  * t) * t) * t) * t) % self.TURNAS) * self.AS2R
        except Exception as e:
            raise

    def __d_mhb2000(self, t):
        """ Mean elongation of the Moon from the Sun (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean elongation of the Moon from the Sun
        """
        try:
            return ((   1072260.70369     \
                  + (1602961601.2090      \
                  + (        -6.3706      \
                  + (         0.006593    \
                  + (        -0.00003169) \
                  * t) * t) * t) * t) % self.TURNAS) * self.AS2R
        except Exception as e:
            raise

    def __om_iers2003(self, t):
        """ Mean longitude of the ascending node of the Moon
            (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Mean longitude of the ascending node of the
                         Moon
        """
        try:
            return ((    450160.398036    \
                  + (  -6962890.5431      \
                  + (         7.4722      \
                  + (         0.007702    \
                  + (        -0.00005939) \
                  * t) * t) * t) * t) % self.TURNAS) * self.AS2R
        except Exception as e:
            raise

    def __l_mhb2000(self, t):
        """ Mean anomaly of the Moon (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean anomaly of the Moon
        """
        try:
            return (2.35555598 + 8328.6914269554 * t) % self.PI2
        except Exception as e:
            raise

    def __f_mhb2000(self, t):
        """ Mean longitude of the Moon minus that of the ascending node
            (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean longitude of the Moon minus that of the
                         ascending node
        """
        try:
            return (1.627905234 + 8433.466158131 * t) % self.PI2
        except Exception as e:
            raise

    def __d_mhb2000_2(self, t):
        """ Mean elongation of the Moon from the Sun (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean elongation of the Moon from the Sun
        """
        try:
            return (5.198466741 + 7771.3771468121 * t) % self.PI2
        except Exception as e:
            raise

    def __om_mhb2000(self, t):
        """ Mean longitude of the ascending node of the Moon (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Mean longitude of the ascending node of the
                         Moon
        """
        try:
            return (2.18243920 - 33.757045 * t) % self.PI2
        except Exception as e:
            raise

    def __pa_iers2003(self, t):
        """ General accumulated precession in longitude (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : General accumulated precession in longitude
        """
        try:
            return (0.024381750 + 0.00000538691 * t) * t
        except Exception as e:
            raise

    def __lme_iers2003(self, t):
        """ Mercury longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Mercury longitudes
        """
        try:
            return (4.402608842 + 2608.7903141574 * t) % self.PI2
        except Exception as e:
            raise

    def __lve_iers2003(self, t):
        """ Venus longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Venus longitudes
        """
        try:
            return (3.176146697 + 1021.3285546211 * t) % self.PI2
        except Exception as e:
            raise

    def __lea_iers2003(self, t):
        """ Earth longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Earth longitudes
        """
        try:
            return (1.753470314 + 628.3075849991 * t) % self.PI2
        except Exception as e:
            raise

    def __lma_iers2003(self, t):
        """ Mars longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Mars longitudes
        """
        try:
            return (6.203480913 + 334.0612426700 * t) % self.PI2
        except Exception as e:
            raise

    def __lju_iers2003(self, t):
        """ Jupiter longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Jupiter longitudes
        """
        try:
            return (0.599546497 + 52.9690962641 * t) % self.PI2
        except Exception as e:
            raise

    def __lsa_iers2003(self, t):
        """ Saturn longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Saturn longitudes
        """
        try:
            return (0.874016757 + 21.3299104960 * t) % self.PI2
        except Exception as e:
            raise

    def __lur_iers2003(self, t):
        """ Uranus longitudes (IERS 2003)

        :param  float t: ユリウス世紀数
        :return float  : Uranus longitudes
        """
        try:
            return (5.481293872 + 7.4781598567 * t) % self.PI2
        except Exception as e:
            raise

    def __lne_mhb2000(self, t):
        """ Neptune longitude (MHB2000)

        :param  float t: ユリウス世紀数
        :return float  : Neptune longitude
        """
        try:
            return (5.321159000 + 3.8127774000 * t) % self.PI2
        except Exception as e:
            raise

    def __display(self):
        """ Display """
        try:
            print((
                "  [{} TT]\n"
                "  DeltaPsi = {} rad\n"
                "           = {} °\n"
                "           = {} ″\n"
                "  DeltaEps = {} rad\n"
                "           = {} °\n"
                "           = {} ″"
            ).format(
                self.tt.strftime("%Y-%m-%d %H:%M:%S.%f"),
                self.dpsi, self.dpsi_d, self.dpsi_s,
                self.deps, self.deps_d, self.deps_s
            ))
        except Exception as e:
            raise


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

参考までに、スクリプト中に出現する用語の意味は以下のとおり。

1
2
3
4
5
6
7
8
9
10
11
* Mean anomaly of the Moon                                    (平均近点角:月)
* Mean longitude of the Moon minus that of the ascending node (平均黄経:月(昇交点を引いたもの))
* Mean longitude of the Moon's ascending node                 (平均黄経:月の昇交点)
* Mean longitude of Mercury                                   (平均黄経:水星)
* Mean longitude of Venus                                     (平均黄経:金星)
* Mean longitude of Earth                                     (平均黄経:地球)
* Mean longitude of Mars                                      (平均黄経:火星)
* Mean longitude of Jupiter                                   (平均黄経:木星)
* Mean longitude of Saturn                                    (平均黄経:土星)
* Mean longitude of Uranus                                    (平均黄経:天王星)
* General accumulated precession in longitude                 (一般歳差:J2000.0 からの蓄積)

3. Python スクリプトの実行

nut_ls.txt, nut_pl.txt を Python スクリプトと同じディレクトリ内に配置し、 YYYYMMDD or YYYYMMDDHHMMSS or YYYYMMDDHHMMSSUUUUUU で地球時(TT) を指定して実行する。(引数なしでシステム時刻を地球時とみなす)

1
2
3
4
5
6
7
8
$ ./nutation_model.py 20180503123456
  [2018-05-03 12:34:56.000000 TT]
  DeltaPsi = -7.052841833012105e-05 rad
           = -0.0040409807060490495 °
           = -14.547530541776577 ″
  DeltaEps = -2.9031836929708544e-05 rad
           = -0.0016634017275843416 °
           = -5.98824621930363 ″

4. 計算結果の検証

国立天文台・暦象年表のツールで計算した値と比較してみたが、一致することが確認できた。

5. 参考サイト


章動の計算は、天体位置やこよみを正確に計算する際に必要になってきます。

以上。

Comments