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

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

10月4日の弾道ミサイルは，本来グアムを標的とする射程5,000kmの火星12型相当の中距離弾道弾なので，日本を対象とした射程1,500kmのノドン改良型とは違う。が，面倒なので，前回求めた解の推進加速度を7割程度に落とし加速角度を30°にして，日本を狙ったときの到達時間が適当な値となる解があるかを確認してみた。

g = 0.0098; R = 6350; τ = 60; p = 0.75; a = 0.032; s =  30 Degree;

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, 360}]

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

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

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

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, 299},  PlotRange -> {-4, 7}]

(f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 200}]

　86.6527

FindRoot[D[f[t], t] == 0, {t, 300}]

　　{t -> 199.993}

{d[τ]/f[τ], f[τ+1] - f[τ]}

　　{4.14703, 0.971161}

{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2], 

 Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}

　　{4.25923, 12.5271}

{NIntegrate[R d[t]/f[t]^2 , {t, 0, 200}], 

 NIntegrate[R d[t]/f[t]^2 , {t, 0, 360}]}

　　{651.9, 1305.61}

ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, T}],  f[T] - R}, {T, 0, 360}]

図：弾道ミサイルの軌道（横軸 水平距離 km，縦軸 高度 km）

60秒加速後の速度がマッハ 12.5となり，360秒（これは初期値として与えた）で1300 km離れた日本に到達する。最高高度は90 km弱となる。このあたりは推進加速度とその角度を調整すればなんとでもなる。いずれにせよ，カップ麺ができる3分以内に対応する必要があるのだけれど，Ｊアラートにインプットする情報分析や対応措置は時間内に間に合うのだろうか。