弾道計算
参考数値:
第一宇宙速度 (地表すれすれを円軌道で運動) 約 7.9 km/s (= 28,400 km/h)
第二宇宙速度 (地球脱出速度) 約 11.2 km/s（40,300 km/h）
第三宇宙速度 (太陽系脱出速度) 約 16.7 km/s (60,100 km/h)

Download script from .\cannonball_round.py

import sys
import numpy as np
from numpy import sqrt, exp, sin, cos, tan, arctan, pi
from matplotlib import pyplot as plt


"""
弾道計算
参考数値:
第一宇宙速度 (地表すれすれを円軌道で運動) 約 7.9 km/s (= 28,400 km/h)
第二宇宙速度 (地球脱出速度) 約 11.2 km/s（40,300 km/h）
第三宇宙速度 (太陽系脱出速度) 約 16.7 km/s (60,100 km/h)
"""


#===================================
# 定数
#===================================
# 円周率
pi = 3.14159265358979323846

# 万有引力定数
G = 6.67430e-11 # m^3/kg/s^2

# 太陽の質量
Ms = 1.99e30 # kg

# 地球の質量、半径、公転半径
Me = 5.972e24 # kg
Re = 12756.274e3 * 0.5 # m
Rrev = 1.50e11 # m

# 重力加速度 g = G*Me/Re^2
g = 9.8 # m/s^2

# 砲弾 (cannonball) の質量 計算結果には影響しない
mc = 1.0 # kg

# 射出速度、角度
v0 = 100.0 # m/s
Q0 = 45.0 # 度

# 計算時間間隔
dt = 1

# 最大ステップ数
nmax = 100


#===================================
# グラフ設定
#===================================
fontsize = 20
# グラフの更新時間
tsleep = 0.1 # s


#============================
# 起動時引数で初期値を更新
#============================
argv = sys.argv
if len(argv) >= 2:
    v0 = float(argv[1])
if len(argv) >= 3:
    Q0 = float(argv[2])
if len(argv) >= 4:
    dt = float(argv[3])
if len(argv) >= 5:
    nmax = int(argv[4])
if len(argv) >= 6:
    tsleep = float(argv[5])


def usage():
    print("")
    print("使い方: python cannonball_round.py 初速度(m/s) 射出角度(度) 時間間隔(s) 繰り返し回数 描画更新時間(s)")
    print("例1: python cannonball_round.py")
    print(" 初期値 v0 = 100.0 m/s, Q0 = 45度, dt = 1 sで計算")
    print("例2: python cannonball_flat.py 10000 20 10")
    print(" v0 = 10000.0 m/s, Q0 = 20度, dt = 10 sで計算")
    print("")


# 質点に係る力のベクトル Fx, Fy を返す
def force(x, y):
    global g, mc

    r = sqrt(x * x + y * y)
    f = -mc * g * (Re / r)**2
    return f * x / r, f * y / r

def main():
    vx0 = v0 * cos(Q0 * pi / 180.0)
    vy0 = v0 * sin(Q0 * pi / 180.0)
    print("g(R0={:12.4g} m)={:10.4f} m/s^2".format(Re, G * Me / Re / Re))
    print("g={:10.4f} m/s^2".format(g))
    print("vx0 = {:10.4f} m/s".format(vx0))
    print("vy0 = {:10.4f} m/s".format(vy0))
    print("dt = {} s".format(dt))

# 宇宙速度の計算
# 地表での自転速度
    vrot = 2.0 * pi * Re / 24.0 / 60.0 / 60.0
# 地球の公転速度
    vrev = 2.0 * pi * Rrev / 365.24 / 24.0 / 60.0 / 60.0
# 第一宇宙速度: 遠心力と重力のつり合い: mc * v1^2 / Re = G * Me * mc / Re^2
    v1 = sqrt(G * Me / Re)
# 第二宇宙速度: EU + 0.5 * mc * v2 * v2 >= 0: v2 = sqrt(2.0 * G * Me / Re) = sqrt(2.0) * v1
    v2 = sqrt(2.0 * G * Me / Re)
# 地球の公転からの脱出速度: 太陽と地球の重力より初期運動エネルギーが大きい
#     -G * Ms * mc / Rrev + 0.5 * mc * (vrev3 - vrev)^2 >= 0
#     vrev3 = -vrev + sqrt(2.0 * G * Ms / Rrev)
    vrev3 = sqrt(2.0 * G * Ms / Rrev) - vrev
# 第三宇宙速度: 地球の重力を振り払った後にvrev3になる速度
#     -G * Me * mc / Re + 0.5 * mc * v3^2 >= 0.5 * mc * vrev3^2
#     v3 = sqrt(2.0 * G * Ms / Rrev + vrev3^2)
    v3 = sqrt(2.0 * G * Me / Re + vrev3 * vrev3)
    print("")
    print("宇宙速度")
    print("　地表における自転速度: vrot = {:14.4g} m/s".format(vrot))
    print("　地球の公転速度　　　: vrev = {:14.4g} m/s".format(vrev))
    print("　第一宇宙速度　　　　: v1 = {:14.4g} m/s".format(v1))
    print("　第二宇宙速度　　　　: v2 = {:14.4g} m/s".format(v2))
    print("　公転脱出速度　　　　: vrev3 = {:14.4g} m/s".format(vrev3))
    print("　第三宇宙速度　　　　: v3 = {:14.4g} m/s".format(v3))

# 初期状態におけるエネルギー保存則関係の計算
# 重力のポテンシャルエネルギー
    EU = -G * Me * mc / Re
# 運動エネルギー
    EK = 0.5 * mc * v0 * v0
# 全エネルギー
    Etot = EU + EK
# 最大到達高度 hmax: Etot = -G * Me * mc / (Re + hmax) => hmax = -Re - G * Me * mc / Etot
    hmax = -Re - G * Me * mc / Etot
    print("")
    print("初期状態")
    print("　重力のポテンシャルエネルギー: EU = {:14.4g} J".format(EU))
    print("　運動エネルギー　　　　　　　: EK = {:14.4g} J".format(EK))
    print("　全エネルギー　　　　　　　　: Etot= {:14.4g} J".format(Etot))
    print("最大到達高度: hmax = {:14.4g} m".format(hmax))

# 初期値を配列に設定
    t = [0.0]
    vx = [vx0]
    vy = [vy0]
    x = [0.0]
    y = [Re]

# グラフの作成
    fig = plt.figure(figsize = (8, 8))
    ax1 = fig.add_subplot(1, 1, 1)
    ax1.set_xlabel("x (m)", fontsize = fontsize)
    ax1.set_ylabel("y (m)", fontsize = fontsize)
    ax1.set_aspect('equal')
    ax1.tick_params(labelsize = fontsize)
    plotdata, = ax1.plot(x, y, linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 2.0)
# 地表を描画
    Qg = range(361)
    xg = [Re * sin(iQg * pi / 180.0) for iQg in Qg]
    yg = [Re * cos(iQg * pi / 180.0) for iQg in Qg]
    ax1.plot(xg, yg, linestyle = '-', linewidth = 0.5, color = 'red')

# 計算ループ開始
    print("")
    print("{:8}: {:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}\t{:10}"
                .format('t(s)', 'vx(m/s)', 'vy(m/s)', 'x(m)', 'y(m)', 'Ekin(J)', 'EU(J)', 'Etot(J)'))
    for i in range(1, nmax):
        tt = i * dt
# m(vx[i] - vx[i-1]) / dt = 0
# m(vy[i] - vy[i-1]) / dt = -mg
        fx, fy = force(x[i-1], y[i-1])
        vxt = vx[i-1] + dt * fx / mc
        vyt = vy[i-1] + dt * fy / mc
# (x[i] - x[i-1]) / dt = vx[i-1]
# (y[i] - y[i-1]) / dt = vy[i-1]
        xt = x[i-1] + dt * vx[i-1]
        yt = y[i-1] + dt * vy[i-1]

        vx.append(vxt)
        vy.append(vyt)
        x.append(xt)
        y.append(yt)
# エネルギー保存則の確認: 計算誤差があるので、計算結果では全エネルギーは一定にならない
        Ekin = 0.5 * mc * (vxt * vxt + vyt * vyt)
        EU = -G * Me * mc / sqrt(xt*xt + yt*yt)
        Etot = Ekin + EU
        print("{:8g}: {:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}"
                .format(tt, vxt, vyt, xt, yt, Ekin, EU, Etot))

# グラフの表示データを更新
        plotdata.set_data(x, y)
# グラフ表示を更新
        ax1.set_xlim([min(x), max(x)])
        ax1.set_ylim([min(y), max(y)])
        plt.tight_layout()
        plt.pause(tsleep)

# 地表に激突したら計算ループ終了
        r = sqrt(xt * xt + yt * yt)
        if r < Re:
            break


    print("プログラムを終了するにはENTERキーを押してください")
    input()

    usage()
    print("")


if __name__ == "__main__":
    main()