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)
|