マイコン宇宙講座-恒星の歳差計算

恒星の固有運動から歳差を考慮して新しい恒星の位置を計算します。地球は自転軸に沿って回転していますが、いつも同じではなくコマの首振り運動のようにずれを起こします。この周期は約２万６千年で１周します。この運動を歳差運動といい、地球上から恒星の位置を観測すると変化していることがわかります。また、恒星はいつも同じ位置にいるわけではなく、赤経および赤緯方向に固有運動を持ちこれも位置が変化していきます。したがって恒星の位置を正しく計算するには、歳差および固有運動を含めて計算する必要があります。厳密にはこれだけでは不足なのですが。

計算は恒星の位置を角度に換算してから指定します。例題のスピカの赤経および赤緯を角度に変換してみましょう。

赤経α(1950)：13h22m33.3s → 15 * (13 + 22 / 60 + 33.3 / 3600) = 200.63875
赤緯δ(1950)：-10°54′03″ → -(10 + 54 / 60 + 03 / 3600) = −10.90083…

サブルーチン precession

※サブルーチンはlib.pyに記述してください。

# SAISA
def precession(r1, d1, ty, p1, p2, K):
    p1 /= K[6]
    p2 /= K[6]
    ep = 1950.0
    t0 = (ep - 1900.0) / 100.0
    t1 = (ty - ep) / 100.0
    t2 = t1 * t1
    t3 = t2 * t1

    # 天体の位置計算 長沢工著 53頁　式3-11より
    x1 = (2304.25 + 1.396 * t0) * t1 + 0.302 * t2 + 0.018 * t3
    x2 = x1 + 0.791 * t2 + 0.001 * t3
    x3 = (2004.682 - 0.853 * t0) * t1 - 0.426 * t2 - 0.042 * t3

    x1 /= K[6]
    x2 /= K[6]
    x3 /= K[6]

    r2 = r1 + p1 * t1
    d2 = d1 + p2 * t1

    ss = math.cos(d2) * math.sin(r2 + x1)
    cc = math.cos(x3) * math.cos(d2) * math.cos(r2 + x1) - math.sin(x3) * math.sin(d2)

    r3 = quadrant(ss, cc)

    cc = ss / math.sin(r3)
    ss = math.cos(x3) * math.sin(d2) + math.sin(x3) * math.cos(d2) * math.cos(r2 + x1)
    dc = math.atan(ss / cc)
    ra = r3 + x2

    return ra, dc

サブルーチン quadrant

# SYOGEN
def quadrant(ss, cc):
    tt = math.atan(ss / cc)
    if cc < 0:
        tt += math.pi
    else:
        if ss < 0:
            tt += (2.0 * math.pi)

    return tt

メインルーチン m26.py

# m26.py
# マイコン宇宙講座
# 2-6 恒星の歳差計算プログラム
import math
import lib


print()
yy = int(input('計算を始めたい年 ? '))
ra = float(input('恒星の赤経 ? ')) / lib.K[3]
dc = float(input('恒星の赤緯 ? ')) / lib.K[3]
pr = float(input('赤経方向の固有運動量 ? '))
pd = float(input('赤緯方向の固有運動量 ? '))
la = float(input('緯度 ? ')) / lib.K[3]

print()
print('  YEAR   R.A.     DECL.   DIRECTION(N to E)')
print('             。       。         。')

r1 = ra
d1 = dc
for y in range(yy, 1700, 100):
    ty = y + 0.25
    ra, dc = lib.precession(r1, d1, ty, pr, pd, lib.K)

    cc = math.sin(dc) / math.cos(la)
    ss = math.sqrt(1.0 - cc * cc)
    al = math.atan(ss / cc)
    if al < 0:
        al += math.pi
    al *= lib.K[3]
    ra *= lib.K[3]
    dc *= lib.K[3]

    print('  %4d   %7.2f  %7.2f     %6.2f' % (ty, ra, dc, al))

print()

例題 スピカの300年から1600年までの位置の変化を求めてみよう。

計算を始めたい年 ? 300
恒星の赤経 ? 200.63875
恒星の赤緯 ? -10.90083
赤経方向の固有運動量 ? -4.4
赤緯方向の固有運動量 ? -3.3
緯度 ? 50

  YEAR   R.A.     DECL.   DIRECTION(N to E)
             。       。         。
   300    179.41    -1.86      92.89
   400    180.68    -2.42      93.77
   500    181.95    -2.98      94.64
   600    183.22    -3.54      95.52
   700    184.49    -4.10      96.39
   800    185.77    -4.66      97.26
   900    187.04    -5.22      98.13
  1000    188.32    -5.77      99.00
  1100    189.60    -6.32      99.87
  1200    190.89    -6.87     100.73
  1300    192.18    -7.42     101.59
  1400    193.47    -7.97     102.45
  1500    194.76    -8.51     103.31
  1600    196.06    -9.05     104.16

57