マイコン宇宙講座-太陽系天体の高度・方位角計算

惑星の位置は赤経・赤緯以外にも高度と方位角で表すことができます。赤経・赤緯から高度・方位角を求めます。計算に必要なものは、赤経・赤緯および恒星時、緯度です。計算式については、「天体の位置計算 増補版」29頁の式2−7を使用します。実際の計算のしかたを同書30頁に掲載されていますので、電卓なりプログラムを組むなりして計算してみるとよいでしょう。高度・方位角の計算はサブルーチン altitude_directionで行います。

サブルーチン altitude_direction

※サブルーチンはlib.pyに記述してください。

# KODO_HOUI
def altitude_direction(jd, lg, la, ra, dc):
    r1 = jd
    r2 = lg

    r3 = siderealtime_date(r1, r2)
    s1 = r3 - ra
    cc = -math.sin(dc) * math.cos(la) + math.cos(dc) * math.sin(la) * math.cos(s1)
    ss = math.cos(dc) * math.sin(s1)

    tt = quadrant(ss, cc)

    al = tt
    cc = cc / math.cos(al)
    ss = math.sin(dc) * math.sin(la) + math.cos(dc) * math.cos(la) * math.cos(s1)
    hi = math.atan(ss / cc)

    return al, hi

メインルーチン m311.py

# m311.py
# マイコン宇宙講座
# 3-11 太陽系天体の高度・方位角計算プログラム
import math
import lib


M2PI = 2.0 * math.pi

lg = 139.745 / lib.K[3]
la = 35.654 / lib.K[3]
obs = '東京'

lib.PL[3] = 'Sun'

x = [[0 for i in range(10)] for j in range(4)]
f = [0, 0, 0, 0]
q = [0, 0, 0, 0]

print()
std = input('DATE AND TIME(JST) ? ')
dy, dt = std.split(',')
dy = float(dy)
dt = float(dt)

jd, yy, mm, dd, hh, ms, ss = lib.mjd(dy, dt)

print()
print('  Date= %5d 年 %2d 月 %2d 日  Time= %2d 時 %2d 分 %2.1f 秒' % (yy, mm, dd, hh, ms, ss))
print('         MJD=%12.2f ET                Place=%-s' % (jd, obs))
print(' Planets        R.A.   (Date)   Decl.         方位      高度       Distance')
print(' ------------------------------------------(南より西)----------------------')
print('                h   m            。  ,           。        。            AU')

t1 = jd - 33281.92334
t1 = t1 * (2.737909288e-5 + 1.260132857e-17 * t1)
t2 = t1 * t1

ty = yy + (mm - 1) / 12.0 + dd / 365.0
if yy < 0:
    ty = yy - (1 - ((mm - 1) / 12.0 + dd / 365.0))

for pn in range(1, 10):
    e, m, p, n, i, a, rd = lib.mean_elements(pn, t1, t2)
    f, q = lib.eph_const(p, n, i, lib.K)

    ec = e
    mo = M2PI * (m / (M2PI) - int(m / (M2PI)))

    ss, cc, ff = lib.kepler(mo, ec)

    b = a * math.sqrt(1.0 - e * e)
    for n in range(1, 4):
        f[n] = a * f[n]
        q[n] = b * q[n]

    # 惑星座標の計算
    for n in range(1, 4):
        x[n][pn] = ff * f[n] + ss * q[n]

for n in range(1, 4):
    x[n][3] = - (x[n][3])

for pn in range(1, 10):
    if pn != 3:
        for n in range(1, 4):
            x[n][pn] = x[n][pn] + x[n][3]

# 惑星の赤経・赤緯の計算
for pn in range(1, 10):
    cc = x[1][pn]
    ss = x[2][pn]

    tt = lib.quadrant(ss, cc)

    ra = tt
    cc = x[1][pn] / math.cos(ra)
    ss = x[3][pn]
    dc = math.atan(ss / cc)
    r1 = ra
    d1 = dc

    ra1, dc1 = lib.precession(r1, d1, ty, 0, 0, lib.K)

    ra = ra1
    dc = dc1

    al, hi = lib.altitude_direction(jd, lg, la, ra, dc)

    ra *= lib.K[3] / 15.0
    dc *= lib.K[3]
    ds = math.sqrt(x[1][pn]**2 + x[2][pn]**2 + x[3][pn]**2)

    rh = int(ra)
    rm = 60.0 * (ra - rh)
    dh = int(abs(dc))
    dm = 60.0 * (abs(dc) - dh)
    dh = lib.sgn(dc) * dh
    al *= lib.K[3]
    hi *= lib.K[3]

    if hi < 0:
        # hind = '[Can not see]'
        hind = '---.-     ---.-'
        print(' %-7s      %02d  %05.2f       %+03d  %04.1f       %-15s      %8.5f' % (lib.PL[pn], rh, rm, dh, dm, hind, ds))
    else:
        print(' %-7s      %02d  %05.2f       %+03d  %04.1f       %5.1f     %5.1f      %8.5f' % (lib.PL[pn], rh, rm, dh, dm, al, hi, ds))

print(' --------------------------------------------------------------------------')
print()

例題 1981年１月１日６時30分における惑星の高度と方位角を求めてみよう。

DATE AND TIME(JST) ? 19810101,063000

  Date=  1981 年  1 月  1 日  Time=  6 時 30 分 0.0 秒
         MJD=    44604.90 ET                Place=東京
 Planets        R.A.   (Date)   Decl.         方位      高度       Distance
 ------------------------------------------(南より西)----------------------
                h   m            。  ,           。        。            AU
 MERCURY      18  47.08       -24  47.2       ---.-     ---.-       1.43796
 VENUS        17  03.95       -21  56.9       310.1      13.5       1.51182
 Sun          18  45.21       -23  02.1       ---.-     ---.-       0.98330
 MARS         20  12.99       -21  05.4       ---.-     ---.-       2.26756
 JUPITER      12  36.42       -02  30.1        21.3      49.8       5.33757
 SATURN       12  40.13       -01  45.7        20.3      50.8       9.46859
 URANUS       15  47.45       -19  44.1       324.0      25.9      19.50573
 NEPTUNE      17  26.96       -21  54.6       306.2       9.8      31.20136
 PLUTO        13  56.75       +06  13.7       347.0      60.0      30.29481
 --------------------------------------------------------------------------

57