目次

1.はじめに
 2.円制限三体問題-CR3BP
3.双円制限四体問題-BCR4RP

1. はじめに

地球-月領域では、円制限三体問題や双円制限4体問題などの仮定を用いることで、力学的な平衡点（ラグランジュ点）やその周りの周期軌道を用いた効率的な軌道設計の研究が行われてきた。本記事ではこれらの運動方程式の導出を行います。私は特に正規化の部分で混乱した経験があるので、誰かの助けに慣れば幸いです。

2. 円制限三体問題-CR3BP

2.1. 物理量の正規化

以下では宇宙機は地球及び月の運動に影響を及ぼさないものとし，地球と月は二体の重心を中心とした円運動をしているものと仮定し，月と地球が固定されている座標を考える. まず，地球と月の円運動の角速度，地球と月の質量の和，地球と月の間の距離が全て 1 になるように物理量の正規化を行う. つまり，この座標系においての時間，質量，長さをそれぞれ $\tilde{t}, \tilde{m}, \tilde{r}$ とし，対応する実際の時間，質量，長さをそれぞれ $t_{\mathrm{real}}[\mathrm{s}], m_{\mathrm{real}}[\mathrm{kg}], r_{\text {real }}[\mathrm{m}]$ とすると，各物理量は以下の様に正規化される．

$$
\tilde{m}=\dfrac{m_{\mathrm{real}}}{m_e+m_m},\quad \tilde{t}={\omega_{em}}{t_{real}},\quad \tilde{r}=\dfrac{r_{real}}{r_{em}}.
$$

正規化された座標系における月の質量 $\mu$ は，

$$
\mu=\dfrac{m_m}{m_e+m_m}
$$

‼️: 後ほど出てくるが、月と地球が固定された回転座標系での地球と月の位置は，２天体の重心を原点として，地球 $(-\mu,0)$ , 月 $(1-\mu,0)$ となることに注意する。

2.2 CR3BPの運動方程式(慣性系)

宇宙機の慣性系における（正規化していない）運動方程式は次式で表される．なお，地球と月の重力定数をそれぞれ $\mu_e=Gm_e,\quad\mu_m=Gm_m$ とし，原点は地球と月の重心とする．このとき、地球と月はその重心周りの円軌道を角速度ωemで回るという仮定を用いると、詳しい運動方程式が導出できる（省略）。

$$
\mathbf{\ddot{r}}=\dfrac{\partial U}{\partial \mathbf{r}}
$$

$$
U=-\left(\dfrac{\mu_e}{|\mathbf{r}-\mathbf{r_e}|}+\dfrac{\mu_m}{|\mathbf{r}-\mathbf{r_m}|}\right)
$$

2.3 CR3BPの運動方程式(回転系)

地球と月の重心を原点とし，角速度 $\bf{\omega_{em}}$ で回転する系 $(x, y, z)$ で運動を考える．角速度 $\bf{\omega_{em}}$ は地球と月の重心回りの公転角速度であり，その大きさは以下の式で与えられる．

$$
\omega_{em}=\sqrt{\dfrac{Gm_e+Gm_m}{r_{em}^3}}=\sqrt{\dfrac{\mu_e+\mu_m}{r_{em}^3}}
$$

慣性系での宇宙機の運動方程式の左辺（時間微分） $\mathbf{\ddot{r}}$ を回転系で表現すると，

$$
\bf{\ddot{r}}=\dfrac{D^2\bf{r}}{Dt^2}+2\bf{\omega_{em}}\times\dfrac{D\bf{r}}{Dt}+\bf{\omega_{em}}\times(\bf{\omega_{em}}\times\bf{r})
$$

ここで， $\dfrac{D}{Dt}$ は回転系での微分である．右辺の各項を成分表示すると，

$$
\dfrac{D^2\bf{r}}{Dt^2}=\left(\begin{array}{c}
\ddot{x}\
\ddot{y}\
\ddot{z}
\end{array}\right)
$$

$$
2\bf{\omega_{em}}\times\dfrac{D\mathbf{r}}{Dt}=2\omega_{em}\left(\begin{array}{ll}
-\dot{y}\
\dot{x}\
0
\end{array}\right)
$$

であるため，回転系での運動方程式を成分表示すると以下の様になる．

$$
\ddot{x}-2\omega_{em}\dot{y}-\omega_{em}^2x=-\dfrac{\partial U}{\partial x},
$$

$$
\ddot{y}+2\omega_{em}\dot{x}-\omega_{em}^2y=-\dfrac{\partial U}{\partial y},
$$

$$
\ddot{z}=-\dfrac{\partial U}{\partial z}.
$$

左辺の第三項は，位置による項であり，右辺に回してポテンシャル $U$ と合わせて考えると以下のようになる．

$$
\ddot{x}-2\omega_{em}\dot{y}=-\dfrac{\partial}{\partial x}\left(\dfrac{1}{2}\omega_{em}^2x^2-U\right),
$$

$$
\ddot{y}+2\omega_{em}\dot{x}=-\dfrac{\partial}{\partial y}\left(\dfrac{1}{2}\omega_{em}^2y^2-U\right),
$$

$$
\ddot{z}=-\dfrac{\partial U}{\partial z}.
$$

ここで，以下の擬似ポテンシャル $U^*$ を考える．

$$
U^*\equiv\dfrac{1}{2}\omega_{em}^2(x^2+y^2)-U
$$

$$
=\dfrac{1}{2}\omega_{em}^2(x^2+y^2)+\left(\dfrac{\mu_e}{|\mathbf{r}-\mathbf{r_e}|}+\dfrac{\mu_m}{|\mathbf{r}-\mathbf{r_m}|}\right)
$$

2.4 無次元化したCR3BPの運動方程式(回転系)

それぞれの演算を無次元化すると，( $y,z$ は同様)

$$
\tilde{x} =\frac{x}{r_{em}},
$$

$$
\dot{\tilde{x}} =\frac{dx}{dt}=\frac{1}{\omega_{em}r_{em}}\frac{dx}{dt},
$$

$$
\ddot{\tilde{x}} =\frac{1}{\omega_{em}^2r_{em}}\frac{d^2x}{dt^2}.
$$

以上より，運動方程式を無次元化すると，

$$
\omega_{em}^2r_{em}\ddot{\tilde{x}}-2\omega_{em}^2r_{em}\dot{\tilde{y}}=\dfrac{1}{r_{em}}\dfrac{\partial U^*}{\partial \tilde{x}},
$$

$$
\omega_{em}^2r_{em}\ddot{\tilde{y}}+2\omega_{em}^2r_{em}\dot{\tilde{x}}=\dfrac{1}{r_{em}}\dfrac{\partial U^*}{\partial \tilde{y}},
$$

$$
\omega_{em}^2r_{em}\ddot{\tilde{z}}=\dfrac{1}{r_{em}}\dfrac{\partial U^*}{\partial \tilde{z}},
$$

$$
U^*=\dfrac{1}{2}r_{em}^2\omega_{em}^2(\tilde{x}^2+\tilde{y}^2)+G\dfrac{m_e+m_m}{r_{em}}\dfrac{1-\mu}{\sqrt{(\tilde{x}+\mu)^2+\tilde{y}^2}}+G\dfrac{m_e+m_m}{r_{em}}\dfrac{\mu}{\sqrt{(\tilde{x}-1+\mu)^2+\tilde{y}^2}}
$$

$$
=\dfrac{1}{2}r_{em}^2\omega_{em}^2(\tilde{x}^2+\tilde{y}^2)+r_{em}^2\omega_{em}^2\dfrac{1-\mu}{\sqrt{(\tilde{x}+\mu)^2+\tilde{y}^2}}+r_{em}^2\omega_{em}^2\dfrac{\mu}{\sqrt{(\tilde{x}-1+\mu)^2+\tilde{y}^2}}
$$

最後に、運動方程式の両辺を$r_{em}\omega_{em}^2$で割ると、以下のように正規化したCR3BPの運動方程式を求めることができる。

$$
\ddot{\tilde{x}}-2\dot{\tilde{y}}=\dfrac{\partial U^*}{\partial \tilde{x}},
$$

$$
\ddot{\tilde{y}}+2\dot{\tilde{x}}=\dfrac{\partial U^*}{\partial \tilde{y}},
$$

$$
\ddot{\tilde{z}}=\dfrac{\partial U^*}{\partial \tilde{z}},
$$

$$
U^*=\dfrac{1}{2}(\tilde{x}^2+\tilde{y}^2)+\dfrac{1-\mu}{\sqrt{(\tilde{x}+\mu)^2+\tilde{y}^2}}+\dfrac{\mu}{\sqrt{(\tilde{x}-1+\mu)^2+\tilde{y}^2}}
$$

☕コーヒーブレイク（なぜ回転座標系で表すか？）

慣性座標系の運動方程式では、地球月重心、地球中心、月中心のいずれかが用いられるが、いずれの場合でも惑星の位置が時変となります。回転座標系にすることで地球と月がx軸に固定されるため、ラグランジュ点も固定され、様々な軌道設計が行いやすくなります。（Tube Dynamics, Rove Dynamics などの研究があります）

⚠️注意

質量はkg、時間はsを利用することが一般的ですが、距離はm派とkm派がいます。混合してしまわないように注意しましょう。また論文を読むときも、どちらか確認しましょう。

3. 双円制限四体問題-BCR4RP

本記事では、月と地球がその重心B1を中心に円軌道を描くとする。さらに、原点はB1と太陽の重心（系の重心と同じ）とする。また、導出ではなく結果のみの記載になっています。時間があれば導出も書きます....

3.1 物理量の正規化

長さ：Sun-B1 の距離 $A_S$
時間：B1の角速度 $\omega_S$
質量：系の質量 $m_s+m_e+m_m$

で正規化する。

つまり，この座標系においての時間，質量，長さをそれぞれ $\tilde{t}, \tilde{m}, \tilde{r}$ とし，対応する実際の時間，質量，長さをそれぞれ $t_{\mathrm{real}}[\mathrm{s}], m_{\mathrm{real}}[\mathrm{kg}], r_{\text {real }}[\mathrm{m}]$ とすると，各物理量は以下の様に正規化される．

$$
\tilde{m}=\dfrac{m_{\mathrm{real}}}{m_s+m_e+m_m},\quad \tilde{t}={\omega_{S}}{t_{real}},\quad \tilde{r}=\dfrac{r_{real}}{r_{Sun-B1}}.
$$

正規化行った宇宙機の位置は $q_{\mathrm{sc}} = (x, y)^{\top}$ とする。太陽の位置は $q_{\mathrm{s}} = (-\mu_{\mathrm{s}}, 0)^\top$ ，地球と月の重心B1は, $q_{\mathrm{s}} = (1-\mu_{\mathrm{s}}, 0)^\top$ となり，地球と月の位置はそれぞれ以下のように表される．

\begin{align}
{\boldsymbol{q}}_{\mathrm{e}} & = \begin{bmatrix}1-\mu_{\mathrm{s}}-r_{\mathrm{em}} \mu_{\mathrm{m}} \cos \left({\theta}_{\mathrm{m}}\right)\\ \qquad \quad  -r_{\mathrm{em}} \mu_{\mathrm{m}} \sin \left({\theta}_{\mathrm{m}}\right)\end{bmatrix} \\
{\boldsymbol{q}}_{\mathrm{m}} & =\begin{bmatrix}1-\mu_{\mathrm{s}}+r_{\mathrm{em}}\left(1-\mu_{\mathrm{m}}\right) \cos \left({\theta}_{\mathrm{m}}\right)\\ \qquad \quad \quad r_{\mathrm{em}}\left(1-\mu_{\mathrm{m}}\right) \sin \left({\theta}_{\mathrm{m}}\right)\end{bmatrix}
\end{align}

ここで，正規化した各パラメータは以下のような値を用いた¹

$\mu_{\mathrm{S}}=\left(m_{\mathrm{E}}+m_{\mathrm{M}}\right) /\left(m_{\mathrm{S}}+m_{\mathrm{E}}+m_{\mathrm{M}}\right)=3.02319 \times 10^{-6}$ nondim
$r_{\mathrm{em}}= R_{em}/A_S = 3.844e8(m)/1.4968e11(m)=2.56955\times10^{-3}$ nondim
$\omega_S = 1.99640e-7$ rad/s
$\omega_{em}=2.66498e-6$ rad/s
$\omega_{\mathrm{m}}=\omega_{em}/\omega_{S}=13.3489$
$\theta_{\mathrm{m}}(t)=\left(\omega_{\mathrm{m}}-1\right) {t}+{\theta}_{\mathrm{m} 0}(0)$ nondim

3.1 Sun-B1の回転座標系で見たBCR4BPの運動方程式

以上より， $S-B_{1}$ 回転座標系での宇宙機の運動方程式は，

\begin{align}\left\{
\begin{array}{l}
    \ddot{x}= 2\dot{y} + \gamma^*_x  \\
    \ddot{y}= -2\dot{x} + \gamma^*_y
\end{array}
\right.
\end{align}

と表される．ここで $\mu_m = m_m/(m_e+m_m)=0.0122$ とすると， $\gamma^*$ は $S-B_1$ 回転座標系で以下のように定義される擬似ポテンシャルである．

\begin{align}
    \gamma^* = \dfrac{1}{2}(x^2+y^2)+\dfrac{1-\mu_s}{|\boldsymbol{q}_s - \boldsymbol{q}_{sc}|} + \dfrac{\mu_s(1-\mu_m)}{|\boldsymbol{q}_e - \boldsymbol{q}_{sc}|} + \dfrac{\mu_s\mu_m}{|\boldsymbol{q}_m - \boldsymbol{q}_{sc}|}
\end{align}

エネルギーは時変となることに注意する．なお月の重力が無視できる $(\mu_{\mathrm{m}}=0)$ なら，この運動方程式は，太陽-地球円制限三体問題に一致する．

Kaori Onozaki, Hiroaki Yoshimura, and Shane D Ross. Tube dynamics and low energy earth– moon transfers in the 4-body system. Advances in Space Research, Vol. 60, No. 10, pp. 2117–2132, 2017. ↩