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}]
0 件のコメント:
コメントを投稿