mk-mode BLOG

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

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

Python - 平均黄道傾斜角の計算!

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

こんばんは。

以前、 Ruby で平均黄道傾斜角の計算を実装しました。

今回は、同様のことを Python で実現してみました。

0. 前提条件

  • Python 3.6.5 での作業を想定。

1. 計算方法

計算方法等については、過去記事を参照。

2. Python スクリプトの作成

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) を指定して実行する。(引数なしでシステム時刻を地球時とみなす)

1
2
3
4
5
$ ./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