Python - 平均黄道傾斜角の計算！
Updated:
以前、 Ruby で平均黄道傾斜角の計算を実装しました。
今回は、同様のことを Python で実現してみました。
0. 前提条件
- Python 3.6.5 での作業を想定。
1. 計算方法
計算方法等については、過去記事を参照。
2. Python スクリプトの作成
File: mean_obliquity_ecliptic.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
#! /usr/local/bin/python3.6
"""
平均黄道傾斜角の計算
date name version
2018.04.09 mk-mode.com 1.00 新規作成
Copyright(C) 2018 mk-mode.com All Rights Reserved.
---
引数 : 日時(TT（地球時）)
書式：YYYYMMDD or YYYYMMDDHHMMSS or YYYYMMDDHHMMSSUUUUUU
無指定なら現在(システム日時)を地球時とみなす。
"""
from datetime import datetime
import re
import sys
import traceback
class MeanObliquityEcliptic:
def __init__(self):
self.__get_arg()
def exec(self):
""" 実行 """
try:
self.__calc_jd()
self.__calc_t()
self.__calc_eps_a()
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]):
tt = sys.argv[1].ljust(20, "0")
else:
sys.exit(0)
try:
self.tt = datetime.strptime(tt, "%Y%m%d%H%M%S%f")
except ValueError as e:
print("Invalid date!")
sys.exit(0)
except Exception as e:
raise
def __calc_jd(self):
""" ユリウス日の計算
* 地球時 self.tt のユリウス日を計算し、self.jd に設定
"""
year, month, day = self.tt.year, self.tt.month, self.tt.day
hour, minute, second = self.tt.hour, self.tt.minute, self.tt.second
second += self.tt.microsecond
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
self.jd = d + t
except Exception as e:
raise
def __calc_t(self):
""" ユリウス世紀数の計算
* ユリウス日 self.jd のユリウス世紀数を計算し、 self.t に設定
"""
try:
self.t = (self.jd - 2451545) / 36525
except Exception as e:
raise
def __calc_eps_a(self):
""" 黄道傾斜角(εA)の計算
* ユリウス世紀数 self.t から黄道傾斜角 ε （単位: rad）を計算し、
self.eps_a に設定
* 以下の計算式により求める。
ε = 84381.406 - 46.836769 * T - 0.0001831 T^2 + 0.00200340 T^3
- 5.76 * 10^(-7) * T^4 - 4.34 * 10^(-8) * T^5
ここで、 T は J2000.0 からの経過時間を 36525 日単位で表したユリウス
世紀数で、 T = (JD - 2451545) / 36525 である。
"""
t = self.t
try:
self.eps_a = (84381.406 \
+ ( -46.836769 \
+ ( -0.0001831 \
+ ( 0.00200340 \
+ ( -5.76e-7 \
+ ( -4.34e-8) \
* t) * t) * t) * t) * t) / 3600
except Exception as e:
raise
def __display(self):
""" 結果表示 """
try:
print((
" 地球時(TT): {}\n"
" ユリウス日(JD): {}\n"
" ユリウス世紀数(JC): {}\n"
"平均黄道傾斜角(eps_a): {} °"
).format(
datetime.strftime(self.tt, "%Y-%m-%d %H:%M:%S.%f"),
self.jd, self.t, self.eps_a
))
except Exception as e:
raise
if __name__ == '__main__':
try:
obj = MeanObliquityEcliptic()
obj.exec()
except Exception as e:
traceback.print_exc()
sys.exit(1)
- Gist - Python script to calc a Mean Obliquity of the Ecliptic.
__calc_eps_a()メソッド内の t のべき乗計算部分を特殊な記述にしているのは、演算コストのかかる乗算回数を減らして高速化を図るため。（例えば、1 + 2 * t + 3 * t**2 + 4 * t**3だと乗算回数が「6回」だが、1 + (2 + (3 + 4 * t) * t) * tと書き換えることで乗算回数が「3回」になる）
3. Python スクリプトの実行
YYYYMMDD or YYYYMMDDHHMMSS or YYYYMMDDHHMMSSUUUUUU で地球時(TT) を指定して実行する。（引数なしでシステム時刻を地球時とみなす）
$ ./mean_obliquity_ecliptic.py 20180503123456123456
地球時(TT): 2018-05-03 12:34:56.000000
ユリウス日(JD): 2458242.0242592595
ユリウス世紀数(JC): 0.18335453139656285
平均黄道傾斜角(eps_a): 23.436893964544698 °
4. 計算結果の検証
国立天文台・暦象年表のツールで計算した値と比較してみたが、かなりの精度で一致することが確認できた。
以上。
Comments