楕円の内接円，曲率円

# 楕円の曲率円

曲率円とは，楕円の長径端においてその楕円に内接する円のことである。

## 1.　楕円の内接円

説明を始めるにあたり，まずは楕円に 2 点で内接する円について調べる。

楕円の長径と短径を a, b，中心座標を (0, 0) とする。
楕円に 2 点で接する円の半径を r，中心座標を (d, 0) とする。

```julia
using SymPy

@syms a, b, r, d, x, y

eq1 = x^2/a^2  + y^2/b^2 - 1
eq2 = (x - d)^2 + y^2 - r^2;
```

eq1, eq2 の連立方程式を解くと,4対の交点座標 (x, y) が得られる。

```julia
res = solve([eq1, eq2], (x, y));
```

そのうちの x > 0 である 2 対は以下のようになる。分子のルートの中身はいわゆる判別式で，この値が 0 のときには x は重解で，楕円と円が第一象限の一点で交わることになる。

```julia
res[1][1] |> simplify |> display
res[1][2] |> simplify |> display
```

$\frac{a \left(a d + \sqrt{- a^{2} b^{2} + a^{2} r^{2} + b^{4} + b^{2} d^{2} - b^{2} r^{2}}\right)}{a^{2} - b^{2}}$

$- \sqrt{- \frac{b^{2} \left(- a^{4} + a^{2} b^{2} + a^{2} d^{2} + a^{2} r^{2} + 2 a d \sqrt{- a^{2} b^{2} + a^{2} r^{2} + b^{4} + b^{2} d^{2} - b^{2} r^{2}} + b^{2} d^{2} - b^{2} r^{2}\right)}{a^{4} - 2 a^{2} b^{2} + b^{4}}}$

```julia
res[4][1] |> simplify |> display
res[4][2] |> simplify |> display
```

$\frac{a \left(a d - \sqrt{- a^{2} b^{2} + a^{2} r^{2} + b^{4} + b^{2} d^{2} - b^{2} r^{2}}\right)}{a^{2} - b^{2}}$

$\sqrt{\frac{b^{2} \left(a^{4} - a^{2} b^{2} - a^{2} d^{2} - a^{2} r^{2} + 2 a d \sqrt{- a^{2} b^{2} + a^{2} r^{2} + b^{4} + b^{2} d^{2} - b^{2} r^{2}} - b^{2} d^{2} + b^{2} r^{2}\right)}{a^{4} - 2 a^{2} b^{2} + b^{4}}}$

判別式が 0 になる条件を求めるには，上述のようにするより以下のほうが直接的であろう。

eq2 を y についてとき，eq1 に代入すると x についての二次式になる。x の 2 次，1 次，0 次の係数c2, c2, c0 を取り出し，判別式 D = c^2 - 4c1*c0 を求め， D = 0 を d について解く。

```julia
eq = eq1(y => sqrt(r^2 - (x - d)^2))*a^2*b^2 |> expand
```

$- a^{2} b^{2} - a^{2} d^{2} + 2 a^{2} d x + a^{2} r^{2} - a^{2} x^{2} + b^{2} x^{2}$

```julia
c2 = collect(eq, x).coeff(x, 2)
```

$- a^{2} + b^{2}$

```julia
c1 = collect(eq, x).coeff(x, 1)
```

$2 a^{2} d$

```julia
c0 =collect(eq, x).coeff(x, 0) |> simplify
```

$a^{2} \left(- b^{2} - d^{2} + r^{2}\right)$

```julia
D = c1^2 - 4c2*c0
```

$4 a^{4} d^{2} - a^{2} \left(- 4 a^{2} + 4 b^{2}\right) \left(- b^{2} - d^{2} + r^{2}\right)$

```julia
res2 = solve(D, d)
```

$\left[ \begin{array}{r}- \frac{\sqrt{\left(a - b\right) \left(a + b\right) \left(b - r\right) \left(b + r\right)}}{b}\\,,\frac{\sqrt{\left(a - b\right) \left(a + b\right) \left(b - r\right) \left(b + r\right)}}{b}\end{array} \right]$

d = sqrt((a^2 - b^2)*(b^2 - r^2)) が得られる。

```julia
@syms a, b, r, d, x, y

eq1 = x^2/a^2  + y^2/b^2 - 1
eq2 = (x - d)^2 + y^2 - r^2;
res = solve([eq1, eq2], (x, y))
#(a, b, r) = (3, 1, 0.8)
d = sqrt((a^2 - b^2)*(b^2 - r^2))/b
d(a => 3, b => 1, r => 0.9) |> println
(res[1][1](a => 3, b => 1, r => 0.9))(d=> 1.69705627484771) |> N
```

    1.23288280059380

$\frac{\frac{9 d^{2}}{4} + 0.924662100445346 d \sqrt{0.657894736842105 d^{2} - 1}}{2 d}$

    a = 3;  b = 1;  r = 0.9
    d = 1.23288;  (x, y) = (1.38699, 0.886707)

![output_16_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/bce8cf2c-cba1-88d9-ac88-c7bdd5837e03.png)


a, b, r から d, x, y を求める関数は以下の通り。

```julia
function f(a, b, r)
    d = sqrt((a^2 - b^2)*(b^2 - r^2))/b
    x = a^2*d/(a^2 - b^2)
    y = sqrt(r^2 - (x - d)^2)
    return (d, x, y)
end;

f(3, 1, 0.9)
```

    (1.232882800593795, 1.3869931506680195, 0.8867073925484099)

判別式を r で解くと，以下のようになる。

```julia
res3 = solve(D, r)
```

$\left[ \begin{array}{r}b \sqrt{\frac{a^{2} - b^{2} - d^{2}}{\left(a - b\right) \left(a + b\right)}}\\,,- b \sqrt{- \frac{- a^{2} + b^{2} + d^{2}}{\left(a - b\right) \left(a + b\right)}}\end{array} \right]$

a, b, d を与えて r, x, y を求める関数は以下の通り。

```julia
function g(a, b, d)
    r = b*sqrt((a^2 - b^2 - d^2)/((a^2 - b^2)))
    x = a^2*d/(a^2 - b^2)
    y = sqrt(r^2 - (x - d)^2)
    return (r, x, y)
end;
g(3, 1, 1.23288280059380)
```

    (0.8999999999999991, 1.3869931506680249, 0.8867073925484089)

## 2.　曲率円

円の半径を小さくしていくと長径の端で 1 点で接するようになる。このときの円が，曲率円である。
曲率円の半径を r とすると中心は (a - r, 0) である。

```julia
include("julia-source.txt");

using SymPy

@syms a, b, r, d, x, y

res3 = solve(sqrt((a^2 - b^2)*(b^2 - r^2)) - (a - r), r)
```

$\left[ \begin{array}{r}\frac{a - \sqrt{\left(b - 1\right) \left(b + 1\right)} \left(a - b\right) \left(a + b\right)}{a^{2} - b^{2} + 1},\\
\frac{a + \sqrt{\left(b - 1\right) \left(b + 1\right)} \left(a - b\right) \left(a + b\right)}{a^{2} - b^{2} + 1}\end{array} \right]$

d = a - r となるときの半径を求める。

```julia
solve(sqrt((a^2 - b^2)*(b^2 - r^2))/b - (a - r), r)
```

$\left[ \begin{array}{r}\frac{b^{2}}{a}\end{array} \right]$

a, b を与えて曲率円の半径を求める関数は以下の通り，

```julia
h(a, b) = b^2/a
h(3, 1)
```

    0.3333333333333333



    a = 3;  b = 1;  r = 0.9
    d = 1.23288;  (x, y) = (1.38699, 0.886707);  曲率円の半径(R) = 0.333333

![output_29_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/be9549d1-431d-3011-0db0-0ed77e4c3c6e.png)