マイコン宇宙講座-月食シミュレーション

月食は月が地球の影に入る現象のことです。この地球の影に入るということはどういうことでしょうか。月は地球の周りを公転しています。そのとき、月、地球、太陽というように一直線に並んだときに、太陽の光に照らされた地球には影ができます。その影に入ると、月食という現象が起きます。必ずしも毎回起きる現象ではありません。比較的レアな現象といってよいでしょう。ですから月食が起きると、ニュースになったりします。ただし、月食は起きるものの日本からは見えないこともあります。

部分月食

月食には2つのパターンがあり、部分月食、皆既月食があります。部分月食は月が影の中に入ることはなく、部分的に入る現象のことです。それに対して皆既月食は、月が完全に影の中に入ってしまう現象です。面白いのは、月の99%が影の中に入ったとしても、残りの1%が入っていませんので部分月食となります。なかなかシビアですね。

日本で見ることができる2022年の月食は、11月8日の皆既月食のみです。それ以降も部分日食や皆既月食があります。国立天文台天文情報センター暦計算室のサイトに日月食データーベースがありますので確認してみるとよいでしょう。

サブルーチン sun

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

# SUN
def sun(ta, c, g, j, l, m, n, v):
    # SUN-THETA
    sx =    - 1.000e-5 * ta * math.sin(g + l)
    sx = sx - 1.000e-5 * math.cos(g - l - j)
    sx = sx - 1.400e-5 * math.sin(2 * g - l)
    sx = sx - 3.000e-5 * ta * math.sin(g - l)
    sx = sx - 3.900e-5 * math.sin(n - l)
    sx = sx - 4.000e-5 * math.cos(l)
    sx = sx + 4.200e-5 * math.sin(2 * g + l)
    sx = sx - 2.080e-4 * ta * math.sin(l)
    sx = sx + 3.334e-3 * math.sin(g + l)
    sx = sx + 9.999e-3 * math.sin(g - l)
    sx = sx + 0.39793 * math.sin(l)

    # SUN-RHO
    sy =      2.7e-5 * math.sin(2 * g - 2 * v)
    sy = sy - 3.3e-5 * math.sin(g - j)
    sy = sy + 8.4e-5 * ta * math.cos(g)
    sy = sy - 1.4e-4 * math.cos(2 * g)
    sy = sy - 0.033503 * math.cos(g)
    sy = sy + 1.000421

    # SUN-PHI
    sz =     -1.700e-5 * math.cos(2 * g - 2 * v)
    sz = sz - 1.900e-5 * math.sin(g - v)
    sz = sz + 2.400e-5 * math.sin(4 * g - 8 * m + 3 * j)
    sz = sz - 2.500e-5 * math.cos(g - j)
    sz = sz + 3.000e-5 * math.sin(c - l)
    sz = sz + 4.600e-5 * ta * math.sin(2 * l)
    sz = sz + 6.801e-5 * math.sin(2 * g)
    sz = sz - 7.900e-5 * math.sin(n)
    sz = sz - 8.000e-5 * ta * math.sin(g)
    sz = sz - 9.500e-5
    sz = sz - 3.460e-4 * math.sin(g + 2 * l)
    sz = sz - 1.038e-3 * math.sin(g - 2 * l)
    sz = sz + 0.032116 * math.sin(g)
    sz = sz - 0.041295 * math.sin(2 * l)

    return sx, sy, sz

メインルーチン m59.py

# m59.py
# マイコン宇宙講座
# 5-9 月食シミュレーション・プログラム
from PIL import Image, ImageDraw, ImageFont
from tkinter.constants import SOLID
import tkinter as tk
import math
import lib


M2PI = 2.0 * math.pi

day = [0 for i in range(20)]
hour = [0 for i in range(20)]
minute = [0 for i in range(20)]
xpos = [0 for i in range(20)]
ypos = [0 for i in range(20)]
shadowsize = [0 for i in range(20)]
moonsize = [0 for i in range(20)]

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 = 30.0

i = 0
ii = 0
yz = 0
mz = 0

print()
print('  Date ( JST )                       - X -     - Y -       本影   月の大きさ')
print(' ---------------------------------------------------------------------------')
print('       年   月   日    時   分          。       。          。       ,')

while True:
    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)

    mr = ra
    md = dc
    rl = dl

    sx, sy, sz = lib.sun(ta, c, g, j, l, m, n, v)

    r0 = l
    r1 = sx
    r2 = sy
    r3 = sz

    ra, dc, dl = lib.positions(r0, r1, r2, r3)

    ra += math.pi
    dc = -(dc)
    if ra > M2PI:
        ra -= M2PI

    # PRINT POSITIONS
    x = -(mr - ra) * lib.K[3] * 60.0
    y = (md - dc) * lib.K[3] * 60.0
    rl = 6378.16 * rl
    rr = 0.76 * math.atan(6378.16 / rl) * lib.K[3] * 60.0
    ll = math.atan(1738.0 / rl) * lib.K[3] * 60.0

    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)) + 0.1
    hh = int(hh)
    if int(ms + 0.1) == 60:
        hh = hh + 1
        if hh == 24:
            hh = 0
        ms = 0

    if (x ** 2 + y ** 2) > 16000.0 and (ii == 0):
        jd += (iv / 1440.0)
        continue

    if (x ** 2 + y ** 2) > 16000.0:
        break
    else:
        while True:
            if int(ii / 2) != ii / 2:
                break

            i = i + 1
            day[i] = dd
            hour[i] = hh
            minute[i] = ms
            xpos[i] = x
            ypos[i] = y
            shadowsize[i] = rr
            moonsize[i] = ll
            break

        sfmt = '                 %2d    %2d   %2d        %6.1f   %6.1f       %5.1f    %5.1f' % (dd, hh, ms, x, y, rr, ll)
        if yy != yz:
            sfmt = '   %4d    %2d   %2d    %2d   %2d        %6.1f   %6.1f       %5.1f    %5.1f' % (yy, mm, dd, hh, ms, x, y, rr, ll)
        else:
            if mm != mz:
                sfmt = '           %2d   %2d    %2d   %2d        %6.1f   %6.1f       %5.1f    %5.1f' % (mm, dd, hh, ms, x, y, rr, ll)
        print(sfmt)

        yz = yy
        mz = mm
        ii += 1

        if ii < 18:
            jd += iv / 1440.0
            continue
    #break
print(' ---------------------------------------------------------------------------')
print()
std = input('PRESS ANY KEY')

# GESYOKU PLOT
root = tk.Tk()
root.resizable(False, False)
root.geometry('642x440')
root.title('マイコン宇宙講座 - 月食シミュレーションプログラム')

img = Image.new('RGB', (640, 400), (0, 0, 0))
draw = ImageDraw.Draw(img)

canvas = tk.Canvas(root, width=640, height=400, bg='black')
canvas.pack(anchor=tk.NW)

# Color
moon_color = (240, 210, 0)
shadow_color = (107, 2, 2)

# フォント名は実行環境に合わせて変更すること
# TrueTypeの等幅フォント名を指定する
fonts = ImageFont.truetype('TakaoPGothic.ttf', 15)

if ypos[1] < 0:
    sfmt = '%4d 年  %2d 月  %2d 日の月食' % (yy, mm, day[1])
    draw.text((214, 16), '本影の中をとおる月の姿', font=fonts)
    draw.text((200, 32), sfmt, font=fonts)
else:
    sfmt = '%4d 年  %2d 月  %2d 日の月食' % (yy, mm, day[1])
    draw.text((214, 352), '本影の中をとおる月の姿', font=fonts)
    draw.text((200, 368), sfmt, font=fonts)

rr = shadowsize[1]
for m in range(360):
    x = rr * math.cos(m / lib.K[3]) / 1.5
    y = rr * math.sin(m / lib.K[3]) / 1.5
    x = int(80 + x) * 4
    y = int(50 + y) * 4
    draw.rectangle((x - 2, y - 2, x + 2, y + 2), fill=shadow_color)

for i in range(1, 10):
    if xpos[i] == 0 and ypos[i] == 0:
        break
    xpos[i] = xpos[i] / 1.5 + 80
    ypos[i] = 50 - ypos[i] / 1.5

    x = xpos[i] / 2 - 2
    y = ypos[i] / 4

    if x < 0 or x > 75:
        continue
    if y < 0 or y > 24:
        continue

    sfmt = '%2d時' % (hour[i])
    draw.text((x * 8, y * 16), sfmt, font=fonts)
    if abs(minute[i]) > 1:
        sfmt = '%2d分' % (minute[i])
        draw.text((x * 8, y * 16 + 16), sfmt, font=fonts)

    for m in range(0, 360, 4):
        x = moonsize[i] * math.cos(m / lib.K[3]) / 1.5
        y = moonsize[i] * math.sin(m / lib.K[3]) / 1.5
        x = int(xpos[i] + x)
        y = int(ypos[i] + y)

        if x < 0 or x > 159:
            continue
        if y < 0 or y > 99:
            continue

        x = x * 4
        y = y * 4
        draw.rectangle((x - 2, y - 2, x + 2, y + 2), fill=moon_color)

# 描画した月食の画像を保存
img.save('./images/gesyoku.png')

# 月食表示
photo = tk.PhotoImage(file='./images/gesyoku.png')
canvas.create_image(1, 1, image=photo, anchor=tk.NW)

button = tk.Button(root, text='閉じる', width=6, relief=SOLID, cursor='hand1', command=root.destroy)
button.place(x=558, y=407)

root.mainloop()

例題 1979年9月6日の月食を表示してみよう。時刻は17時とします。

DATE AND TIME(JST) ? 19790906,170000

  Date ( JST )                       - X -     - Y -       本影   月の大きさ
 ---------------------------------------------------------------------------
       年   月   日    時   分          。       。          。       ,
   1979     9    6    17    0          93.2    -55.5        46.7     16.7
            9    6    17   30          75.9    -50.3        46.7     16.7
            9    6    18    0          58.7    -45.2        46.7     16.7
            9    6    18   30          41.4    -40.0        46.7     16.7
            9    6    19    0          24.1    -34.8        46.7     16.7
            9    6    19   30           6.9    -29.6        46.7     16.7
            9    6    20    0         -10.4    -24.4        46.7     16.7
            9    6    20   30         -27.6    -19.2        46.7     16.7
            9    6    21    0         -44.8    -14.0        46.6     16.7
            9    6    21   30         -62.0     -8.7        46.6     16.7
            9    6    22    0         -79.3     -3.5        46.6     16.7
            9    6    22   30         -96.5      1.8        46.6     16.7
            9    6    23    0        -113.7      7.1        46.6     16.7
 ---------------------------------------------------------------------------
1979年9月6日の月食の様子
1979年9月6日の月食の様子