to geomland.isec_linecircle :l :k
;01 line (x1,y1)-(x2,y2) circle (xr,yr,r)
;02 (x-x1)/(x2-x1)=(y-y1)/(y2-y1)
;03 a=y2-y1 b=x1-x2 c=x2*y1-x1*y2
;04 a*x+b*y+c=0
make local "a :l.final.y-:l.initial.y
make local "b :l.initial.x-:l.final.x
make local "c :l.final.x*:l.initial.y-:l.initial.x*:l.final.y
;05 if a=0 => b*y+c=0
if :a=0
[
if :b=0 [ output set ]
make local "Y (-:c/:b)
make local "D :k.radius*:k.radius - (:Y-:k.center.y) * (:Y-:k.center.y)
if :D<0 [output set]
if :D=0 [output point :k.center.x :Y]
make "D sqrt :D
make local "X1 :k.center.x+:D
make local "X2 :k.center.x-:D
output (set point :X1 :Y point :X2 :Y)
]
;06 if a<>0
;07 a*x=-c-b*y
;08 (x-xr)^2+(y-yr)^2-r^2=0
;09 (a*x-a*xr)^2+(a*y-a*yr)^2-a^2*r^2=0
;10 (-c-b*y-a*xr)^2+(a*y-a*yr)^2-a^2*r^2=0
;11 b*b*y*y+2*b*y*(c+a*xr)+(c+a*xr)^2+a*a*y*y-2*a*a*yr*y+a*a*yr*yr-a*a*r*r = 0
;11 (a*a+b*b)*y^2+2*(b*c+a*b*xr-a*a*yr)*y+(c+a*xr)^2+a^2*(yr-r)*(yr+r) = 0
;12 A=a^2+b^2<>0 B=b*(c+a*xr)-a^2*yr C=(c+a*xr)^2+a*a*(yr*yr-r*r)
;13 A*y*y+2*B*y+C=0
make local "A :a*:a+:b*:b
make local "B :b * (:c+:a*:k.center.x) - :a*:a*:k.center.y
make local "C (:c+:a*:k.center.x) * (:c+:a*:k.center.x) + :a*:a * (:k.center.y-:k.radius) * (:k.center.y+:k.radius)
make local "D :B*:B-:A*:C
if :D<0 [ output set ]
if :D=0
[
make local "Y (-:B/:A)
make local "X (-:b*:Y-:c) / :a
output point :X :Y
]
make "D sqrt :D
make local "Y1 (-:B+:D) / :A
make local "Y2 (-:B-:D) / :A
make local "X1 (-:b*:Y1-:c) / :a
make local "X2 (-:b*:Y2-:c) / :a
output (set point :X1 :Y1 point :X2 :Y2)
end
|