2つの円の交点を求める

## 2つの円の位置関係
2つの円の位置関係は、

-  離れている
-  外接する
-  2点で交わる
-  内接する
-  一方の円の内部にある

の５パターン。

<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/f030e5c6-36a3-dc9e-2b40-a1f360436156.png"> <img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/74a4ceee-b73e-2d4b-3e9b-46cc9c8a3767.png"> <img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/b3779135-29ec-d4bd-a3be-85d7c5f37d60.png"> <img width="223" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/68d6c636-5c80-385d-da31-bf9d504c88c3.png"> <img width="223" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/62942548-927e-ea01-dfd2-afc94b89404f.png"> 

位置関係は、2つの円の半径と、2つの円の中心間の距離を比べることで判定することができる。
2つの円の半径を $r_1,r_2$、2つの円の中心間の距離を $d$ とすると、

-  離れている $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots \ \ d > r_1+r_2$
-  外接する $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots \ \ d = r_1+r_2$
-  2点で交わる $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots \ \ |r_1-r_2| < d < r_1+r_2$
-  内接する $\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots \ \ d = |r_1+r_2|$
-  一方の円の内部にある $\ \cdots \ \ d < |r_1+r_2|$

となる。

## サンプルプログラム
#### 2つの円を選択して、その位置関係を調べる。

```scheme
;; 2つの円を選択して、その位置関係を調べる。
(defun c:PositionOfTwoCircles (/ circle1 circle2 n)
    (setq circle1 (editor:EntselCircle "1つ目の円を選択")
          circle2 (editor:EntselCircle "2つ目の円を選択")
    )
    
    ;; 共通接線の数を調べる
    (setq n (circle:CountCommonTangents circle1 circle2))
    (alert
        (cond
            ((= n 0) "一方の円の内部にある")
            ((= n 1) "内接する")
            ((= n 2) "2点で交わる")
            ((= n 3) "外接する")
            ((= n 4) "離れている")
        )
    )
)
```

### Subfunctions

```scheme
;; == Subfunctions ==
;; 図形のDXF定義データの値を取得
(defun entity:GetDxf (g e)
    (cond
        ((= (type e) 'ENAME) (cdr (assoc g (entget e))))
        ((listp e) (cdr (assoc g e))) 
    )
)
;; 図形のタイプを取得
(defun entity:GetType (e) (entity:GetDxf 0 e))
;; 円図形をentsel 選択
(defun editor:EntselCircle (msg / circle)
    ;; 円を選択
    (while (null circle)
        (setvar 'ERRNO 0)
        (setq circle (car (entsel (strcat "\n" msg " :\n"))))
        (cond
            ;; 空振り
            ((= 7 (getvar 'ERRNO))
             (princ "\n... 空振り! 再選択 ...\n")
             (setq circle nil)
            )
            ;; 空Enter
            ((= 52 (getvar 'ERRNO))
             (princ "\n... 空Enter! 再選択 ...\n")
             (setq circle nil)
            )
            ;; 図形が選択された
            (circle
             (if (/= "CIRCLE" (entity:GetType circle))
                 (progn
                     (princ "\n... 円ではありません! 再選択 ...\n")
                     (setq circle nil)
                 )
                 (princ "\n... 円が選択されました。\n")
             )
            )
        )
    )
    circle
)
;; 円の中心を取得
(defun circle:GetCenterPoint (e) (entity:GetDxf 10 e))
;; 円の半径を取得
(defun circle:GetRadius (e) (entity:GetDxf 40 e))
;; 共通接線の本数を返す
(defun circle:CountCommonTangents (circle1 circle2 / r1 r2 o1 o2 d)
    (setq r1 (circle:GetRadius circle1)
          o1 (circle:GetCenterPoint circle1)
          r2 (circle:GetRadius circle2)
          o2 (circle:GetCenterPoint circle2)
          d  (distance o1 o2)
    )
    (cond
        ((< d (abs (- r1 r2))) 0)                       ;_一方の円の内部にある
        ((equal d (abs (- r1 r2)) 1e-10) 1)             ;_内接する
        ((and (< (abs (- r1 r2)) d) (< d (+ r1 r2))) 2) ;_2点で交わる
        ((equal d (+ r1 r2) 1e-10) 3)                   ;_外接する
        ((> d (+ r1 r2)) 4)                             ;_離れている
        (T nil)
    )
)
```

## 2つの円の交点を通る直線の方程式を求める
2つの円の位置関係が、2点で交わる場合の、2つの交点を通る直線の方程式を求める。

中心の座標が $(x_1,y_1)$ 半径 $r_1$ の円 $(x-x_1)^2+(y-y_1)^2=r_1^2 \ \cdots \ ①$ と、
中心の座標が $(x_2,y_2)$ 半径 $r_2$ の円 $(x-x_2)^2+(y-y_2)^2=r_2^2 \ \cdots \ ②$
の位置関係が、2点で交わる場合の、2つの交点を通る直線の方程式を求める。
<img width="600" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/384860/1477a101-0347-d6e6-6350-f8d82fe03042.png">

式①②の連立方程式として求める。
円の方程式 ①から②を引く、
$$(x-x_1)^2-(x-x_2)^2+(y-y_1)^2-(y-y_2)^2=r_1^2-r_2^2$$
展開して整理すると、
$$(2x_2-2x_1)x+(2y_2-2y_1)y+(x_1^2-x_2^2+y_1^2-y_2^2+r_2^2-r_1^2)=0$$
さらに展開して整理すると、
$$2(x_2-x_1)x+2(y_2-y_1)y+\\{(x_1+x_2)(x_1-x_2)+(y_1+y_2)(y_1-y_2)+(r_2+r_1)(r_2-r_1)\\}=0$$
となり、この式は

- $a=2(x_2-x_1)$
- $b=2(y_2-y_1)$
- $c=\\{(x_1+x_2)(x_1-x_2)+(y_1+y_2)(y_1-y_2)+(r_2+r_1)(r_2-r_1)\\}$

とした[直線の方程式](https://qiita.com/tydesign/items/4d3909a826363d8fc29f)の一般形 $ax+by+c=0$ の形になっている。
2つの円の交点を通る直線は1つしかないので、この式がその直線の方程式となる。

#### 【例】中心の座標が (3,2) 半径 5 の円と中心の座標が (9,4) 半径 3 の円の2つの交点を通る直線の方程式を求める。
この2つの円の方程式はそれぞれ、
$$(x-3)^2+(y-2)^2=25$$
$$(x-9)^2+(y-4)^2=9$$
上記公式、
$$2(x_2-x_1)x+2(y_2-y_1)y+\\{(x_1+x_2)(x_1-x_2)+(y_1+y_2)(y_1-y_2)+(r_2+r_1)(r_2-r_1)\\}=0$$
に値を代入すると、
$$2(9-3)x+2(4-2)y+\\{(3+9)(3-9)+(2+4)(2-4)+(3+5)(3-5)\\}=0$$
$$12x+4y-100=0$$
よって、
$$3x+y-25=0$$
これが円の交点を通る直線の方程式です。

## 2つの円の交点を求める
上記の2つの円の交点を通る直線の方程式を求める公式を使って、その直線の方程式を求める。次にその直線とどちらかの円との交点を、[円と直線の交点を求める](https://qiita.com/tydesign/items/36b42465b3f5086bd0c5)公式を使って求める。

#### 【例】中心の座標が (3,2) 半径 5 の円と中心の座標が (9,4) 半径 3 の円の2つの交点を求める。
交点を通る直線の方程式は、上記により、
$$3x+y-25=0$$
この直線と、円 $(x-3)^2+(y-2)^2=25$ との交点を求める。

中心の座標が $(x_0,y_0)$ 半径 $r$ の円 $(x-x_0)^2+(y-y_0)^2=r^2$と、直線 $ax+by+c=0$ の交点を $A,B$ とすると、その座標は、
$$A=\left(\frac{aD-b\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+x_0,\frac{bD+a\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+y_0\right)$$
$$B=\left(\frac{aD+b\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+x_0,\frac{bD-a\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+y_0\right)$$
$$ただし、D=\mid ax_0+by_0+c\mid$$

この公式を使って求める。
まず $D$ を求める
$$D=\mid ax_0+by_0+c\mid=\mid 3\cdot3+2\cdot1-25\mid=\mid-14\mid=14$$
点 $A$ の $x$ 座標を求める
$$x_A=\frac{aD-b\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+x_0=\frac{3\cdot14-\sqrt{(3^2+1^2)5^2-14^2}}{3^2+1^2}+3$$
$$=\frac{42-\sqrt{54}}{10}+3=\frac{72-3\sqrt{6}}{10}\fallingdotseq 6.465153$$
点 $A$ の $y$ 座標を求める
$$y_A=\frac{bD+a\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+y_0=\frac{1\cdot14+3\sqrt{(3^2+1^2)5^2-14^2}}{3^2+1^2}+2$$
$$=\frac{14+3\sqrt{54}}{10}+2=\frac{34+9\sqrt{6}}{10}\fallingdotseq 5.604541$$
点 $B$ の $x$ 座標を求める
$$x_B=\frac{aD+b\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+x_0=\frac{3\cdot14+\sqrt{(3^2+1^2)5^2-14^2}}{3^2+1^2}+3$$
$$=\frac{42+\sqrt{54}}{10}+3=\frac{72+3\sqrt{6}}{10}\fallingdotseq 7.934847$$
点 $B$ の $y$ 座標を求める
$$y_B=\frac{bD-a\sqrt{(a^2+b^2)r^2-D^2}}{a^2+b^2}+y_0=\frac{1\cdot14-3\sqrt{(3^2+1^2)5^2-14^2}}{3^2+1^2}+2$$
$$=\frac{14-3\sqrt{54}}{10}+2=\frac{34-9\sqrt{6}}{10}\fallingdotseq 1.195459$$

したがって、
$$A=\left(\frac{72-3\sqrt{6}}{10}, \frac{34+9\sqrt{6}}{10} \right)\fallingdotseq (6.465153,5.604541)$$
$$B=\left(\frac{72+3\sqrt{6}}{10}, \frac{34-9\sqrt{6}}{10} \right)\fallingdotseq (7.934847,1.195459)$$
となる。