楕円の周長と一方の半径から残された半径を求める python

# 概要
周長 L
長半径 a
短半径 b
の楕円で、Lとaまたはbが分かっているときに、残された半径を求める。
![ellipes.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2778069/083bb39f-4291-66fb-0f7d-e44035f0030a.png)

# 環境
python 3.10.6
sympy==1.11.1
numpy==1.21.5

科学技術計算ライブラリであるscipyをインストールしておく
```bash
pip install scipy
```

# コード


```python
import scipy as sp
import numpy as np

def get_ellipse_another_radius(L,r):
  """
  楕円の周長と半径rから残されたほうの半径を求める
  L:楕円の周長
  r:長半径または短半径
  """
  k = 0.9 # 離心率初期値
  dk = 0.001  # デルタk
  # rが長半径か短半径か判定する
  if 2*np.pi*r > L:
    # rが長半径の場合
    a = r
    # 楕円の周の長さの公式
    def F(k):return 4*a*sp.special.ellipe(k**2)-L 
    # F(k)の微分
    def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
    # ニュートン・ラプソン法で離心率を求める
    for i in range(10):
      k = k - F(k)/F_diff(k)
    # 離心率の式から短半径を求める
    b = a*np.lib.scimath.sqrt(1-k**2)
    return b
  else:
    # rが短半径の場合
    b = r
    # 楕円の周の長さの公式
    def F(k):return 4*b/np.lib.scimath.sqrt(1-k**2)*sp.special.ellipe(k**2)-L 
    # F(k)の微分
    def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
    # ニュートン・ラプソン法で離心率を求める
    for i in range(10):
      k = k - F(k)/F_diff(k)
    # 離心率の式から長半径を求める
    a = b/np.lib.scimath.sqrt(1-k**2)
    return a
```

# 解説
楕円の基本公式
```math
第二種完全楕円積分\\
E(k)=\int_{0}^{\frac{\pi}{2}}\sqrt{1-k^2\sin^2t}dt\\
楕円の周長の公式\\
L=4aE(k)\\
離心率\\
k=\sqrt{1-\frac{b^2}{a^2}}\hspace{20px},\hspace{20px}(b<a)
```
一方の半径rが長半径aか短半径bかで計算が変わるため、まずはrが長半径か短半径かを判定する。
それには下の図を見ると分かりやすい。
![ellipes2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2778069/789117d5-c542-7b75-5ca0-b94052b8b7e0.png)

2πr > Lのときはrが長半径である。
まずはrが長半径の場合について求める。
## rが長半径の場合
```python
# rが長半径の場合
a = r
```
離心率は変形すると
```math
b=a\sqrt{1-k^2}
```
```python
# 離心率の式から短半径を求める
b = a*np.lib.scimath.sqrt(1-k**2)
```
となり、aの値は分かっているので離心率kが分かればbが求まる。

離心率kを求めるにはニュートン・ラフソン法を用いる
```math
ニュートン・ラフソン法\\
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
```
詳しくは解説しないが、f(x)=0のときのxの近似値を反復計算によって得ることができる手法である。

楕円の周長の公式を次のようにする
```math
f(k)=4aE(k)-L=0
```
```python
# 楕円の周の長さの公式
def F(k):return 4*a*sp.special.ellipe(k**2)-L
```
第二種完全楕円積分はscipy.special.ellipeで計算できる。※m=k^2 となっているので、k^2 を代入するところだけは注意
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html

離心率kをニュートン・ラフソン法に当てはめると下記になる。
```math
k_{n+1} = k_n-\frac{f(k_n)}{f'(k_n)}\\
k_0 = 0.9とする\\

```
計算を繰り返していくと離心率kの近似値が求まる。
今回は１０回繰り返してみる
```python
k = 0.9 # 離心率初期値
dk = 0.001  # デルタk
```
```python
# F(k)の微分
def F_diff(k):return (F(k+dk)-F(k-dk))/(2*dk)
# ニュートン・ラフロソン法で離心率を求める
for i in range(10):
  k = k - F(k)/F_diff(k)
```

kの近似値が得られたのでbを求めることができる
```python
# 離心率の式から短半径を求める
b = a*np.lib.scimath.sqrt(1-k**2)
```

## rが短半径の場合
```python
# rが長半径の場合
b = r
```
離心率を変形して
```math
a=\frac{b}{\sqrt{1-k^2}}
```
aを、楕円の周長の公式に代入して長半径の場合と同じようにする
```math
f(k)=\frac{4b}{\sqrt{1-k^2}}E(k)-L=0
```
あとは同じくニュートン・ラフソン法を用いて離心率kを計算し、長半径aが求まる。

以下省略

# まとめ
楕円の周長と一方の半径rから残された半径を求める方法について記述しました。
- rが長半径なのか短半径なのか判別する。
- ニュートン・ラフソン法に楕円の周長の公式を当てはめて、離心率kの近似値を求める。
- 離心率の公式からもう一方の半径を求める。


# 参考文献
- [”楕円の円周と一方の径からもう一方の径を求める。”をカシオの高精度計算サイトkeisan.casio.jpに自作式としてUP.](https://sci.tea-nifty.com/blog/2020/01/post-9c35a7.html)
- [scipy.special.ellipe](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipe.html)
- [微分・積分の数学　ニュートン・ラフソン法　［前編］](https://gihyo.jp/dev/serial/01/java-calculation/0067#:~:text=%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E3%83%BB%E3%83%A9%E3%83%95%E3%82%BD%E3%83%B3%E6%B3%95%5B1%5D,%E5%80%A4%E3%82%92%E5%BE%97%E3%82%89%E3%82%8C%E3%81%BE%E3%81%99%E3%80%82)
- [ウィキペディア 周長](https://ja.wikipedia.org/wiki/%E5%91%A8%E9%95%B7#:~:text=%E6%A5%95%E5%86%86%E3%81%AE%E5%91%A8%E9%95%B7,-%E3%80%8C%E6%A5%95%E5%86%86%23%E6%A5%95%E5%86%86%E3%81%AE&text=%E6%A5%95%E5%86%86%20x2%2Fa2,%E6%96%B9%E3%81%AE%E4%B8%80%E4%BE%8B%E3%82%92%E7%A4%BA%E3%81%99%E3%80%82)