# 状況を再現してみる
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()

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

## 弾道計算関数
 
# 固定値
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)

# 打揚角度 [°]
deg0 = 89.6
# 初速 [km] 
vo = 48400

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

タグ ：

Python

プログラミング

自作自演

卵

弾道計算