マイコン宇宙講座-月の高度・方位角計算
天体の位置は、赤経・赤緯で表す赤道座標系と高度・方位角で表す地平座標系の２つになります。観測では望遠鏡は赤道儀架台にマウントされていますので、赤道座標系が用いられますが、一般的には地平座標系のほうがわかりやすいものとなります。たとえば、北の方角から時計回りで180°、地平線から高度が45°と表現したほうが、目的の天体を見つけやすいからです。
とはいえ、月くらいに大きな天体になると、地平上に出ていればすぐに見つけられますので、あまり座標系にこだわらなくてもよいのですが、肉眼や双眼鏡で見える比較的明るい天体は、やはり高度・方位角を知っておくと、あそこに見えている天体は？という疑問について簡単に知ることができます。
赤道座標系から地平座標系の変換は、「天体の位置計算 増補版」に計算式が掲載されていますので、関数電卓があれば簡単に求めることができます。プログラムは、「マイコン宇宙講座」第3章の太陽系天体の高度・方位角計算で呼び出したサブルーチン altitude_directionを利用して、月の視位置から求めています。このサブルーチンは月や太陽系天体以外の恒星の高度・方位角計算でも使用することができます。
メインルーチン m56.py
プログラムへの入力は日付と時刻および計算間隔（分単位）で行います。
# m56.py
# マイコン宇宙講座
# 月の高度・方位計算プログラム
import math
import lib
M2PI = 2.0 * math.pi
lg = 139.745 / lib.K[3]
la = 35.654 / lib.K[3]
obs = 'Tokyo'
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)
iv = input('INTERVAL(UNIT OF MINUTE) ? ')
iv = float(iv)
ii = 0
yz = 0
mt = 0
print()
print('%-s (経度:%5.3f 緯度:%3.3f)' % (obs, lg * lib.K[3], la * lib.K[3]))
print(' Date (JST) R.A. Decl. 方位 高度 距離')
print('--------------------------------------------------------------------------------')
print(' 年 月 日 時 分 h m 。 , 。 。 Km')
while ii < 32:
t0 = (jd - 15019.5) / 36525.0
a, b, c, d, e, g, j, l, m, n, v, w = lib.argument(t0)
mx, my, mz = lib.moon(t0, a, b, d, e, n, g, w)
r0 = c
r1 = mx
r2 = my
r3 = mz
ra, dc, dl = lib.positions(r0, r1, r2, r3)
dl1 = dl
if ra < 0:
ra = ra + M2PI
r1 = la
r2 = lib.geocentric_distance(r1)
r0 = r2 / 6378.16
r1 = la
r2 = 0
r3 = lib.geocentric_latitude(r1, r2)
fi = r3
r1 = jd
r2 = lg
r3 = lib.siderealtime_date(r1, r2)
s1 = r3
x0 = r0 * math.cos(fi) * math.cos(s1)
y0 = r0 * math.cos(fi) * math.sin(s1)
z0 = r0 * math.sin(fi)
xx = dl * math.cos(dc) * math.cos(ra)
yy = dl * math.cos(dc) * math.sin(ra)
zz = dl * math.sin(dc)
x = xx - x0
y = yy - y0
z = zz - z0
dl2 = math.sqrt(x * x + y * y + z * z)
ss = y
cc = x
tt = lib.quadrant(ss, cc)
ra = tt
dc = math.atan(z / (y / math.sin(ra)))
al, hi = lib.altitude_direction(jd, lg, la, ra, dc)
# MJD to YY-MM-DD
jd = jd + 0.375
yy, mm, dd = lib.jdate(jd, lib.T)
jd = jd - 0.375
hh = 24.0 * (dd - int(dd))
dd = int(dd)
ms = 60.0 * (hh - int(hh))
hh = int(hh)
ms = int(ms + 0.1)
if ms == 60:
hh = hh + 1
ms = 0
while True:
if hi < -0.17:
break
# Print Position
if ra < 0:
ra = ra + M2PI
if ra > M2PI:
ra = ra - M2PI
ra = ra * lib.K[3] / 15.0
dc = dc * lib.K[3]
dl2 = 6378.16 * dl2
rh = int(ra)
rm = 60.0 * (ra - rh)
sg = lib.sgn(dc)
dc = abs(dc)
dh = int(dc)
dm = 60.0 * (dc - dh)
s1 = '+'
if sg < 0:
s1 = '-'
hi = hi * lib.K[3]
al = al * lib.K[3]
if yy != yz:
sfmt = ' %4d %2d %2d %2d %2d %2d %5.2f %1s%2d %4.1f %5.1f %4.1f %6d' % (yy, mm, dd, hh, ms, rh, rm, s1, dh, dm, al, hi, dl2)
else:
if mm != mt:
sfmt = ' %2d %2d %2d %2d %2d %5.2f %1s%2d %4.1f %5.1f %4.1f %6d' % (mm, dd, hh, ms, rh, rm, s1, dh, dm, al, hi, dl2)
else:
sfmt = ' %2d %2d %2d %2d %5.2f %1s%2d %4.1f %5.1f %4.1f %6d' % (dd, hh, ms, rh, rm, s1, dh, dm, al, hi, dl2)
print(sfmt)
yz = yy
mt = mm
break
ii = ii + 1
if ii < 32:
jd = jd + iv / 1440
print('--------------------------------------------------------------------------------')
print()例題 1981年９月13日16時から30分間隔で月の高度および方位角を求めてみよう。
DATE AND TIME(JST) ? 19810913,160000
INTERVAL(UNIT OF MINUTE) ? 30
Tokyo (経度:139.745 緯度:35.654)
Date (JST) R.A. Decl. 方位 高度 距離
--------------------------------------------------------------------------------
年 月 日 時 分 h m 。 , 。 。 Km
1981 9 13 17 0 22 54.37 -11 22.4 278.4 -7.9 375201
13 17 30 22 55.47 -11 18.1 282.4 -2.1 374462
13 18 0 22 56.53 -11 13.7 286.6 3.7 373732
13 18 30 22 57.53 -11 9.3 291.0 9.3 373019
13 19 0 22 58.49 -11 4.8 295.7 14.8 372335
13 19 30 22 59.39 -11 0.2 300.8 20.1 371688
13 20 0 23 0.26 -10 55.5 306.4 25.1 371088
13 20 30 23 1.08 -10 50.6 312.5 29.7 370543
13 21 0 23 1.87 -10 45.7 319.4 33.9 370060
13 21 30 23 2.63 -10 40.5 327.1 37.6 369646
13 22 0 23 3.37 -10 35.2 335.6 40.5 369307
13 22 30 23 4.09 -10 29.8 344.8 42.6 369046
13 23 0 23 4.80 -10 24.2 354.6 43.8 368866
13 23 30 23 5.50 -10 18.4 4.6 43.9 368769
13 24 0 23 6.21 -10 12.5 14.5 43.0 368754
14 0 30 23 6.93 -10 6.4 23.9 41.2 368821
14 1 0 23 7.67 -10 0.2 32.6 38.4 368968
14 1 30 23 8.42 - 9 53.9 40.5 35.0 369189
14 2 0 23 9.21 - 9 47.4 47.6 30.9 369480
14 2 30 23 10.03 - 9 40.9 54.0 26.4 369835
14 3 0 23 10.89 - 9 34.2 59.8 21.5 370247
14 3 30 23 11.79 - 9 27.5 65.1 16.3 370707
14 4 0 23 12.74 - 9 20.7 70.0 10.9 371206
14 4 30 23 13.74 - 9 13.9 74.6 5.4 371736
14 5 0 23 14.79 - 9 7.1 79.0 -0.3 372285
14 5 30 23 15.89 - 9 0.2 83.2 -6.1 372844
--------------------------------------------------------------------------------