to geomland.isec_circlecircle :a :b
;01 (a.x-X)^2 + (a.y-Y)^2 = a.r^2
;02 (b.x-X)^2 + (b.y-Y)^2 = b.r^2
;03=01-02 ((a.x-X)-(b.x-X))*((a.x-X)+(b.x-X)) + ((a.y-Y)-(b.y-Y))*((a.y-Y)+(b.y-Y)) = a.r^2-b.r^2
;04 (a.x-X-b.x+X)*(a.x-X+b.x-X) + (a.y-Y-b.y+Y)*(a.y-Y+b.y-Y) = a.r^2-b.r^2
;05 (a.x-b.x)*(a.x+b.x-2*X) + (a.y-b.y)*(a.y+b.y-2*Y) = a.r^2-b.r^2
;06 (a.x-b.x)*(a.x+b.x)-2*(a.x-b.x)*X + (a.y-b.y)*(a.y+b.y)-2*(a.y-b.y)*Y = a.r^2-b.r^2
;06 -2*(a.x-b.x)*X -2*(a.y-b.y)*Y + a.x^2-b.x^2 + a.y^2-b.y^2 = a.r^2-b.r^2
;08 A=-2*(a.x-b.x) B=-2*(a.y-b.y) C=+a.r^2 -b.r^2 -a.x^2 +b.x^2 -a.y^2 +b.y^2
make local "A -2 * (:a.center.x-:b.center.x)
make local "B -2 * (:a.center.y-:b.center.y)
make local "C :a.radius*:a.radius - :b.radius*:b.radius - :a.center.x*:a.center.x + :b.center.x*:b.center.x - :a.center.y*:a.center.y + :b.center.y*:b.center.y
if :A<>0
[
;09 A*X+B*Y=C
;10 X=(C-B*Y)/A= C/A-(B/A)*Y (if A not 0)
;11(10,01) (a.x-C/A+(B/A)*Y)^2 + (a.y-Y)^2 = a.r^2
;12 (a.x-C/A)^2 +2(a.x-C/A)*(B/A)*Y + ((B/A)*Y)^2 + (a.y^2-2*a.y*Y+Y^2 = a.r^2
;13 ((B/A)^2+1)*Y^2+(2(a.x-C/A)*(B/A)-2*a.y)*Y+((a.x-C/A)^2+a.y^2-a.r^2)=0
;14 P=((B/A)^2+1) Q=((a.x-C/A)*(B/A)-a.y) R=((a.x-C/A)^2+a.y^2-a.r^2)
; (multiply back with A^2
;14a P=B^2+A^2 Q=(a.x*A-C)*B-a.y*A^2 R=(a.x*A-C)^2+(a.y^2-a.r^2)*A^2
make local "P :B*:B+:A*:A
make local "Q (:A*:a.center.x-:C) * :B-:a.center.y*:A*:A
make local "R (:A*:a.center.x-:C) * (:A*:a.center.x-:C) +:A*:A * (:a.center.y-:a.radius) * (:a.center.y+:a.radius)
;15 P*Y^2 + 2*Q*Y +R = 0
;16 D=sqrt(Q^2-P*R) (if Q^2>=P*R)
make local "D :Q*:Q-:P*:R
; D<0 => no intersection
if :D<0 [output set]
; D=0 => 1 point
if :D=0
[
make local "Y1 (-:Q/:P)
make local "X1 (:C-:B*:Y1)/:A
output (set point :X1 :Y1)
]
;16a D=sqrt(D)
make "D sqrt :D
;17 Y1 = (-Q+D)/P X1 = (C-B*Y1)/A
;18 Y2 = (-Q-D)/P X2 = (C-B*Y2)/A
make local "Y1 (-:Q+:D)/:P
make local "X1 (:C-:B*:Y1)/:A
make local "Y2 (-:Q-:D)/:P
make local "X2 (:C-:B*:Y2)/:A
output (set point :X1 :Y1 point :X2 :Y2)
]
[
;09 A*X+B*Y=C (A=0)
;10 B*Y=C
if :B=0
[
if :C=0
[ output :a ]
[ output set ]
]
make local "Y :C/:B
;11 (a.x-X)^2 = a.r^2-(a.y-Y)^2
;12 (a.x-X)^2 = (a.r-a.y+Y)*(a.r+a.y-Y)
;13 (a.x-X) = +- sqrt[(a.r-a.y+Y)*(a.r+a.y-Y)]
;14 X = a.x -+ sqrt[(a.r-a.y+Y)*(a.r+a.y-Y)]
make local "X1 :a.center.x-(sqrt (:a.radius-:a.center.y+:Y)*(:a.radius+:a.center.y-:Y))
make local "X2 2*:a.center.x-:X1
output (set point :X1 :Y point :X2 :Y)
]
end
|