マイコン宇宙講座-月の位置計算
最も身近な天体である月は、暦や詩われていたり、月見うどん、月見酒といった食べ物や飲み物の名称、中秋の名月といったイベントに登場してきます。また、月の表面上の模様が餅つきをするうさぎに例えられるなど、月に関する情報は絶えません。天文現象でも月食や日食に加えて、名称としてストロベリームーン、ブルームーン、スーパームーンなどがあります。これだけ話題豊富な天体は珍しいといえます。
さて、太古から知られてきた月ですが、実際の位置計算では、ロマン的な要素はなく複雑な動きから、他の天体と較べて位置を求めるのが大変な天体なのです。
サブルーチン argument
※サブルーチンはlib.pyに記述してください。
def argument(ta):
tb = ta * ta
a = 296.104608 + 477000.0 * ta + 198.849108 * ta + 9.192e-3 * tb
a = fnr(a)
b = 11.250889 + 483120.0 * ta + 82.02515 * ta - 3.211e-3 * tb
b = fnr(b)
c = 270.434164 + 480960.0 * ta + 307.883142 * ta - 1.133e-3 * tb
c = fnr(c)
d = 350.737486 + 444960.0 * ta + 307.114217 * ta - 1.436e-3 * tb
d = fnr(d)
e = 98.998753 + 35640.0 * ta + 359.372886 * ta
e = fnr(e)
g = 358.475833 + 35999.04975 * ta - 1.5e-4 * tb
g = fnr(g)
j = 225.444651 + 2880.0 * ta + 154.906654 * ta
j = fnr(j)
l = 279.696678 + 36000.76892 * ta + 3.03e-4 * tb
l = fnr(l)
m = 319.529425 + 19080.0 * ta + 59.8585 * ta + 1.81e-4 * tb
m = fnr(m)
n = 259.183275 - 1800.0 * ta - 134.142008 * ta + 2.078e-3 * tb
n = fnr(n)
v = 212.603219 + 58320.0 * ta + 197.803875 * ta + 1.286e-3 * tb
v = fnr(v)
w = 342.767053 + 58320.0 * ta + 199.211911 * ta + 3.1e-4 * tb
w = fnr(w)
return a, b, c, d, e, g, j, l, m, n, v, wサブルーチン moon
※サブルーチンはlib.pyに記述してください。
def moon(ta, a, b, d, e, n, g, w):
# MOON-THETA
mx = -1.01e-3 * math.sin(a + b - 4 * d)
mx = mx - 1.02e-3 * math.sin(a - b - 4 * d - n)
mx = mx - 1.03e-3 * ta * math.sin(a - b - n)
mx = mx - 1.07e-3 * math.sin(a - g - b - 2 * d - n)
mx = mx - 1.21e-3 * math.sin(2.0 * a - b - 4 * d - n)
mx = mx + 1.30e-3 * math.sin(3 * a + b + n)
mx = mx - 1.31e-3 * math.sin(a + b - n)
mx = mx + 1.36e-3 * math.sin(a + b - d + n)
mx = mx - 1.45e-3 * math.sin(g + b)
mx = mx - 1.49e-3 * math.sin(a + g - b - 2 * d)
mx = mx + 1.57e-3 * math.sin(g - b + d - n)
mx = mx - 1.59e-3 * math.sin(g - b)
mx = mx + 1.84e-3 * math.sin(a - g + b - 2 * d + n)
mx = mx - 1.94e-3 * math.sin(b - 2 * d - n)
mx = mx - 1.96e-3 * math.sin(g - b + 2 * d - n)
mx = mx + 2.00e-3 * math.sin(b - d)
mx = mx - 2.05e-3 * math.sin(a + g - b)
mx = mx + 2.35e-3 * math.sin(a - g - b)
mx = mx + 2.46e-3 * math.sin(a - 3 * b - n)
mx = mx - 2.62e-3 * math.sin(2 * a + b - 2 * d)
mx = mx - 2.83e-3 * math.sin(a + g + b - 2 * d)
mx = mx - 3.39e-3 * math.sin(g - b - 2 * d - n)
mx = mx + 3.45e-3 * math.sin(a - b + n)
mx = mx - 3.47e-3 * math.sin(g - b + 2 * d)
mx = mx - 3.83e-3 * math.sin(b + d + n)
mx = mx - 4.11e-3 * math.sin(a + g + b + n)
mx = mx - 4.42e-3 * math.sin(2 * a - b - 2 * d - n)
mx = mx + 4.49e-3 * math.sin(a - b + 2 * d)
mx = mx - 4.56e-3 * math.sin(3 * b - 2 * d + n)
mx = mx + 4.66e-3 * math.sin(a + b + 2 * d + n)
mx = mx + 4.9e-3 * math.sin(2 * a - b)
mx = mx + 5.61e-3 * math.sin(2 * a + b)
mx = mx + 5.64e-3 * math.sin(a - g + b + n)
mx = mx - 6.38e-3 * math.sin(a + g - b - n)
mx = mx - 7.13e-3 * math.sin(a + g - b - 2 * d - n)
mx = mx - 9.29e-3 * math.sin(g + b - 2 * d)
mx = mx - 9.47e-3 * math.sin(2 * a - b - n)
mx = mx + 9.65e-3 * math.sin(a - g - b - n)
mx = mx + 9.70e-3 * math.sin(b + 2 * d)
mx = mx + 0.01064 * math.sin(b - d + n)
mx = mx - 0.0125 * ta * math.sin(b + n)
mx = mx - 0.01434 * math.sin(g + b - 2 * d + n)
mx = mx - 0.01652 * math.sin(a + g + b - 2 * d + n)
mx = mx - 0.01868 * math.sin(2 * a + b - 2 * d + n)
mx = mx + 0.02700 * math.sin(2 * a + b + n)
mx = mx - 0.02994 * math.sin(a - b - 2 * d)
mx = mx - 0.03759 * math.sin(g + b + n)
mx = mx - 0.03982 * math.sin(g - b - n)
mx = mx + 0.04732 * math.sin(b + 2 * d + n)
mx = mx - 0.04771 * math.sin(b - n)
mx = mx - 0.06505 * math.sin(a + b - 2 * d)
mx = mx + 0.13622 * math.sin(a + b)
mx = mx - 0.14511 * math.sin(a - b - 2 * d - n)
mx = mx - 0.18354 * math.sin(b - 2 * d)
mx = mx - 0.20017 * math.sin(b - 2 * d + n)
mx = mx - 0.38899 * math.sin(a + b - 2 * d + n)
mx = mx + 0.40248 * math.sin(a - b)
mx = mx + 0.65973 * math.sin(a + b + n)
mx = mx + 1.96763 * math.sin(a - b - n)
mx = mx + 4.95372 * math.sin(b)
mx = mx + 23.89684 * math.sin(b + n)
# MOON-RHO
my = 0.05491 * math.cos(2 * a + g)
my = my + 0.06290 * math.cos(a + d)
my = my - 0.06444 * math.cos(4 * d)
my = my - 0.06652 * math.cos(2 * a - g)
my = my - 0.07369 * math.cos(g - 4 * d)
my = my + 0.08119 * math.cos(a - 3 * d)
my = my - 0.09261 * math.cos(a + 4 * g)
my = my + 0.10177 * math.cos(a - 2 * b + 2 * d)
my = my + 0.10225 * math.cos(a + g + 2 * d)
my = my - 0.10243 * math.cos(a + 2 * g - 2 * d)
my = my - 0.12291 * math.cos(2 * b)
my = my - 0.12291 * math.cos(2 * a - 2 * b)
my = my - 0.12428 * math.cos(a + g - 4 * d)
my = my - 0.14986 * math.cos(3 * a)
my = my - 0.16070 * math.cos(a - g + 2 * d)
my = my - 0.16949 * math.cos(a - d)
my = my + 0.17697 * math.cos(a + 2 * b - 2 * d)
my = my - 0.18815 * math.cos(2 * a - 4 * d)
my = my - 0.19482 * math.cos(2 * g - 2 * d)
my = my + 0.22383 * math.cos(2 * b - 2 * d)
my = my + 0.22594 * math.cos(3 * a - 2 * d)
my = my + 0.24454 * math.cos(2 * a + g - 2 * d)
my = my - 0.31717 * math.cos(g + d)
my = my - 0.36333 * math.cos(a - 4 * d)
my = my + 0.47999 * math.cos(a - g - 2 * d)
my = my + 0.63844 * math.cos(g + 2 * d)
my = my + 0.86170 * math.cos(g)
my = my + 1.50534 * math.cos(a - 2 * b)
my = my - 1.67417 * math.cos(a + 2 * d)
my = my + 1.99463 * math.cos(a + g)
my = my + 2.07579 * math.cos(d)
my = my - 2.45500 * math.cos(a - g)
my = my - 2.74067 * math.cos(a + g - 2 * d)
my = my - 3.83002 * math.cos(g - 2 * d)
my = my - 5.37817 * math.cos(2 * a)
my = my + 6.60763 * math.cos(2 * a - 2 * d)
my = my - 53.97626 * math.cos(2 * d)
my = my - 68.62152 * math.cos(a - 2 * d)
my = my - 395.13669 * math.cos(a)
my = my + 3649.33705
# MOON-PHI
mz = -1.00e-3 * math.sin(a - g - 2 * b - 2 * n)
mz = mz - 1.00e-3 * math.sin(a + g - 4 * d)
mz = mz + 1.00e-3 * math.sin(2 * a - g)
mz = mz + 1.02e-3 * math.sin(a - g + 2 * d)
mz = mz - 1.06e-3 * math.sin(2 * a - 2 * b - n)
mz = mz - 1.06e-3 * math.sin(2 * a + n)
mz = mz - 1.09e-3 * math.sin(a + 2 * b - 2 * d)
mz = mz - 1.10e-3 * math.sin(2 * b - d + 2 * n)
mz = mz + 1.12e-3 * math.sin(4 * d)
mz = mz - 1.22e-3 * math.sin(2 * a - n)
mz = mz - 1.22e-3 * math.sin(2 * a + 2 * b + n)
mz = mz + 1.49e-3 * math.sin(g + 2 * b - 2 * d + 2 * n)
mz = mz - 1.57e-3 * math.sin(2 * a - 4 * d)
mz = mz + 1.71e-3 * math.sin(a + g + 2 * b - 2 * d + 2 * n)
mz = mz - 1.75e-3 * math.sin(2 * a + g - 2 * d)
mz = mz - 1.90e-3 * math.sin(2 * g - 2 * d)
mz = mz + 1.93e-3 * math.cos(a + 16 * e - 18 * w)
mz = mz + 1.94e-3 * math.sin(2 * a + 2 * b - 2 * d + 2 * n)
mz = mz + 2.01e-3 * math.sin(g - 2 * d - n)
mz = mz + 2.01e-3 * math.sin(g + 2 * b - 2 * d + n)
mz = mz - 2.07e-3 * math.sin(a + 2 * g - 2 * d)
mz = mz - 2.10e-3 * math.sin(2 * g)
mz = mz - 2.13e-3 * math.sin(2 * d - n)
mz = mz - 2.13e-3 * math.sin(2 * b + 2 * d + n)
mz = mz - 2.15e-3 * math.sin(3 * a - 2 * d)
mz = mz - 2.47e-3 * math.sin(a - 4 * d)
mz = mz - 2.53e-3 * math.sin(a - 2 * b + 2 * d)
mz = mz + 2.79e-3 * ta * math.sin(2 * b + 2 * n)
mz = mz - 2.80e-3 * math.sin(2 * a + 2 * b + 2 * n)
mz = mz + 3.12e-3 * math.sin(3 * a)
mz = mz - 3.17e-3 * math.sin(a + 2 * b)
mz = mz - 3.50e-3 * math.sin(a + 16 * e - 18 * w)
mz = mz + 3.90e-3 * math.sin(g + 2 * b + 2 * n)
mz = mz + 4.13e-3 * math.sin(g - 2 * b - 2 * n)
mz = mz - 4.90e-3 * math.sin(2 * n)
mz = mz - 4.91e-3 * math.sin(2 * b + 2 * d + 2 * n)
mz = mz + 5.04e-3 * math.sin(g + d)
mz = mz + 5.16e-3 * math.sin(a - d)
mz = mz - 6.21e-3 * math.sin(g + 2 * d)
mz = mz + 6.48e-3 * math.sin(a - 2 * b - 2 * d - n)
mz = mz + 6.48e-3 * math.sin(a - 2 * d + n)
mz = mz + 7.00e-3 * math.sin(a - g - 2 * d)
mz = mz + 0.01122 * math.sin(a + 2 * d)
mz = mz + 0.01410 * math.sin(a - 2 * d - n)
mz = mz + 0.01410 * math.sin(a + 2 * b - 2 * d + n)
mz = mz + 0.01424 * math.sin(a - 2 * b)
mz = mz + 0.01506 * math.sin(a - 2 * b - 2 * d - 2 * n)
mz = mz - 0.01567 * math.sin(2 * b - 2 * d)
mz = mz + 0.02077 * math.sin(2 * b - 2 * d + 2 * n)
mz = mz - 0.02527 * math.sin(a + g)
mz = mz - 0.02952 * math.sin(a - n)
mz = mz - 0.02952 * math.sin(a + 2 * b + n)
mz = mz - 0.03487 * math.sin(d)
mz = mz + 0.03684 * math.sin(a - g)
mz = mz - 0.03983 * math.sin(2 * d + n)
mz = mz + 0.03983 * math.sin(2 * b - 2 * d + n)
mz = mz + 0.04037 * math.sin(a + 2 * b - 2 * d + 2 * n)
mz = mz + 0.04221 * math.sin(2 * a)
mz = mz - 0.04273 * math.sin(g - 2 * d)
mz = mz - 0.05566 * math.sin(2 * a - 2 * d)
mz = mz - 0.05697 * math.sin(a + g - 2 * d)
mz = mz - 0.06846 * math.sin(a + 2 * b + 2 * n)
mz = mz - 0.08724 * math.sin(a - 2 * b - n)
mz = mz - 0.08724 * math.sin(a + n)
mz = mz - 0.11463 * math.sin(2 * b)
mz = mz - 0.18647 * math.sin(g)
mz = mz - 0.20417 * math.sin(a - 2 * b - 2 * n)
mz = mz + 0.59616 * math.sin(2 * d)
mz = mz + 1.07142 * math.sin(n)
mz = mz - 1.07447 * math.sin(2 * b + n)
mz = mz - 1.28658 * math.sin(a - 2 * d)
mz = mz - 2.47970 * math.sin(2 * b + 2 * n)
mz = mz + 6.32962 * math.sin(a)
return mx, my, mzサブルーチン positions
※サブルーチンはlib.pyに記述してください。
def positions(r0, r1, r2, r3):
ss = r3 / math.sqrt(r2 - r1 * r1)
cc = math.sqrt(1.0 - ss * ss)
ra = r0 + math.atan(ss / cc)
ss = r1 / math.sqrt(r2)
cc = math.sqrt(1.0 - ss * ss)
dc = math.atan(ss / cc)
dl = math.sqrt(r2)
return ra, dc, dlメインルーチン m51.py
プログラムの入力は日付と時刻、計算間隔を時間単位入力してください。１日おきであれば24、30分おきであれば0.5というようになります。
# m51.py
# マイコン宇宙講座
# 月の位置計算プログラム
import lib
import math
M2PI = 2.0 * math.pi
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 HOUR) ? ')
iv = float(iv)
ii = 0
yz = 0
md = 0
print()
print(' Date (JST) R.A. (Date) Decl. Distance')
print(' -----------------------------------------------------------------')
print(' 年 月 日 時 分 h m 。 , Km')
while(ii < 32):
ta = (jd - 15019.5) / 36525.0
a, b, c, d, e, g, j, l, m, n, v, w = lib.argument(ta)
mx, my, mz = lib.moon(ta, a, b, d, e, n, g, w)
r0 = c
r1 = mx
r2 = my
r3 = mz
ra, dc, dl = lib.positions(r0, r1, r2, r3)
if ra < 0:
ra = ra + M2PI
ra = ra * lib.K[3] / 15.0
dc = dc * lib.K[3]
dl = 6378.16 * dl
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)
if int(ms + 0.1) == 60:
hh = hh + 1
ms = 0
rh = int(ra)
rm = 60.0 * (ra - rh)
sg = lib.sgn(dc)
dc = abs(dc)
dh = int(dc)
dm = 60.0 * (dc - dh)
s = '+'
if sg < 0:
s = '-'
if yy != yz:
print(' %4d %2d %2d %2d %2d %2d %5.2f %1s%2d %5.1f %6.0f' % (yy, mm, dd, hh, ms, rh, rm, s, dh, dm, dl))
else:
if mm != md:
print(' %2d %2d %2d %2d %2d %05.2f %1s%2d %5.1f %6.0f' % (mm, dd, hh, ms, rh, rm, s, dh, dm, dl))
else:
print(' %2d %2d %2d %2d %5.2f %1s%2d %5.1f %6.0f' % (dd, hh, ms, rh, rm, s, dh, dm, dl))
yz = yy
md = mm
ii = ii + 1
if ii < 32:
jd = jd + iv / 24.0
else:
break
print(' -----------------------------------------------------------------')
print()例題 1981年9月1日21時から1日毎の月の位置を求めてみましょう。
DATE AND TIME(JST) ? 19810901,210000
INTERVAL(UNIT OF HOUR) ? 24
Date (JST) R.A. (Date) Decl. Distance
-----------------------------------------------------------------
年 月 日 時 分 h m 。 , Km
1981 9 1 21 0 12 55.63 - 0 37.1 396711
2 21 0 13 41.48 - 5 0.7 400025
3 21 0 14 27.15 - 9 8.3 402559
4 21 0 15 13.28 -12 51.8 404080
5 21 0 16 0.39 -16 3.2 404420
6 21 0 16 48.90 -18 35.1 403490
7 21 0 17 39.02 -20 19.9 401295
8 21 0 18 30.73 -21 10.6 397939
9 21 0 19 23.78 -21 1.2 393627
10 21 0 20 17.70 -19 48.0 388658
11 21 0 21 11.98 -17 30.8 383403
12 21 0 22 6.19 -14 13.6 378272
13 21 0 23 0.10 -10 4.8 373666
14 21 0 23 53.76 - 5 17.3 369930
15 21 0 0 47.46 - 0 7.4 367306
16 21 0 1 41.68 + 5 6.1 365906
17 21 0 2 36.91 +10 3.3 365710
18 21 0 3 33.53 +14 24.6 366587
19 21 0 4 31.65 +17 52.4 368334
20 21 0 5 30.94 +20 12.3 370724
21 21 0 6 30.66 +21 15.9 373547
22 21 0 7 29.78 +21 1.1 376637
23 21 0 8 27.31 +19 32.5 379884
24 21 0 9 22.57 +17 0.2 383225
25 21 0 10 15.27 +13 37.2 386624
26 21 0 11 5.56 + 9 38.0 390046
27 21 0 11 53.84 + 5 16.5 393436
28 21 0 12 40.65 + 0 45.6 396697
29 21 0 13 26.61 - 3 43.1 399690
30 21 0 14 12.32 - 7 59.2 402239
10 1 21 0 14 58.35 -11 53.7 404144
2 21 0 15 45.18 -15 17.9 405209
-----------------------------------------------------------------