虚数が実数から現れる時（惑星の軌道編）

　[前回](https://qiita.com/hibit/items/2735ea92236dcb8bebac)に引き続き、Pythonを用いて、実数のみを用いた簡単な関係式をシミュレートすると不思議な結果が出てくること及び、虚数を使うとそれがうまく説明できるようなケースを紹介することで、より皆様に虚数に親しみを持ってもらおうと思います。

# 惑星の運動方程式

　突然ですが、万有引力による惑星の運動方程式は以下のようなものになります。

```math 
\displaylines{
\vec{F} = m\vec{a} = -\frac{GMm}{r^3}\vec{r} \\
\vec{F}:力 \\
m:惑星の質量 \\
\vec{a}:惑星の加速度 \\
G:万有引力定数 \\
M:恒星の質量 \\
m:惑星の質量 \\
\vec{r}:惑星の位置 \\
r:\vec{r}の長さ}
```
　解析的に解くには2回積分をする必要があり、かなりめんどくさいのですが、そこはPythonに任せて力任せにシミュレートすればなんともないです。また、仮定は可能な限り単純化します。$GM=1$となるような都合のいい恒星（質点）を考えます。すると以下のようになります。

```math 
\displaylines{
m\vec{a} = -m\frac{\vec{r}}{r^3} \\
\vec{a} = \frac{\vec{r}}{r^3} \\}
```

　これをプログラムでシミュレートしていきます。

## 自由落下

　初期条件として、$\vec{r}=(100,0),\vec{v}=(0,0)$を考えます。恒星の中心から100mの位置にある惑星が、恒星に向かって自由落下していくようなイメージです。ソースコードは以下のようになります（import等は省略）。

```python 
size = 1500

x, y, vx, vy, ax, ay = np.zeros((6, size))
x[0] = 100
y[0] = 0
vx[0] = 0.0
vy[0] = 0.0

for i in range(size-1):
    r = np.linalg.norm((x[i], y[i]))
    if r < 10:
        print("too close!")
        break
    ax[i] = -x[i]/r**3
    ay[i] = -y[i]/r**3
    vx[i+1] = vx[i] + ax[i]
    vy[i+1] = vy[i] + ay[i]
    x[i+1] = x[i] + vx[i]
    y[i+1] = y[i] + vy[i]


plt.figure(figsize=(10, 10))
plt.scatter(x, y, s=10, c=range(size), cmap=cm.coolwarm)
plt.xlim(-400, 400)
plt.ylim(-400, 400)
plt.grid()
plt.show()
```

　`x,y`,`vx,vy`,`ax,ay`がそれぞれ座標、速さ、加速度を表していて、運動方程式をシミュレートしています。`r`が0に近すぎる場合、計算結果がおかしくなってしまうので（現実世界でも、その前に地表に激突したり、電子の反発力に阻まれたりして、0に近い距離には近づけません）、ある閾値より小さくなった計算を打ち切ることにします。

　結果は以下のようになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/76dcb02b-c529-3dc7-5964-26bee2030435.png)

　青がスタート位置、赤が終点になります。

## 円軌道

　次に、グラフで言うと上向きに方向に初速度を加えます。うまく値を設定すると、真円の軌道に乗せることができます。`vy[0]=0.1`がそのような値となります。（`size`もうまく描画範囲になるように調整しています）

```python
x, y, vx, vy, ax, ay = np.zeros((6, size))
x[0] = 100
y[0] = 0
vx[0] = 0.0
vy[0] = 0.1 # ここの部分
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/9030ed8e-998c-d29f-fb13-fdf8ab55d370.png)

　ちなみに、`vy[0] = 0.125`だと以下のような楕円になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/f767b118-ce5b-2d57-5064-72c0f9ba5b98.png)

## 円でなくなる

　しかし、ある速度を超えると、円軌道ではなくなってしまいます。

　`vy[0] = 0.1414`だと以下のような感じになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/5f5c0a3c-9a6a-60b6-bae6-64f86103cd04.png)

　ネタバレをするとこれは放物線になります。

　`vy[0] = 0.2`だと以下のような感じになります。（軌道の特徴を見やすくするため、グラフの中心と尺度を変えています）

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/f240dc88-21fb-78bc-b5f2-813514b87343.png)

　明らかに双曲線ですね。

　現実の世界の天体も、（おおまかに近似すれば）このように楕円軌道（惑星、衛星）や双曲線運動（彗星等はこれの場合もあります）を行うことが知られています。

# 楕円曲線のパラメータ

　この挙動を定量的に解析するためには、先の運動方程式を解く必要がありますが、そのあたりは先人の知恵があるので、要点だけを拝借します。楕円軌道の各パラメータは以下のようになります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/32d45510-3efb-6673-753d-3f2dd074ba8d.png)
（https://www.eee.kagoshima-u.ac.jp/~watanabe-lab/misc/%E6%A5%95%E5%86%86%E3%81%A8%E5%B9%B3%E5%9D%87.pdf より引用）

　半通経 $l$ は以下のように計算されます。（なお、以下のパラメータ群は、__恒星の質量$M$を都合のいい場合に設定し、垂直方向にのみ初速度ある場合に成り立つ特殊値__であることに注意してください）

```math 
l = r^2v^2
```

　離心率$e$はどうなるでしょうか。これは計算が少しめんどくさいのですが、結論を述べると、

```math 
e=\sqrt{1+r^2v^4-2rv^2}
```

　になります。要は、$v=\pm\frac{1}{\sqrt{r}}$の時に最小値0を取り、$v=\pm\sqrt{\frac{2}{r}}$の時に1となり、後は増え続けるだけになります。

　離心率は、

- 0の時は真円に
- 1未満の時は楕円に
- 1の時に放物線に
- 1より大きい時は双曲線

　になるので、シミュレーションの結果と一致していることがわかります。

　余談になりますが、以下のケースは高校物理の範囲でも計算することができます。

- 真円...遠心力と重力が釣り合う時

```math 
\displaylines{
\frac{m}{r^2} = \frac{mv^2}{r} \\
v = \pm\frac{1}{\sqrt{r}}}
```

- 放物線...無限遠で運動エネルギー＋位置エネルギーが0になる時

```math 
\displaylines{
0 = \frac{1}{2}mv^2 - \frac{m}{r} \\
v = \pm\sqrt{\frac{2}{r}}}
```

　さて、話を戻して、楕円のパラメータをもう少し詳しく見てみます。

　軌道長半径$a$は、

```math 
\displaylines{
a = \frac{l}{1-e^2} \\
=\frac{r^2v^2}{1-(1+r^2v^4-2rv^2)} \\
=\frac{r^2v^2}{2rv^2-r^2v^4} \\
=\frac{r}{2-rv^2}}
```

 になります。

　$a=f(v)$をグラフにするとこんな感じです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/03c47a3b-8fd5-1d22-0ef2-585236fffcbd.png)

　`vy[0]=0.1`で半径100の真円になります。ちなみに、0.1以下の値だと、原点を「左側の」焦点とした楕円運動に、0.1以上だと、原点を「右側の」焦点とした楕円運動になります。軌道長半径は大きくなり続けていきますが、`vy[0]`が$\frac{\sqrt{2}}{10}$ でついに無限半径の楕円軌道、つまり放物線に至ります。

　そこからが少し奇妙です。**軌道長半径がマイナス**となってしまいました。マイナスの半径、というのはあり得るのでしょうか。そして実際の運動は双曲線となっているのですが、これらの物理的関係を見出すことはできるのでしょうか。

# 周期

　このあたりを深堀りするために、円軌道のもう1つの要素である**周期**を考えてみたいと思います。

　ケプラーの法則より、惑星の公転周期$P$は、

```math 
\frac{a^3}{P^2} = \frac{G(M+m)}{4\pi^2}
```
　となりますが、今回は都合の良い物理定数を採用しているので、

```math  
P = 2\pi{\sqrt{(\frac{r}{2-rv^2})}}^3
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/5e71a169-585f-b309-8ef8-3c7dacda6e74.png)


　こんな感じになります。

　`vy[0] = 0.1``size=62831`の時、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/b9e8d78a-6f99-5b5c-8a5a-4879998312a5.png)

　`vy[0] = 0.12``size=15109`の時、

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/288125/ef2a8270-6a91-1931-275e-e16a8d46241c.png)

　という感じに、周期を含めて正確にシミュレートできていることがわかります。

# 虚数

　ここでやっとメインテーマである虚数に移ります。$v$がある値を超えると、軌道長半径$a$がマイナスになりました。

　ということは、計算上は、その平方根に比例する周期は**虚数の周期**となります。観測上では双曲線運動となっていますが。

　虚数を用いれば、$v$がある値を超えると運動の種類が変わるという「場合分け」は必要なく、「惑星は$v$の値に関わらず円運動をし、その周期は軌道長半径の3乗の平方根に比例する」という**1つの記述で対応できる**ことになります。

　前回の漸化式の例では「基本的には指数関数で、ある場合に振動が現れるが、それは虚数を含めて考えたら指数関数の一種である」という整理ができましたが、今回は「基本的には円軌道で、ある場合に双曲線が現れるが、それは虚数を含めて考えたら円運動の一種である」という整理ができました。また、その「境目」には、直線が現れたり、放物線が現れたりします。

　もちろんこれは同じ構造を持っていて、虚数の指数が振動的な挙動を示す、またその逆も然りということは、オイラーの公式から自然に推論できます。

# 余談

　軌跡を計算する際に、角運動量 $\vec{L}=\vec{r} \times \vec{v}$ を計算すると、角運動量保存則が（楕円や放物線においても）成り立っていることがわかります。

```python
    rvec = np.array([x[i],y[i],0])
    vvec = np.array([vx[i],vy[i],0])
    momentum = np.cross(rvec,vvec)
    print(momentum)
```

　以下は、`vy[0]=0.1414` の時。

```python
[ 0.          0.         14.14]
[ 0.         -0.         14.14001414]
[ 0.         -0.         14.14002828]
[ 0.         -0.         14.14004242]
[ 0.         -0.         14.14005656]
[ 0.         -0.         14.1400707]
[ 0.         -0.         14.14008484]
[ 0.         -0.         14.14009898]
[ 0.         -0.         14.14011312]
[ 0.         -0.         14.14012726]...
```