弾道ミサイルの軌道（６）

弾道ミサイルの軌道（４）からの続き

H3ロケットの発射がフェイルセーフで中止された次の日，なんだか間の悪いことに北朝鮮のICBM級ミサイルが平壌近傍から発射され，北海道渡島大島西200kmに落下した。

飛行時間は66分，飛行距離は900km，最高高度は5700km，落下速度はマッハ20程度と推測される。ペイロードを調整すれば，射程が14,000kmでアメリカ本土が射程に入ると報道された。

それでは，以前作ったのプログラムで確認してみよう。一生懸命Juliaのフォルダを探していたが，Mathematicaのコードだった。昨年11月のコードのパラメータで，飛行時間と飛行距離を微調整した結果が次の通り。

g = 0.0098; R = 6350; τ= 86; p = 0.75; 

a = 0.0446; s = 86.5 Degree; T = 3960; 

(* g = 0.0098; R = 6350; τ = 87; p = 0.75; a

 = 0.0446; s = 45 Degree; T = 5100; *)

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - 

 g R^2/r[t]^2 + fr[t, τ], r[0] == R, r'[0] == 0, 

 h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h},

 {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, 

{t, 0, T}, PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}], 

  f[tt] - R}, {tt, 0, T}]

燃焼時間τが86秒，燃料重量比が0.75，燃焼加速度が 0.0446km/s^2，投射角が86.5度である。これで飛行時間を与えると，飛行距離と最高高度と落下時速度が概ね再現できる。万有引力は距離の2乗に反比例し，コリオリ力や空気抵抗は無視，地球の形を考慮して2次元極座標で質量が変化する1段ロケットの微分方程式を解いている。

燃焼時間τを87秒にして，投射角を45度にすれば，飛行時間が5100秒で飛行距離は14,600km，落下時速度はマッハ24になる。このときペイロードはほとんど変化させていない。燃料重量比は同じで燃焼時間を1秒のばしただけだ。

図：投射角86度のロフテッド軌道と投射角45度の弾道軌道