マイコン宇宙講座-基本サブルーチンの作成

日時を指定するとユリウス日を求めます。時刻系は世界時UTから現在はUTCに変わっています。ですが、実際は世界時と変わりません。日本標準時刻と世界時の関係は次のようになります。

UTC=JST-9

世界時は正確な呼称を協定世界時といいます。これはグリニッジ子午線を基準として時刻を表し、天体の位置計算では協定世界時０時における位置を求めることになります。

基本サブルーチンの作成

天文計算において基本となる時刻の計算に使用するユリウス日を求めるサブルーチンを作成します。プログラムはマイコン宇宙講座のものと少し変わっています。当時は8ビット時代であることを踏まえると演算桁数が多いものは避けるべきでしたが、現在ではModify Julian DayではなくJulian Dayとして計算することができます。

サブルーチン ib.pyです。今後記事中にサブルーチンと書かれたプログラムはlib.pyへ記述してください。

サブルーチンは以下のようにメインプログラムから呼び出します。

import lib

特定のサブルーチンを呼び出すには、サブルーチン名の先頭にlibを付けて次のように記述します。

jd = lib.mjd(yy, mm, dd)

サブルーチン lib.py ヘッダー部分

#  マイコン宇宙講座　サブルーチンライブラリ
# Python版
import math


# 月に対応させるためT[0]には0を入れてある。
T = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]

# N-BASICの配列要素に対応
K = [0, 0.01720209895, 0, 180.0 / math.pi, 0.91743695, 0.39788118,
     (180.0 / math.pi) * 3600.0]

# 惑星名
PL = ['', 'MERCURY', 'VENUS', 'EARTH', 'MARS', 'JUPITER', 'SATURN', 'URANUS', 'NEPTUNE', 'PLUTO']


def fno(x):
    return int(x + 0.5)


def fna(x):
    return int(x * 10.0 + 0.5) / 10.0


def fnb(x):
    return int(x * 100.0 + 0.5) / 100.0


def fnc(x):
    return int(x * 1000.0 + 0.5) / 1000.0


def fnr(x):
    return 2.0 * math.pi * (x / 360.0 - int(x / 360.0))


# N-BASICのSGN関数
def sgn(n):
    if n < 0:
        x = -1
    else:
        x = 1
    return x


# N-BASICのTAB関数
def tab(n):
    s = ''
    for i in range(n):
        s += ' '

    return s

サブルーチン mjd

# Modify Julian Day
def mjd(dy, dt):
    s1 = sgn(dy)
    dy = abs(dy)

    yy = int(dy / 10000.0)
    mm = int(dy / 100.0) % 100
    dd = dy % 100
    hh = int(dt / 10000.0)
    ms = int(dt / 100.0) % 100
    ss = dt % 100.0

    if s1 < 0:
        yy *= s1
    else:
        if yy < 100:
            yy += 1900

    jd = julian(yy, mm, dd)

    # ユリウス日に時刻を小数点にした値を加えてJSTからUTに変更
    jd += hh / 24.0 + ms / 1440.0 + ss / 86400.0
    jd -= 0.375

    # ユリウス日(JD)を修正ユリウス日(MJD)にする
    jd -= 2400000.5

    return jd, yy, mm, dd, hh, ms, ss

サブルーチン julian

# Julian Day
# 計算式は「新装改訂版 天文計算入門」(厚生社 1996 長谷川一郎著)より引用
def julian(yy, mm, dd):
    p5 = yy + (mm - 1.0) / 12.0 + dd / 365.25
    if mm < 3:
        mm += 12.0
        yy -= 1.0

    # yy < 0
    jd = sgn(yy) * int(abs(yy) * 365.25) + int(30.59 * (mm - 2)) + dd + 1721085.5

    # yy <= 1582.77
    if p5 >= 0:
        jd = int(yy * 365.25) + int(30.59 * (mm - 2)) + dd + 1721086.5

    # yy >= 1582.78
    if p5 >= 1582.78:
        jd = int(yy * 365.25) + int(yy / 400.0) - int(yy / 100.0) + int(30.59 * (mm - 2)) + dd + 1721088.5

    # N-BASICでは変数がすべてグローバル変数となるため、変更した場合はもとに戻す処理が必要
    # Pythonでは引数を参照渡ししない限り下記の処理は不要。
    yy *= sgn(yy)
    if mm > 12:
        mm -= 12.0
        yy += 1.0

    return jd

サブルーチン jdate

# from Julian Day to Date
def jdate(jd, T):
    jj = jd
    yy = int(2.7379093e-3 * jd + 1858.877)
    mm = 1.0
    dd = 0

    jd = julian(yy, mm, dd) - 2400000.5
    r2 = jj - jd

    # うるう年の判定
    if (yy % 4 == 0 and yy % 100 != 0) or yy % 400 == 0:
        T[2] = 29

    r1 = 0
    m = 1
    while m < 13:
        if int(r2) - r1 - T[m] <= 0:
            break
        r1 += T[m]
        m += 1

    mm = m
    dd = r2 - r1
    T[2] = 28
    jd = jj

    if mm == 13:
        yy += 1.0
        mm -= 12.0

    return yy, mm, dd

次回は、メインプログラムJST⇔MJD変換プログラム、年月日算出プログラム、万年カレンダープログラムなどを掲載していきます。

57