弾道ミサイルの軌道（１）からの続き
そこで，各方向の加速度を0.069 km/s^2 (7G)，加速時間を60 sにしたところ，vt = 4.07 km/s，vr = 3.58 km/s，v0 = 5.42 km/s（マッハ15.9）が得られた。このときの，最高高度が1160 km，到達距離が4270 km である。これでよいかと思ったが，調べてみると，そもそも推進剤の質量が全質量に対して非常に大きな割合を占めることがわかった。やり直し。
質量が変化する場合の極座標の運動方程式では，$m \ddot{\bm{r}}\ $の項に$\ \dot{m} \dot{\bm{r}}\ $が加わることになる。これを極座標にすれば，$\dot{m} (\dot{r} \bm{e}_r + r \dot{\theta} \bm{e}_\theta)\ $である。そこで，運動方程式の各成分は次の通り。
$\ddot{r} - \dfrac{h^2}{r^3} + \dfrac{\dot{m}}{m} \dot{r} = \frac{1}{m} F_r = -g \dfrac{R^2}{r^2} + \alpha H(\tau- t)$
$\dfrac{1}{r }(\dot{h} + \dfrac{\dot{m}}{m} h) = \frac{1}{m} F_\theta = 0 + \beta H(\tau-t) $
なお，$H(x) = 1 (x>0) ; =0 (x<0)$はヘヴィサイドの階段関数である。$\alpha, \beta$は，加速開始時刻 $t=0$から加速終了時刻 $t=\tau$ まで動径方向と角度方向に加わる加速度を表わす。
出発前の弾道ミサイルの全質量を$m_0$，全質量に対する推進剤の割合を $p$とすると，加速中（$0 \le t \le \tau$）の単位時間当たりの噴射質量は，$\dfrac{p m_0}{\tau}$となる。そこで，時刻 $t$における弾道ミサイルの質量は，$m(t) = m_0 - \dfrac{p m_0}{\tau} t = m_0 (1 - p\dfrac{t}{\tau})$，$\dot{m}(t) = - m_0 \dfrac{p}{\tau}$。
したがって，$\dfrac{\dot{m}}{m} = - \dfrac{p}{\tau - p t} \quad (0 \le t \le \tau)$ であり，$t>\tau$ では $\ \dfrac{\dot{m}}{m} = 0\ $となる。
なお，最高高度は，$\dot{r}(t)=0$となる時刻 $t_p$に対応する$r(t_p)$で与えられ，到達距離は，$\int_0^T \dfrac{R h(t)}{r(t)^2}dt$で得られる。ただし，$T$は到達時間である。
準備ができたので再計算してみる。全重量に対する推進剤の比率は，p=0.75（3/4）と仮定した。加速時間 60s，加速度 0.0446 km/s^2（4.6 G）で到達時間1320 sが再現される。最高高度は，t=682sのとき 1060 km。加速終了時の速度は，vt = 4.69 km/s，vr = 3.32 km/s，v0 = 5.75 km/s （マッハ16.9）である。また，到達距離は4930 kmであり，6-7%の誤差で報告値に一致した。
g = 0.0098; R = 6350; τ = 60; p = 0.75; a= 0.0446; s = Pi/4fr[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, 1320}]f[t_] := r[t] /. sol[[1, 1]]d[t_] := h[t] /. sol[[1, 2]]Plot[{6350, f[t]}, {t, 0, 1320}]Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, 1319}, PlotRange -> {-4, 5}](f[t] - R) /. FindRoot[D[f[t], t] == 0, {t, 660}]{d[τ]/f[τ], f[τ+1] - f[τ]}{Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2],Sqrt[(d[τ]/f[τ])^2 + (f[τ+1] - f[τ])^2]/.34}{NIntegrate[R d[t]/f[t]^2 , {t, 0, τ}],NIntegrate[R d[t]/f[t]^2 , {t, 0, 1320}]}
図１：軌道半径の時間変化
図２：速度vr，速度vtの大圏射影，速度vt の時間変化
推進時間 60s ，加速度 4.55G，加速角 45° で，高度 1059 km，到達距離 4933 km
推進時間 90s ，加速度 3.22G，加速角 49° で，高度 1038 km，到達距離 4813 km
推進時間120s ，加速度 2.57G，加速角 52.5° で，高度 1015 km，到達距離 4697 km
というわけで，報告値の 1-2%の範囲に収めることも可能な定性的モデルができた。
（付録）最後のパラメタでは，津軽海峡に到達する発射後420 s，距離1368km，高度811kmの速度は3.9 km/s であり，30秒ほどで津軽海峡（幅130 km）を通過することになる。
弾道ミサイルの軌道（３）に続く
