マイコン宇宙講座-太陽と惑星の出没時刻計算
日の出、日の入り時刻に加えて惑星の出没時刻を求めます。N-BASICコードをPythonへ書き直すにあたり、最も苦労しました。というのは処理がGOTO文が使われており、GOTO文がないPythonはループ処理を使わず得ませんでした。しかもループは3重になっており、一番内側から一番外側へジャンプするという処理になっており、それをどう実装するか悩みました。結局、フラグで判断して処理しています。GOTO文があるPHPのほうが実装しやすかったのです。こういうときはGOTO文があるとよいのですね。
さて、出没時刻は0.1分まで求めていますが、実際は1分単位で構いません。秒単位まで必要はないと思います。お天気ニュースでも日の出、日の入りは1分単位ですしね。秒単位の観測は地平線近くでは大気による影響が強すぎて計算した時刻にきちんと観測できません。
サブルーチン appear
対象天体の出没時刻と方位角を計算して返します。
※サブルーチンはlib.pyに記述してください。
# 出没時刻および方位角の計算
def appear(jd, lg, la, ra, dc, ds, rd, sg, hd, K):
if ds > 0:
pp = math.atan(4.26e-5 / ds)
rd = math.atan(rd / (1.5e+8 * ds))
else:
rd = 0.0
pp = 0.0
rr = 0.0102
cc = -math.tan(dc) * math.tan(la)
ss = math.sqrt(1 - cc * cc)
h0 = math.atan(ss / cc)
if h0 < 0:
h0 = h0 + math.pi
dh = (rr + rd - pp) / (math.cos(la) * math.cos(dc) * math.sin(h0))
r1 = jd + sg
r2 = lg
r3 = siderealtime_date(r1, r2)
s1 = r3
if hd == 0:
s2 = ra - h0 - dh
ta = s2 - s1
if ta < 0:
ta = ta + 2.0 * math.pi
ta = ta / 1.002738
ta = K[3] * ta / 15.0
if ta > 24:
ta = ta - 24.0
return ta, 0
if hd == 1:
s3 = ra + h0 + dh
td = s3 - s1
if td < 0:
td = td + 2.0 * math.pi
td = td / 1.002738
td = K[3] * td / 15.0
if td > 24:
td = td - 24
cc = math.sin(dc) / math.cos(la)
ss = -math.cos(dc) * math.sin(h0)
al = abs(math.atan(ss / cc))
if dc < 0:
al = -al + math.pi
return td, al
メインルーチン m312.py
# m312.py
# マイコン宇宙講座
# 3-12 太陽と惑星の出没時刻計算プログラム
import math
import lib
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]
ra = [0 for i in range(10)]
dc = [0 for i in range(10)]
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('計算の間隔 ? ')
pn = input('惑星のナンバー ? ')
iv = float(iv)
pn = int(pn)
ja = jd + 365
print()
print(' %-7s(東経 %7.3f 北緯 %6.3f) での %-7s の出没時刻 (JST)' % (obs, lg * lib.K[3], la * lib.K[3], lib.PL[pn]))
print()
print(' Date 出 没 方位')
print(' --------------------------------------------------------(北からの角度)---')
print(' 年 月 日 時 分 時 分 。')
ha = 0
hd = 0
sg = 0
yz = 0
dz = 0
ta1 = 0
td1 = 0
ta = 0
td = 0
al = 0
mz = 0
# 出没時刻の計算
# L300:
while True:
ty = yy + (mm - 1) / 12 + dd / 365
if hd == 0:
ha = ta1
if hd == 1:
ha = td1
# L400:
while True:
sg = 0
if ha > 15:
sg = -1
t1 = jd + sg + ha / 24.0 - 33281.92334
t1 = t1 * (2.737909288e-5 + 1.260132857e-17 * t1)
t2 = t1 * t1
# L440:
while True:
# EARTH
e, m, p, n, i, a, rd = lib.mean_elements(3, t1, t2)
pe = p
nd = n
ic = i
f, q = lib.eph_const(pe, nd, ic, lib.K)
ec = e
mo = m / (2.0 * math.pi)
mo = 2.0 * math.pi * (mo - int(mo))
ss, cc, ff = lib.kepler(mo, ec)
b = a * math.sqrt(1.0 - e * e)
for n in range(1, 4):
x[n][3] = -1 * (ff * a * f[n] + ss * b * q[n])
# 惑星の平均軌道要素の計算
if pn != 3:
e, m, p, n, i, a, rd = lib.mean_elements(pn, t1, t2)
pe = p
nd = n
ic = i
f, q = lib.eph_const(pe, nd, ic, lib.K)
ec = e
mo = m / (2.0 * math.pi)
mo = 2.0 * math.pi * (mo - int(mo))
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][pn] = x[n][pn] + x[n][3]
# L900:
# 惑星の赤経・赤緯の計算
cc = x[1][pn]
ss = x[2][pn]
tt = lib.quadrant(ss, cc)
ra[pn] = tt
cc = x[1][pn] / math.cos(ra[pn])
ss = x[3][pn]
dc[pn] = math.atan(ss / cc)
r1 = ra[pn]
d1 = dc[pn]
rax, dcx = lib.precession(r1, d1, ty, 0, 0, lib.K)
ra[pn] = rax
dc[pn] = dcx
ds = math.sqrt(x[1][pn]**2 + x[2][pn]**2 + x[3][pn]**2)
rdx = rd
result = lib.appear(jd, lg, la, rax, dcx, ds, rdx, sg, hd, lib.K)
flags_0 = False
flags_1 = False
flags_2 = False
if hd == 1:
td = result[0]
al = result[1]
if abs(td - ha) < 2.0e-3:
flags_2 = True
else:
ha = td
continue
else:
if hd == 0:
# 出
ta = result[0]
if abs(ta - ha) > 2.0e-3:
ha = ta
flags_1 = True
else:
if hd == 0:
hd = 1
ha = ta
flags_0 = True
if flags_0 or flags_1 or flags_2:
break
# end of L440
if flags_0 or flags_2:
break
# end of L400
if flags_0:
continue
if flags_2:
ta1 = ta
ta = ta + 9
if ta > 24:
ta = ta - 24
td1 = td
td = td + 9
if td > 24:
td = td - 24
h1 = int(ta)
m1 = 60.0 * (ta - h1)
h2 = int(td)
m2 = 60.0 * (td - h2)
al = al * lib.K[3]
yy, mm, dd = lib.jdate(jd, lib.T)
sfmt = ' %2d %2d %4.1f %2d %4.1f %6.1f' % (dd, h1, m1, h2, m2, al)
if yy != yz:
sfmt = ' %4d %2d %2d %2d %4.1f %2d %4.1f %6.1f' % (yy, mm, dd, h1, m1, h2, m2, al)
if mm != mz:
sfmt = ' %2d %2d %2d %4.1f %2d %4.1f %6.1f' % (mm, dd, h1, m1, h2, m2, al)
print(sfmt)
yz = yy
mz = mm
if jd < ja:
jd += iv
ha = 0
hd = 0
else:
break
# end of L300
print(' -------------------------------------------------------------------------')
print()例題 1981年1月1日９時00分00秒から10日おきの火星の出没時刻を求めてみよう。
DATE AND TIME(JST) ? 19810101,090000
計算の間隔 ? 10
惑星のナンバー ? 4
東京 (東経 139.745 北緯 35.654) での MARS の出没時刻 (JST)
Date 出 没 方位
--------------------------------------------------------(北からの角度)---
年 月 日 時 分 時 分 。
1 1 8 13.1 18 9.6 116.3
11 7 59.8 18 10.6 113.7
21 7 45.0 18 10.3 110.9
31 7 28.6 18 9.9 107.6
2 10 7 11.1 18 9.4 104.1
20 6 52.6 18 8.5 100.4
3 2 6 33.4 18 7.2 96.6
12 6 13.6 18 5.7 92.7
22 5 53.6 18 3.8 88.8
4 1 5 33.4 18 1.6 85.0
11 5 13.5 17 59.3 81.2
21 4 53.8 17 56.7 77.7
5 1 4 34.7 17 54.1 74.3
11 4 16.3 17 51.2 71.2
21 3 58.8 17 48.0 68.4
31 3 42.3 17 44.4 65.9
6 10 3 27.0 17 40.1 63.9
20 3 12.9 17 35.0 62.2
30 3 0.0 17 28.8 61.1
7 10 2 48.4 17 21.4 60.3
20 2 37.8 17 12.4 60.0
30 2 28.1 17 1.8 60.2
8 9 2 19.1 16 49.5 60.8
19 2 10.6 16 35.5 61.8
29 2 2.3 16 19.9 63.1
9 8 1 53.9 16 2.7 64.8
18 1 45.4 15 44.1 66.7
28 1 36.4 15 24.1 68.8
10 8 1 26.9 15 3.0 71.0
18 1 16.8 14 40.9 73.4
28 1 5.9 14 17.9 75.9
11 7 0 54.1 13 54.0 78.4
17 0 41.5 13 29.5 80.9
27 0 27.8 13 4.2 83.3
12 7 0 12.9 12 38.2 85.6
17 23 54.9 12 11.6 87.9
27 23 36.9 11 44.1 89.9
1 6 23 16.8 11 15.7 91.7
-------------------------------------------------------------------------