乙Py先生のプログラミング教室
初学者のためのプログラミング学習サイト

弾道計算

とあるビルに卵が投げつけられたそうなので
検証してみることとしました


解説動画はこちら



さて、この事件のように階下の道から卵を投げて
このような状況になるのか検証してみました。

とあるビルを再現するかのような画像が
twitterにあったのでそちらを参考にすると

卵はビルの三階、柵際10cmあたりの所に着弾

現場のビルの状況は
近場を歩いている人の身長を1.6mと仮定すると
3階の柵まで7.6m , 3階の天井が8.9m

3階の柵の高さを1m、柵の厚みを5cmで
設定してみました

ビルを可視化してみるとこんな感じですかね

# 状況を再現してみる
from matplotlib import pyplot as plt
from math import pi, cos, sin,atan,sqrt
%matplotlib inline

ax = plt.gca()

# 壁
plt.hlines(    0   ,       -2 , 2 , color="gray")
plt.hlines(8.9   , 0.1    , 2 , color="gray")
plt.hlines(7.6   , 0.1    , 0.15 , color="gray")
plt.hlines(6.6   , 0.15 , 2 , color="gray")
plt.vlines(2   , 6.6    , 8.9 , color="gray")
plt.vlines(0.1   ,       0 , 7.6 , color="gray")
plt.vlines(0.15 ,     0 , 7.6 , color="gray")

plt.xlim(-2,10)
plt.ylim(-1,10)
plt.show()
download


ここからは卵を投げることを考えてみましょう。

卵は一般的なLサイズ(重量60g , 断面積12.5cm2)
発射場所はビルの壁際の道から投下するものとします
変数に用意しておきましょう

# '発射時高度 [m] : '
ht = 0
# '質量 [g] : '
mas = 60
# '断面積[cm2] : '
area = 12.5


卵の弾道計算には関数を設定します

これは以前ライフルの弾道計算に用いたものです
[Pythonプログラミング]2km先の敵を狙撃するために弾道計算してみよう


## 弾道計算関数
 
# 固定値
dt = 0.1# 時間刻み[s]
gr = 9.80665 # 重力定数 [m/s2]  初期値 9.80665
cd = 0.3 # 空気抵抗係数(Cd)
airRho =  1.22 # 空気密度 [kg/m3]  初期値 1.22 (1atm , 摂氏15, 湿度50%)
w0 = 0 # 風速(地上) [m/s]
w1 = 0 # 風速(毎100m) [m/s]
 
# 微分計算用関数
def dvx( vx, vy, w, k):
    return -k * (vx - w) * sqrt((vx - w) * (vx - w) + vy*vy)
 
def dvy( vx, vy, w, k, g):
    return -g - k * vy * sqrt((vx - w) * (vx - w) + vy*vy)
 
# 弾道計算関数
def ballistic_calc(mas,area,deg0,vo,ht):
    deg = deg0  / 180 * pi # 角度 [ラジアン]
    _S = area / 10000 # 断面積 [m2]
    _K = 0.5 * _S * cd * airRho / mas *1000 # 抵抗力係数
    _T , _X , _Y , _v = 0 , 0 , ht , vo
    _vx , _vy , _W = _v * cos( deg ) , _v * sin( deg ) , w0 + w1 * (_Y / 100)
    _y_max , _t_max = _Y , 0
    k_x1, k_x2, k_x3, k_x4, k_y1, k_y2, k_y3, k_y4 = 0,0,0,0,0,0,0,0
    vx_2 , vy_2 , x_n2 , y_n2 = 0,0,0,0
    result_x , result_y = [0],[0]
    for i in range(1,101):
        k_x1 = dt * dvx( _vx, _vy, _W, _K )
        k_y1 = dt * dvy( _vx, _vy, _W, _K, gr )
        k_x2 = dt * dvx( _vx + k_x1 / 2, _vy + k_y1 / 2, _W, _K )
        k_y2 = dt * dvy( _vx + k_x1 / 2, _vy + k_y1 / 2, _W, _K, gr )
        k_x3 = dt * dvx( _vx + k_x2 / 2, _vy + k_y2 / 2, _W, _K )
        k_y3 = dt * dvy( _vx + k_x2 / 2, _vy + k_y2 / 2, _W, _K, gr )
        k_x4 = dt * dvx( _vx + k_x3, _vy + k_y3, _W, _K )
        k_y4 = dt * dvy( _vx + k_x3, _vy + k_y3, _W, _K, gr )
        vx_2 = _vx + (k_x1 + 2 * k_x2 + 2 * k_x3 + k_x4) / 6
        vy_2 = _vy + (k_y1 + 2 * k_y2 + 2 * k_y3 + k_y4) / 6
        x_n2 = _X + ( _vx + vx_2 ) / 2 * dt
        y_n2 = _Y + ( _vy + vy_2 ) / 2 * dt
        #print('経過秒数 : {0:.2} , 飛距離 : {1:.05}m , 到達高度 : {2:.05}m'.format(i*dt,x_n2 , y_n2) )
        result_x.append(x_n2)
        result_y.append(y_n2)
        if y_n2<0:  
            #print('break')
            break
        _T = i * dt
        _vx , _vy = vx_2 , vy_2
        _X , _Y = x_n2 , y_n2
        _v = sqrt( _vx * _vx + _vy * _vy )
        deg = atan ( _vy / _vx )
        _W = w0 + w1 * (_Y / 100)
    return result_x , result_y
次に描画用の関数を設定します
def egg_throw_viewer(mas,area,deg0,vo,ht,xlmi,xlmx,ylmi,ylmx):
    ax = plt.gca()
    bx,by = ballistic_calc(mas,area,deg0,vo/3600,ht)
    print('経過秒数 : {0:.2} , 飛距離 : {1:.05}m , 到達最高度 : {2:.05}m'.format((len(bx) - 1)*0.1 , bx[-1] , max(by)))
    print('時速 : {0}km , 角度 : {1:.05}°'.format(vo / 1000 , deg0))

    plt.plot(bx , by , linestyle = "dashed")
    plt.hlines(max(by) , 0 ,1,color="green")

    # 壁
    plt.hlines(    0   ,       -2 , 2 , color="gray")
    plt.hlines(8.9   , 0.1    , 2 , color="gray")
    plt.hlines(7.6   , 0.1    , 0.15 , color="gray")
    plt.hlines(6.6   , 0.15 , 2 , color="gray")
    plt.vlines(2   , 6.6    , 8.9 , color="gray")
    plt.vlines(0.1   ,       0 , 7.6 , color="gray")
    plt.vlines(0.15 ,     0 , 7.6 , color="gray")

    plt.xlim(xlmi,xlmx)
    plt.ylim(ylmi,ylmx)
    plt.show()


用意ができたところで
軽く投げてみましょう。
# 打揚角度 [°]
deg0 = 89.5
# 初速 [km] 
vo = 100000

egg_throw_viewer(mas,area,deg0,vo,ht,-1,1,-1,10)

download-1


めちゃめちゃ上の方に向かって投げますが
角度が浅いとしたの壁に
速度が早すぎると天井や上の壁にぶつかります。


成功パターンはこれでした。
# 打揚角度 [°]
deg0 = 89.6
# 初速 [km] 
vo = 48400

egg_throw_viewer(mas,area,deg0,vo,ht,-1,1,-1,10)
download-2

これならギリ成功できます。


まとめ

卵は一般的なLサイズ(重量60g , 断面積12.5cm2)
壁際から時速 : 48.4km , 角度 : 89.6°で投げると
天井にもぶつからず、柵際10cmあたりの場所に
投げ込むことが出来ました。

1投げで成功させたのは神の領域のコントロール
正に「奇跡」を呼ぶ「軌跡」
だったんでしょうねーーーーーーーーー

きっと1000卵ほど投げて練習したんでしょう
きっとスーパーから卵が無くなりますね
そんなスーパーに近くに住んでいる人が
犯人なんでしょうwwwwwwww

真実は一つです。

それでは。

世界情勢がきな臭くなってきました。
いざという時のために
いろいろ出来るようになっておきましょう。

本日は2km先の敵を狙撃するための
弾道計算です。


解説動画はこちら



本日はライフルの弾道計算です

計算方法は非常に複雑ですが
基本的な考えは放物線運動です。

ライフルの弾は水平に構えたとしても
真横には飛ばず、空気抵抗や
重力の影響を受けてどんどん下に落ちて行きます。

そのため、少し角度をつけてやらないと
遠くには飛びません。

計算の際に必要になってくるのは
・弾の質量(重さ)
・弾の断面積
・銃の角度
・弾の初速
・銃を打つ際の高度
になります。

これを可変に出来るように
簡易シミュレーターを作ってみました。

ソースコードはこちら
from matplotlib import pyplot as plt
import ipywidgets as widgets
from IPython.display import display
from ipywidgets import interact, FloatSlider, IntSlider
from math import pi, cos, sin,atan,sqrt
%matplotlib inline

# 固定値
dt = 0.1# 時間刻み[s]
gr = 9.80665 # 重力定数 [m/s2]  初期値 9.80665
cd = 0.3 # 空気抵抗係数(Cd)
airRho =  1.22 # 空気密度 [kg/m3]  初期値 1.22 (1atm , 摂氏15, 湿度50%)
w0 = 0 # 風速(地上) [m/s]
w1 = 0 # 風速(毎100m) [m/s]

# 微分計算用関数
def dvx( vx, vy, w, k):
    return -k * (vx - w) * sqrt((vx - w) * (vx - w) + vy*vy)

def dvy( vx, vy, w, k, g):
    return -g - k * vy * sqrt((vx - w) * (vx - w) + vy*vy)

# 弾道計算関数
def ballistic_calc(mas,area,deg0,vo,ht):
    deg = deg0  / 180 * pi # 角度 [ラジアン]
    _S = area / 10000 # 断面積 [m2]
    _K = 0.5 * _S * cd * airRho / mas *1000 # 抵抗力係数
    _T , _X , _Y , _v = 0 , 0 , ht , vo
    _vx , _vy , _W = _v * cos( deg ) , _v * sin( deg ) , w0 + w1 * (_Y / 100)
    _y_max , _t_max = _Y , 0
    k_x1, k_x2, k_x3, k_x4, k_y1, k_y2, k_y3, k_y4 = 0,0,0,0,0,0,0,0
    vx_2 , vy_2 , x_n2 , y_n2 = 0,0,0,0
    result_x , result_y = [0],[0]
    for i in range(1,101):
        k_x1 = dt * dvx( _vx, _vy, _W, _K )
        k_y1 = dt * dvy( _vx, _vy, _W, _K, gr )
        k_x2 = dt * dvx( _vx + k_x1 / 2, _vy + k_y1 / 2, _W, _K )
        k_y2 = dt * dvy( _vx + k_x1 / 2, _vy + k_y1 / 2, _W, _K, gr )
        k_x3 = dt * dvx( _vx + k_x2 / 2, _vy + k_y2 / 2, _W, _K )
        k_y3 = dt * dvy( _vx + k_x2 / 2, _vy + k_y2 / 2, _W, _K, gr )
        k_x4 = dt * dvx( _vx + k_x3, _vy + k_y3, _W, _K )
        k_y4 = dt * dvy( _vx + k_x3, _vy + k_y3, _W, _K, gr )
        vx_2 = _vx + (k_x1 + 2 * k_x2 + 2 * k_x3 + k_x4) / 6
        vy_2 = _vy + (k_y1 + 2 * k_y2 + 2 * k_y3 + k_y4) / 6
        x_n2 = _X + ( _vx + vx_2 ) / 2 * dt
        y_n2 = _Y + ( _vy + vy_2 ) / 2 * dt
        #print('経過秒数 : {0:.2} , 飛距離 : {1:.05}m , 到達高度 : {2:.05}m'.format(i*dt,x_n2 , y_n2) )
        result_x.append(x_n2)
        result_y.append(y_n2)
        if y_n2<0: 
            #print('break')
            break
        _T = i * dt
        _vx , _vy = vx_2 , vy_2
        _X , _Y = x_n2 , y_n2
        _v = sqrt( _vx * _vx + _vy * _vy )
        deg = atan ( _vy / _vx )
        _W = w0 + w1 * (_Y / 100)
    return result_x , result_y


Widgetsを使ってシミュレーションするのは
こちらの関数です。

mas = widgets.FloatSlider(value=25.4  , min=5  , max=50, step=0.1 , description='質量 [g] : ')
area = widgets.FloatSlider(value=1.19, min=1 , max=5, step=0.01 , description='断面積[cm2] : ')
deg0 = widgets.FloatSlider(value=2.0  , min=0.1 , max=90  , step=0.1 , description='打揚角度 [°] : ')
vo = widgets.IntSlider(value=890, min=500, max=1200 , step=10 , description='初速 [m/s] : ')
ht = widgets.IntSlider(value=0  , min=0 , max=100  , step=1 , description='発射時高度 [m] : ')

@interact(mas=mas,area=area,deg0=deg0,vo=vo,ht=ht)
def on_value_change(mas,area,deg0,vo,ht):
    ax = plt.gca()
    bx,by = ballistic_calc(mas,area,deg0,vo,ht)
    print('経過秒数 : {0:.2} , 飛距離 : {1:.05}m , 到達最高度 : {2:.05}m'.format((len(bx) - 1)*0.1 , bx[-1] , max(by)))
    plt.plot(bx , by , linestyle = "dashed")
    plt.hlines(0 , max(bx),5,color="red")
    plt.hlines(max(by) , max(bx),5,color="green")
    plt.show()

スクリーンショット 2022-03-12 17.51.27

こんな感じで5つの変数を動かせるようにしています。
デフォルトでは一般的な狙撃にしようされる弾丸
7.62x51mm NATO弾の重量 : 25.4g
初速890 [m/s]
断面積1.19cm2 を入れています。


2km先の敵を打つには
結構な角度をつけないと届かないようですね。

弾や角度条件を変えれば
色々な銃の狙撃をシミュレーション出来るので
狙撃の予定がある人は
参考にしてみて下さい。

本日はスナイパーに関する
ライフルの弾道計算のお話でした。

それでは。




 

このページのトップヘ