マイコン宇宙講座-月の視位置計算
「マイコン宇宙講座-月の位置計算」で求めた月の位置は、地心つまり地球の中心から見た位置になります。ですが、実際は地球の表面、地上から見ています。そのため、地球の半径分だけズレることになります。地球に近い天体ほど、その影響を受けることになります。月は地球から約38万Kmほどであり、非常に近い位置にあるといえます。
近くにあるモノをみることで、このズレを体験することができます。最初に左目だけでみた位置を地心からみた位置とすると、右目だけでみると、位置がズレて見えると思います。視位置はこのズレを補正した位置のことをいいます。視位置は最初に地心から求めた位置に、サブルーチン geocentric_latitude、geocentric_distanceを呼び出して補正計算することで視位置を求めることができます。
メインルーチン m54.py
プログラムは東京における視位置ですが、10行目〜12行目の観測地の経緯度を変更することで、東京以外の場所から視位置を求めることができます。
# m54.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 HOUR) ? ')
iv = float(iv)
ii = 0
yz = 0
mt = 0
print()
print(' Date (JST) 地心での月の位置 距離 %-7sでの月の位置 距離' % (obs))
print(' -------------------------------------------------------------------------------')
print(' 年 月 日 時 分 h m 。 , Km 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)
mr = tt
md = math.atan(z / (y / math.sin(mr)))
# 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
# 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]
dl1 = 6378.16 * dl1
h1 = int(ra)
m1 = 60.0 * (ra - h1)
sg = lib.sgn(dc)
dc = abs(dc)
d1 = int(dc)
f1 = 60.0 * (dc - d1)
s1 = '+'
if sg < 0:
s1 = '-'
mr = mr * lib.K[3] / 15.0
md = md * lib.K[3]
dl2 = 6378.16 * dl2
h2 = int(mr)
m2 = 60.0 * (mr - h2)
sg = lib.sgn(md)
md = abs(md)
d2 = int(md)
f2 = 60.0 * (md - d2)
s2 = '+'
if sg < 0:
s2 = '-'
sfmt = ' %2d %2d %2d %2d %5.2f %s%2d %4.1f %6d %2d %5.2f %s%2d %4.1f %6d' % (dd, hh, ms, h1, m1, s1, d1, f1, dl1, h2, m2, s2, d2, f2, dl2)
if yy != yz:
sfmt = ' %4d %2d %2d %2d %2d %2d %5.2f %s%2d %4.1f %6d %2d %5.2f %s%2d %4.1f %6d' % (yy, mm, dd, hh, ms, h1, m1, s1, d1, f1, dl1, h2, m2, s2, d2, f2, dl2)
if mm != mt:
sfmt = ' %2d %2d %2d %2d %2d %5.2f %s%2d %4.1f %6d %2d %5.2f %s%2d %4.1f %6d' % (mm, dd, hh, ms, h1, m1, s1, d1, f1, dl1, h2, m2, s2, d2, f2, dl2)
print(sfmt)
yz = yy
mt = mm
ii = ii + 1
if ii < 32:
jd = jd + iv / 24.0
print(' -------------------------------------------------------------------------------')
print()例題 1981年９月１日 22時、東京における１ヶ月間の月の視位置を求めてみましょう。
DATE AND TIME(JST) ? 19810901,210000
INTERVAL(UNIT OF HOUR) ? 24
Date (JST) 地心での月の位置 距離 Tokyo での月の位置 距離
-------------------------------------------------------------------------------
年 月 日 時 分 h m 。 , Km h m 。 , Km
9 1 21 0 12 55.63 - 0 37.1 396710 12 52.77 - 1 8.9 398263
2 21 0 13 41.48 - 5 0.7 400025 13 38.52 - 5 31.8 400931
3 21 0 14 27.15 - 9 8.3 402558 14 24.17 - 9 40.0 402793
4 21 0 15 13.28 -12 51.8 404080 15 10.35 -13 24.9 403639
5 21 0 16 0.39 -16 3.2 404420 15 57.62 -16 38.8 403317
6 21 0 16 48.90 -18 35.1 403490 16 46.40 -19 13.5 401756
7 21 0 17 39.02 -20 19.9 401295 17 36.90 -21 1.4 398978
8 21 0 18 30.73 -21 10.6 397938 18 29.13 -21 54.7 395105
9 21 0 19 23.78 -21 1.2 393626 19 22.80 -21 47.2 390365
10 21 0 20 17.70 -19 48.0 388657 20 17.42 -20 34.7 385081
11 21 0 21 11.98 -17 30.8 383403 21 12.43 -18 16.7 379648
12 21 0 22 6.19 -14 13.6 378272 22 7.33 -14 57.4 374501
13 21 0 23 0.10 -10 4.8 373666 23 1.87 -10 45.7 370060
14 21 0 23 53.76 - 5 17.3 369930 23 56.07 - 5 54.9 366678
15 21 0 0 47.46 - 0 7.4 367306 0 50.22 - 0 42.3 364591
16 21 0 1 41.68 + 5 6.1 365906 1 44.77 + 4 32.7 363887
17 21 0 2 36.91 +10 3.3 365710 2 40.20 + 9 29.9 364508
18 21 0 3 33.53 +14 24.6 366587 3 36.86 +13 49.6 366275
19 21 0 4 31.65 +17 52.4 368334 4 34.83 +17 14.4 368931
20 21 0 5 30.94 +20 12.3 370723 5 33.76 +19 31.0 372200
21 21 0 6 30.66 +21 15.9 373546 6 32.94 +20 31.6 375826
22 21 0 7 29.78 +21 1.1 376636 7 31.38 +20 15.0 379606
23 21 0 8 27.31 +19 32.5 379884 8 28.17 +18 46.3 383403
24 21 0 9 22.57 +17 0.2 383224 9 22.70 +16 15.3 387134
25 21 0 10 15.27 +13 37.2 386623 10 14.74 +12 55.0 390751
26 21 0 11 5.56 + 9 38.0 390046 11 4.45 + 8 59.1 394218
27 21 0 11 53.84 + 5 16.5 393435 11 52.22 + 4 41.1 397483
28 21 0 12 40.65 + 0 45.6 396697 12 38.61 + 0 13.4 400460
29 21 0 13 26.61 - 3 43.1 399690 13 24.22 - 4 12.8 403029
30 21 0 14 12.32 - 7 59.2 402238 14 9.65 - 8 27.6 405030
10 1 21 0 14 58.35 -11 53.7 404144 14 55.46 -12 21.8 406289
2 21 0 15 45.18 -15 17.9 405209 15 42.15 -15 47.1 406631
-------------------------------------------------------------------------------