マイコン宇宙講座-彗星の位置予報
明るく、肉眼でも見える彗星は、それだけで華やかな天文現象です。いままでに数々の彗星が天空を賑わせてきました。残念ながら近年は、昔のように華々しい彗星の出現はないのは寂しい気がします。かっては眼視による彗星捜索は観測機器の進歩によって写真に取って代わられ、また、SOHOで撮影された写真から数々の太陽に接近する彗星が発見されています。時代の流れといえばそれまでですが、眼視による発見で、天空を飾る池谷・関彗星のような発見はもうこないのかもしれません。
コンピュータの進化によって、天体の位置計算は大きく変わってきました。今のようにコンピュータがない時代は、三角関数表や対数表を用いて手で計算していました。観測結果から軌道要素の計算、位置予報ならともかく摂動を含む計算は、数ヶ月や年単位の計算時間が必要だったと聞きます。位置予報でも大変な計算ですが。今では、コンピュータによって、昔のように時間をかけることなく瞬時に求めることができます。もちろん、計算だけでなくコンピュータによる天体の自動導入、画像から位置計算や測光といった処理をアマチュアでも行えるようになり、コンピュータはいまやとってはならないものになっています。
サブルーチン Nearly Parabolic
※サブルーチンはlib.pyに記述してください。
# NEARLY PARABOLIC
def nearly_parabolic(tp, t, q, ec, K):
r1 = 1.0 + 9.0 * ec
aa = math.sqrt(0.1 * r1)
bb = 5.0 * (1.0 - ec) / r1
cc = math.sqrt(5.0 * (1.0 + ec) / r1)
b = 1.0
a0 = 0.0
while True:
u = b * aa * K[1] * (tp - t) / (math.sqrt(2.0) * math.pow(q, 1.5))
r2 = 1.0
while True:
v = (u + 2.0 * math.pow(r2, 3) / 3.0) / (1.0 + r2 * r2)
if abs(v - r2) > 1.0e-6:
r2 = v
else:
break
a = bb * v * v
a2 = a * a
a3 = a2 * a
b = 1.0 - 0.017142857 * a2 - 0.003809524 * a3
if abs(a - a0) > 1.0e-6:
a0 = a
else:
break
c = 1.0 + 0.4 * a + 0.21714286 * a2 + 0.12495238 * a3
d = 1.0 - a + 0.2 * a2 + 0.00571429 * a3
qd = q * d
v = cc * c * v
cc = qd * (1.0 - math.pow(v, 2))
ss = 2.0 * qd * v
v = 2.0 * math.atan(v)
return ss, cc, vメインルーチン m61.py
プログラム実行時に求められる日付と時刻については、近日点通過時刻の前後を指定するようにしてください。ハレー彗星のように近日点通過時刻が1986年2月9日の場合、現在の位置を求めようとして2022年7月1日というように入力すると、尾の見かけ上の長さを計算中に平方根の値が範囲外（マイナス値）になってしまうため、プログラムがブレークします。その距離の遠さから計算結果が０より小さくなってしまうからだとは思います。これを回避するために、０より小さくなった場合は、尾の見かけ上の長さを０とするように処理を追加しています。なお、計算は1950年分点で計算されていますので、2022年を指定しても位置は正しくないと思います。
# m61.py
# マイコン宇宙講座
# 6-1 彗星の位置予報プログラム
import math
import lib
M2PI = 2.0 * math.pi
x = [[0 for i in range(10)] for j in range(4)]
f = [0, 0, 0, 0]
q = [0, 0, 0, 0]
fx = [0, 0, 0, 0] # f[n + 3]
qx = [0, 0, 0, 0] # q[n + 3]
sx = [0, 0, 0, 0] # s[n]
ex = [0, 0, 0, 0] # e[n]
tx = [0, 0, 0, 0]
elm = [0 for i in range(11)]
# 彗星の軌道要素
# 他の彗星の場合は、ここのデータを入れ替える
# ただし、離心率が0.8を超える彗星のみ
def OrbitElements():
elm = [0 for i in range(11)]
c = 'COMET P/Halley'
elm[1] = 1986.0 # 近日点通過時刻 T
elm[2] = 2.0
elm[3] = 9.6613
elm[4] = 0.587096 # 近日点距離 q
elm[5] = 0.967267 # 離心率 e
elm[6] = 111.8534 # 近日点引数 ω
elm[7] = 58.1531 # 昇交点黄経 Ω
elm[8] = 162.2378 # 軌道傾斜角 i
elm[9] = 5.2 # 全光度 m1
elm[10] = 0.2 # 尾
return c, elm
print()
std = input('計算を始めたい日と時分(JST) ? ')
dy, dt = std.split(',')
dy = float(dy)
dt = float(dt)
jd, yy, mm, dd, hh, ms, ss = lib.mjd(dy, dt)
jd = 10.0 * int(jd / 10.0)
jd1 = jd
std = input('計算期間(日単位) ? ')
ii = int(std)
std = input('計算間隔(日単位) ? ')
iv = int(std)
yz = 0
mz = 0
jd2 = jd + ii
# 彗星の位置予報
# 彗星の軌道要素の読み込み
c, elm = OrbitElements()
yy = elm[1]
mm = elm[2]
dd = elm[3]
pe = elm[6]
nd = elm[7]
ic = elm[8]
m1 = elm[9]
tl = elm[10]
jd = lib.julian(yy, mm, dd) - 2400000.5
t = jd
pe = pe / lib.K[3]
nd = nd / lib.K[3]
ic = ic / lib.K[3]
f, q = lib.eph_const(pe, nd, ic, lib.K)
for i in range(1, 4):
fx[i] = f[i]
qx[i] = q[i]
no = lib.K[1] / math.pow(elm[4] / (1.0 - elm[5]), 1.5)
a = elm[4] / (1.0 - elm[5])
pd = pow(a, 1.5)
print()
print(' ' + c)
print()
print(' T = %4d %2d %7.4f ET' % (yy, mm, dd))
print(' Peri.= %8.5f e = %10.6f' % (pe * lib.K[3], elm[5]))
print(' Node = %9.5f (1950) a = %9.5f AU' % (nd * lib.K[3], a))
print(' Inc. = %8.5f n = %9.6f' % (ic * lib.K[3], no * lib.K[3]))
print(' q = %9.6f (AU) P = %4.2f 年' % (elm[4], pd))
print()
print('Date (0時ET) R.A. (1950) Decl. Delta r Mag Elong. PA Tail')
print('----------------------------------------------------------------------------------')
print(' 年 月 日 h m 。 , AU AU 等 。 。 。')
# 太陽座標(X,Y,Z)
while True:
t1 = jd1 - 33281.92334
t1 = t1 * (2.737909288e-5 + 1.260132857e-17 * t1)
t2 = t1 * t1
e, m, p, n, i, a, rd = lib.mean_elements(3, t1, t2) # GOSUB 5300
pe = p
nd = n
ic = i
f, q = lib.eph_const(pe, nd, ic, lib.K) # GOSUB 4000
ec = e
mo = m / M2PI
mo = M2PI * (mo - int(mo))
ss, cc, ff = lib.kepler(mo, ec) # GOSUB 8500
b = a * math.sqrt(1.0 - ec * ec)
for n in range(1, 4):
f[n] = a * f[n]
q[n] = b * q[n]
rs = 0
for n in range(1, 4):
x[n][3] = -(ff * f[n] + ss * q[n])
rs = rs + math.pow(x[n][3], 2)
rs = math.sqrt(rs) # 地球と太陽との距離
# 彗星の位置(赤経・赤緯)
dl = 0
while True:
tp = jd1 - 0.005776 * dl
ss, cc, v = lib.nearly_parabolic(tp, t, elm[4], elm[5], lib.K)
rc = 0
ds = 0
for n in range(1, 4):
sx[n] = fx[n] * cc + qx[n] * ss
ex[n] = sx[n] + x[n][3]
rc = rc + math.pow(sx[n], 2)
ds = ds + math.pow(ex[n], 2)
rc = math.sqrt(rc) # 太陽と彗星の距離
ds = math.sqrt(ds) # 地球と彗星の距離
if abs(ds - dl) > 0.001:
dl = ds
else:
break
ss = ex[2]
cc = ex[1]
tt = lib.quadrant(ss, cc) # GOSUB 3500
ra = tt
dc = math.atan(ex[3] * math.cos(ra) / ex[1])
mg = m1 + 5.0 * math.log(ds) / 2.30259 + 15.0 * math.log(rc) / 2.30259
cc = (rs * rs + ds * ds - rc * rc) / (2.0 * rs * ds)
ss = math.sqrt(1.0 - cc * cc)
eg = math.atan(ss / cc)
if eg < 0:
eg = eg + math.pi
dl = 0
for n in range(1, 4):
tx[n] = fx[n] * (rc + tl) * math.cos(v) + qx[n] * (rc + tl) * math.sin(v)
tx[n] = tx[n] + x[n][3]
dl = dl + math.pow(tx[n], 2)
ss = tx[2]
cc = tx[1]
tt = lib.quadrant(ss, cc) # GOSUB 3500
rt = tt
dt = math.atan(tx[3] * math.cos(rt) / tx[1])
dl = math.sqrt(dl)
cc = (math.pow(ds, 2) + math.pow(dl, 2) - math.pow(tl, 2)) / (2.0 * ds * dl)
# (1.0 - cc * cc)の計算結果が０より小さい場合は見かけ上の尾の長さを０とする
if (1.0 - cc * cc) < 0:
cl = 0
else:
ss = math.sqrt(1.0 - cc * cc)
cl = math.atan(ss / cc)
ss = math.cos(dt) * math.sin(rt - ra)
cc = math.cos(dc) * math.sin(dt) - math.sin(dc) * math.cos(dt) * math.cos(rt - ra)
tt = lib.quadrant(ss, cc) # gosub 3500
pa = tt
# 位置予報
jd = jd1
yy, mm, dd = lib.jdate(jd, lib.T)
ra = ra * lib.K[3] / 15.0
dc = dc * lib.K[3]
rh = int(ra)
rm = 60.0 * (ra - rh)
if int(rm + 1.0e-3) == 60:
rh = rh + 1
rm = 0
sg = lib.sgn(dc)
dc = abs(dc)
dh = int(dc)
dm = 60.0 * (dc - dh)
if int(dm + 0.01) == 60:
dc = dc + 1
dm = 0
str_sg = '+'
if sg < 0:
str_sg = '-'
eg = eg * lib.K[3]
cl = cl * lib.K[3]
pa = pa * lib.K[3]
sfmt = ' %2d %2d %5.2f %s%2d %4.1f %6.3f %6.3f %4.1f %5.1f %5.1f %4.1f' % (dd, rh, rm, str_sg, dh, dm, ds, rc, mg, eg, pa, cl)
if yy != yz:
sfmt = '%4d %2d %2d %2d %5.2f %s%2d %4.1f %6.3f %6.3f %4.1f %5.1f %5.1f %4.1f' % (yy, mm, dd, rh, rm, str_sg, dh, dm, ds, rc, mg, eg, pa, cl)
else:
if mm != mz:
sfmt = ' %2d %2d %2d %5.2f %s%2d %4.1f %6.3f %6.3f %4.1f %5.1f %5.1f %4.1f' % (mm, dd, rh, rm, str_sg, dh, dm, ds, rc, mg, eg, pa, cl)
print(sfmt)
yz = yy
mz = mm
jd1 = jd1 + iv
if jd1 < jd2:
continue
else:
break
print('-----------------------------------------------------------------------------------')
print()例題 1985年8月1日から1年間のハレー彗星の位置予報を計算してみよう。なお、時刻は９時とします。９時とするのは、位置予報は世界時０時で発表されるためです。
計算を始めたい日と時分(JST) ? 19850801,090000
計算期間(日単位) ? 365
計算間隔(日単位) ? 10
COMET P/Halley
T = 1986 2 9.6613 ET
Peri.= 111.85340 e = 0.967267
Node = 58.15310 (1950) a = 17.93591 AU
Inc. = 162.23780 n = 0.012975
q = 0.587096 (AU) P = 75.96 年
Date (0時ET) R.A. (1950) Decl. Delta r Mag Elong. PA Tail
----------------------------------------------------------------------------------
年 月 日 h m 。 , AU AU 等 。 。 。
1985 7 24 5 47.18 +18 44.0 3.991 3.198 15.8 33.9 261.7 0.5
8 3 5 53.46 +18 54.6 3.761 3.081 15.4 41.9 264.3 0.6
13 5 59.36 +19 4.0 3.510 2.962 15.0 50.0 266.3 0.8
23 6 4.64 +19 12.8 3.238 2.841 14.6 58.3 267.8 1.0
9 2 6 9.00 +19 21.7 2.950 2.717 14.1 66.9 269.2 1.2
12 6 12.02 +19 31.7 2.646 2.592 13.5 75.9 270.2 1.5
22 6 13.07 +19 44.5 2.332 2.464 12.9 85.4 271.0 1.8
10 2 6 11.16 +20 2.0 2.012 2.333 12.2 95.6 271.4 2.2
12 6 4.60 +20 27.1 1.690 2.199 11.5 107.0 271.4 2.7
22 5 50.34 +21 2.5 1.374 2.063 10.6 120.3 270.4 3.1
11 1 5 22.48 +21 46.7 1.077 1.924 9.6 136.7 267.7 3.2
11 4 29.84 +22 11.9 0.822 1.781 8.5 158.8 260.3 2.3
21 2 59.08 +20 18.9 0.654 1.635 7.5 169.6 91.2 1.5
12 1 1 6.79 +13 43.3 0.632 1.486 6.8 132.1 73.0 7.0
11 23 40.92 + 6 5.0 0.751 1.334 6.5 99.6 67.7 9.3
21 22 49.61 + 0 55.7 0.938 1.180 6.1 75.7 65.5 8.7
31 22 18.39 - 2 14.1 1.139 1.025 5.6 57.2 63.7 7.3
1986 1 10 21 56.71 - 4 23.9 1.323 0.875 4.9 41.4 61.1 5.8
20 21 38.70 - 6 12.5 1.469 0.737 4.1 26.7 55.9 4.2
30 21 21.11 - 8 5.8 1.553 0.632 3.2 12.9 39.4 2.3
2 9 21 2.79 -10 20.3 1.551 0.587 2.7 7.6 308.4 1.5
19 20 44.70 -13 1.2 1.450 0.622 2.9 20.1 268.3 3.9
3 1 20 27.45 -16 12.4 1.269 0.721 3.6 34.5 260.3 6.4
11 20 8.51 -20 19.6 1.039 0.855 4.3 49.7 257.7 8.9
21 19 39.78 -26 33.2 0.787 1.005 4.7 67.4 258.2 11.9
31 18 34.80 -37 13.4 0.549 1.159 4.9 92.3 267.0 14.8
4 10 15 25.04 -47 29.1 0.417 1.314 5.1 131.0 313.7 11.2
20 12 4.69 -32 58.2 0.515 1.466 6.2 147.8 53.4 6.0
30 10 58.28 -19 16.8 0.770 1.615 7.8 130.3 92.9 5.7
5 10 10 35.02 -12 25.4 1.078 1.762 9.1 115.0 103.5 4.8
20 10 26.27 - 8 49.9 1.404 1.905 10.1 102.9 108.4 3.8
30 10 23.93 - 6 50.4 1.732 2.045 11.1 92.5 111.5 3.0
6 9 10 24.96 - 5 44.3 2.058 2.181 11.8 83.0 113.9 2.4
19 10 27.98 - 5 11.1 2.374 2.315 12.5 74.2 116.1 1.9
29 10 32.24 - 4 59.7 2.679 2.446 13.2 65.9 118.3 1.5
7 9 10 37.30 - 5 3.8 2.969 2.575 13.7 57.8 120.6 1.2
19 10 42.86 - 5 19.2 3.241 2.701 14.2 49.9 123.3 1.0
-----------------------------------------------------------------------------------