弾道計算

Download script from .\cannonball_flat.py

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


"""
弾道計算
"""


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

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

# 地球の質量、半径
Me = 5.972e24 # kg
Re = 12756.274e3 * 0.5 # 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_flat.py 初速度(m/s) 射出角度(度) 時間間隔(s) 繰り返し回数 描画更新時間(s)")
    print("例1: python cannonball_flat.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

    return 0.0, -mc * g

def main():
    vx0 = v0 * cos(Q0 * pi / 180.0)
    vy0 = v0 * cos(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))

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

# グラフの作成
    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.tick_params(labelsize = fontsize)
    plotdata, = ax1.plot(x, y, linestyle = '-', linewidth = 0.5, marker = 'o', markersize = 2.0)
# 地表を描画
    xg = [0.0, Re]
    yg = [0.0, 0.0]
    ax1.plot(xg, yg, linestyle = '-', linewidth = 0.5, color = 'red')

# 計算ループ開始
    print("{:8}: {:10}\t{:10}\t{:10}\t{:10}".format('t(s)', 'vx(m/s)', 'vy(m/s)', 'x(m)', 'y(m)'))
    for i in range(1, nmax):
        tt = i * dt
# m(vx[i] - vx[i-1]) / dt = Fx = 0
# m(vy[i] - vy[i-1]) / dt = Fy = -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)
        print("{:8g}: {:10.4g}\t{:10.4g}\t{:10.4g}\t{:10.4g}".format(tt, vxt, vyt, xt, yt))

# グラフの表示データを更新
        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)

# 地表に激突したら計算ループ終了
        if yt <= 0.0:
            break


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

    usage()

if __name__ == "__main__":
    main()