Python script to calculate GMST with IAU1982 theory.

calc_gmst_iau_1982.py

#! /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)