mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages
※2018年9月15日より非力な環境でサーバを運用しているため、各ページの表示に時間がかかる場合があります。ご了承ください。(改良の予定あり)

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

Python - GMST(グリニッジ平均恒星時)の計算(IAU1982理論)!

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

こんばんは。

前回、 Python でグリニッジ平均恒星時等(GMST; Greenwich Mean Sidereal Time)を IAU1982 理論(David Vallado 氏による計算式)を使用して計算しました。

今回は同じ計算を Python で実装してみました。

0. 前提条件

  • Python 3.6.5 での動作を想定。
  • ここでは GMST についての説明はしない。
  • 天文学的な計算については疎いため、誤りがあるかもしれない。

1. Python スクリプトの作成

  • Shebang ストリング(1行目)では、フルパスでコマンド指定している。(当方の慣習
calc_gmst_iau_82.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
#! /usr/local/bin/python3.6
"""
グリニジ平均恒星時 GMST(= Greenwich Mean Sidereal Time)の計算
: IAU1982 版

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

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


class GmstIau82:
    PI2 = math.pi * 2    # => 6.283185307179586
    D2R = math.pi / 180  # => 0.017453292519943295

    def __init__(self):
        self.__get_arg()

    def exec(self):
        try:
            print(self.ut1.strftime("%Y-%m-%d %H:%M:%S UT1"))
            jd_ut1 = self.__gc2jd(self.ut1)
            print("JD(UT1):", jd_ut1)
            gmst = self.__calc_gmst(jd_ut1)
            gmst_d = gmst * 180 / math.pi
            gmst_h = self.__deg2hms(gmst_d)
            print((
                "GMST = {} rad.\n"
                "     = {} deg.\n"
                "     = {}"
            ).format(
                gmst, gmst_d, gmst_h
            ))
        except Exception as e:
            raise

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

    def __gc2jd(self, ut1):
        """ 年月日(グレゴリオ暦) -> ユリウス日(JD) 変換
            : フリーゲルの公式を使用
                JD = int(365.25 * year)
                   + int(year / 400)
                   - int(year / 100)
                   + int(30.59 * (month - 2))
                   + day
                   + 1721088
              ※上記の int(x) は厳密には、x を超えない最大の整数
                (ちなみに、準JDを求めるなら `+ 1721088` が `- 678912` となる)

        :param  datetime ut1: UT1(世界時1)
        :return float     jd: Julian Day
        """
        try:
            year, month, day     = ut1.year, ut1.month, ut1.day
            hour, minute, second = ut1.hour, ut1.minute, ut1.second
            # 1月,2月は前年の13月,14月とする
            if month < 3:
                year  -= 1
                month += 12
            # 日付(整数)部分計算
            jd  = int(365.25 * year) \
                + year // 400 \
                - year // 100 \
                + int(30.59 * (month - 2)) \
                + day \
                + 1721088.5
            # 時間(小数)部分計算
            t  = (second / 3600 \
               + minute / 60 \
               + hour) / 24.0
            return jd + t
        except Exception as e:
            raise

    def __calc_gmst(self, jd_ut1):
        """ GMST(グリニッジ平均恒星時)計算
            : IAU1982理論(by David Vallado)によるもの
                GMST = 18h 41m 50.54841s
                     + 8640184.812866s T + 0.093104s T^2 - 0.0000062s T^3 
                (但し、 T = 2000年1月1日12時(UT1)からのユリウス世紀単位)

        :param  float  jd_ut1: UT1 に対するユリウス日
        :return datetime gmst: グリニッジ平均恒星時(単位:radian)
        """
        try:
            t_ut1= (jd_ut1 - 2451545.0) / 36525
            gmst =  67310.54841 \
                 + (876600.0 * 3600.0 + 8640184.812866 \
                 + (0.093104 \
                 -  6.2e-6 * t_ut1) * t_ut1) * t_ut1
            gmst = (gmst * self.D2R / 240.0) % self.PI2
            if gmst < 0.0:
                gmst += self.PI2
            return gmst
        except Exception as e:
            raise

    def __deg2hms(self, deg):
        """ 99.999° -> 99h99m99s 変換

        :param  float deg: degree
        :return string:    99 h 99 m 99.999 s
        """
        sign = ""
        try:
            h = int(deg / 15)
            _m = (deg - h * 15.0) * 4.0
            m = int(_m)
            s = (_m - m) * 60.0
            if s < 0:
                s *= -1
                sign = "-"
            return "{}{:2d} h {:02d} m {:06.3f} s".format(sign, h, m, s)
        except Exception as e:
            raise


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

2. Python スクリプトの実行

まず、実行権限を付与。

1
$ chmod +x calc_gmst_iau_82.py

そして、コマンドライン引数に UT1(世界時1)を YYYYMMDD[HHMMSS] の書式で指定して実行する。(引数無指定なら、現在時刻を UT1 とみなす)

1
2
3
4
5
6
$ ./calc_gmst_iau_82.py 20180616
2018-06-16 00:00:00 UT1
JD(UT1): 2458285.5
GMST = 4.611451424267976 rad.
     = 264.2167040401474 deg.
     = 17 h 36 m 52.009 s

3. 計算結果の検証

まず、「Ruby - グリニッジ恒星時の計算(IAU2006 理論)!」のスクリプトで計算した結果と比較してみる。
このスクリプトは地球時(TT)を指定して計算するようにしているので正確に比較はできないが、UT1 が近似するように TT を指定して実行して結果を得ると、高い精度で一致することが分かる。

次に、国立天文台・暦象年表のツール「グリニジ恒星時」で計算した値と比較してみる。
時角で2/1000秒の誤差に収まることが確認できた。

また、(当然ながら)前回のブログ記事で紹介した Ruby 版での計算結果と同じになることも確認できた。


以上。

Comments