勝手に添削:続・円の交点を求めるプログラム

三角関数は必要ないので消してみた編。円の中心を結ぶ線と、円の交点を結ぶ線は直行するので、直角三角形にピタゴラスの定理を使ってまず円1の中心から「円の中心を結ぶ線と円の交点を結ぶ線の交点」までの距離を求める。それがk。あとは交点kdirを求めてもう一度ピタゴラスの定理を使ってkdirから直角方向にどれだけずれればいいかを求めている。途中のコメントに書いた数式を参考のこと。a, bってのが円の中心を結んだ線分を交点で分けた2つの線分の長さ。あと円が交差しないときのチェックも必要だと思うけど省いた。mag <= r1 + r2 の時に交点がある…とは限らない(完全に中に含まれるケースがあるから)

import math
size(300, 300)
background(0.3, 0.3, 0.3)
nofill()
stroke(1, 0.1, 0)

def writeEn(posX, posY, R):
    x = -R + posX
    y = -R + posY
    oval(x, y, R * 2, R * 2)

def scale(k ,vec):
    return (k * vec[0], k * vec[1])

def add(v1, v2):
    return (v1[0] + v2[0], v1[1] + v2[1])
    
def circleCross(xc1, yc1, r1, xc2, yc2, r2):
    from math import sqrt
    cross = []
    dir = (float(xc2 - xc1), float(yc2 - yc1))
    sqmag = dir[0] ** 2 + dir[1] ** 2
    r1sq = r1 ** 2
    r2sq = r2 ** 2
    # r1sq == a ** 2 + d ** 2
    # r2sq == b ** 2 + d ** 2
    # sqmag == (a + b) ** 2
    # r1sq - r2sq == a ** 2 - b ** 2 == (a + b) * (a - b)
    ratio = (r1sq - r2sq + sqmag) / sqmag # 2a / (a + b)
    mag = sqrt(sqmag)
    k = ratio * mag / 2
    d = sqrt(r1sq - k ** 2)
    ndir = scale(1 / mag, dir)

    kdir = add(scale(k, ndir), (xc1, yc1))
    writeEn(kdir[0], kdir[1], 2)

    odir = scale(d, (ndir[1], -ndir[0]))
    
    return [
        kdir[0] + odir[0], kdir[1] + odir[1], 
        kdir[0] - odir[0], kdir[1] - odir[1], 
    ]

c1 = [100,200,100]
c2 = [200,100,60]
writeEn(c1[0],c1[1],c1[2])
writeEn(c2[0],c2[1],c2[2])

p1= circleCross(c1[0],c1[1],c1[2],
                c2[0],c2[1],c2[2])
stroke(0.5, 0.1, 0)
writeEn(p1[0], p1[1], 10)
writeEn(p1[2], p1[3], 10)