Line circle intersection: Difference between revisions
SqrtNegInf (talk | contribs) m (→{{header|Raku}}: Fix code: Perl 6 --> Raku) |
(→{{header|REXX}}: added the REXX computer programming language for this task.) |
||
Line 514:
For data set: [(5 2 −2.26) (0.77 2 4) (1 4 0) 4]
Solution(s) is/are: [(4.2615520237084015 2 -1.1671668246843006) (1.13386504516801 2 3.461514141193441)]
</pre>
=={{header|REXX}}==
The formulae used for this REXX version were taken from the MathWorld
webpage: [https://mathworld.wolfram.com/Circle-LineIntersection.html circle line intersection].
<lang rexx>/*REXX program calculates where (or if) a line intersects (or tengents) a cirle. */
/*───────────────────────────────────── line= x1,y1 x2,y2; circle is at 0,0, radius=r*/
parse arg x1 y1 x2 y2 cx cy r . /*obtain optional arguments from the CL*/
if x1=='' | x1=="," then x1= 0 /*Not specified? Then use the default.*/
if y1=='' | y1=="," then y1= -3 /* " " " " " " */
if x2=='' | x2=="," then x2= 0 /* " " " " " " */
if y2=='' | y2=="," then y2= 6 /* " " " " " " */
if cx=='' | cx=="," then cx= 0 /* " " " " " " */
if cy=='' | cy=="," then cy= 0 /* " " " " " " */
if r =='' | r =="," then r = 4 /* " " " " " " */
x_1= x1; x1= x1 + cx
x_2= x2; x2= x2 + cx
y_1= y1; y1= y1 + cy
y_2= y2; y2= y2 + cy
dx= x2 - x1; dy= y2 - y1
dr= sqrt(dx**2 + dy**2)
D= x1 * y2 - x2 * y1
ix1= ( D * dy + sgn(dy) * dx * sqrt(r**2 * dr**2 - D**2) ) / dr**2
ix2= ( D * dy - sgn(dy) * dx * sqrt(r**2 * dr**2 - D**2) ) / dr**2
iy1= (-D * dx + abs(dy) * sqrt(r**2 * dr**2 - D**2) ) / dr**2
iy2= (-D * dx - abs(dy) * sqrt(r**2 * dr**2 - D**2) ) / dr**2
incidence= (r**2 * dr**2 - D**2) / 1
say 'incidence=' incidence
@potla= 'points on the line are: '
if incidence<0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1
say "The line doesn't intersect the circle with radius: " r
end
if incidence=0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1
say "The line is tangent to circle with radius: " r
end
if incidence>0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1
say "The line is secant to circle with radius: " r
end
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
sgn: procedure; if arg(1)<0 then return -1; return 1
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
do j=0 while h>9; m.j= h; h= h%2 + 1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g= (g+x/g) *.5; end /*k*/
numeric digits d; return g/1</lang>
{{out|output|text= when using the default inputs:}}
<pre>
incidence= 1296
points on the line are: (0,-3) and (0,6) are: 0,4
The line is secant to circle with radius: 4
</pre>
|
Revision as of 14:49, 19 March 2020
In plane geometry, a line (or segment) may intersect a circle at 0, 1 or 2 points.
- Task
Implement a method (function, procedure etc.) in your language which takes as parameters:
- the starting point of a line;
- the point where the line ends;
- the center point of a circle;
- the circle's radius; and
- whether the line is a segment or extends to infinity beyond the above points.
The method should return the intersection points (if any) of the circle and the line.
Illustrate your method with some examples (or use the Go examples below).
- References
- See Math Stack Exchange for development of the formulae needed.
- See Wolfram for the formulae needed.
Go
<lang go>package main
import (
"fmt" "math"
)
const eps = 1e-14 // say
type point struct{ x, y float64 }
func (p point) String() string {
// hack to get rid of negative zero // compiler treats 0 and -0 as being same if p.x == 0 { p.x = 0 } if p.y == 0 { p.y = 0 } return fmt.Sprintf("(%g, %g)", p.x, p.y)
}
func sq(x float64) float64 { return x * x }
// Returns the intersection points (if any) of a circle, center 'cp' with radius 'r', // and either an infinite line containing the points 'p1' and 'p2' // or a segment drawn between those points. func intersects(p1, p2, cp point, r float64, segment bool) []point {
var res []point x0, y0 := cp.x, cp.y x1, y1 := p1.x, p1.y x2, y2 := p2.x, p2.y A := y2 - y1 B := x1 - x2 C := x2*y1 - x1*y2 a := sq(A) + sq(B) var b, c float64 var bnz = true if math.Abs(B) >= eps { // if B isn't zero or close to it b = 2 * (A*C + A*B*y0 - sq(B)*x0) c = sq(C) + 2*B*C*y0 - sq(B)*(sq(r)-sq(x0)-sq(y0)) } else { b = 2 * (B*C + A*B*x0 - sq(A)*y0) c = sq(C) + 2*A*C*x0 - sq(A)*(sq(r)-sq(x0)-sq(y0)) bnz = false } d := sq(b) - 4*a*c // discriminant if d < 0 { // line & circle don't intersect return res }
// checks whether a point is within a segment within := func(x, y float64) bool { d1 := math.Sqrt(sq(x2-x1) + sq(y2-y1)) // distance between end-points d2 := math.Sqrt(sq(x-x1) + sq(y-y1)) // distance from point to one end d3 := math.Sqrt(sq(x2-x) + sq(y2-y)) // distance from point to other end delta := d1 - d2 - d3 return math.Abs(delta) < eps // true if delta is less than a small tolerance }
var x, y float64 fx := func() float64 { return -(A*x + C) / B } fy := func() float64 { return -(B*y + C) / A } rxy := func() { if !segment || within(x, y) { res = append(res, point{x, y}) } }
if d == 0 { // line is tangent to circle, so just one intersect at most if bnz { x = -b / (2 * a) y = fx() rxy() } else { y = -b / (2 * a) x = fy() rxy() } } else { // two intersects at most d = math.Sqrt(d) if bnz { x = (-b + d) / (2 * a) y = fx() rxy() x = (-b - d) / (2 * a) y = fx() rxy() } else { y = (-b + d) / (2 * a) x = fy() rxy() y = (-b - d) / (2 * a) x = fy() rxy() } } return res
}
func main() {
cp := point{3, -5} r := 3.0 fmt.Println("The intersection points (if any) between:") fmt.Println("\n A circle, center (3, -5) with radius 3, and:") fmt.Println("\n a line containing the points (-10, 11) and (10, -9) is/are:") fmt.Println(" ", intersects(point{-10, 11}, point{10, -9}, cp, r, false)) fmt.Println("\n a segment starting at (-10, 11) and ending at (-11, 12) is/are") fmt.Println(" ", intersects(point{-10, 11}, point{-11, 12}, cp, r, true)) fmt.Println("\n a horizontal line containing the points (3, -2) and (7, -2) is/are:") fmt.Println(" ", intersects(point{3, -2}, point{7, -2}, cp, r, false)) cp = point{0, 0} r = 4.0 fmt.Println("\n A circle, center (0, 0) with radius 4, and:") fmt.Println("\n a vertical line containing the points (0, -3) and (0, 6) is/are:") fmt.Println(" ", intersects(point{0, -3}, point{0, 6}, cp, r, false)) fmt.Println("\n a vertical segment starting at (0, -3) and ending at (0, 6) is/are:") fmt.Println(" ", intersects(point{0, -3}, point{0, 6}, cp, r, true)) cp = point{4, 2} r = 5.0 fmt.Println("\n A circle, center (4, 2) with radius 5, and:") fmt.Println("\n a line containing the points (6, 3) and (10, 7) is/are:") fmt.Println(" ", intersects(point{6, 3}, point{10, 7}, cp, r, false)) fmt.Println("\n a segment starting at (7, 4) and ending at (11, 8) is/are:") fmt.Println(" ", intersects(point{7, 4}, point{11, 8}, cp, r, true))
}</lang>
- Output:
The intersection points (if any) between: A circle, center (3, -5) with radius 3, and: a line containing the points (-10, 11) and (10, -9) is/are: [(6, -5) (3, -2)] a segment starting at (-10, 11) and ending at (-11, 12) is/are [] a horizontal line containing the points (3, -2) and (7, -2) is/are: [(3, -2)] A circle, center (0, 0) with radius 4, and: a vertical line containing the points (0, -3) and (0, 6) is/are: [(0, 4) (0, -4)] a vertical segment starting at (0, -3) and ending at (0, 6) is/are: [(0, 4)] A circle, center (4, 2) with radius 5, and: a line containing the points (6, 3) and (10, 7) is/are: [(8, 5) (1, -2)] a segment starting at (7, 4) and ending at (11, 8) is/are: [(8, 5)]
Haskell
<lang haskell>import Data.Tuple.Curry
main :: IO () main =
mapM_ putStrLn $ concatMap (("" :) . uncurryN task) [ ((-10, 11), (10, -9), ((3, -5), 3)) , ((-10, 11), (-11, 12), ((3, -5), 3)) , ((3, -2), (7, -2), ((3, -5), 3)) , ((3, -2), (7, -2), ((0, 0), 4)) , ((0, -3), (0, 6), ((0, 0), 4)) , ((6, 3), (10, 7), ((4, 2), 5)) , ((7, 4), (11, 18), ((4, 2), 5)) ]
task :: (Double, Double)
-> (Double, Double) -> ((Double, Double), Double) -> [String]
task pt1 pt2 circle@(pt3@(a3, b3), r) = [line, segment]
where xs = map fun $ lineCircleIntersection pt1 pt2 circle ys = map fun $ segmentCircleIntersection pt1 pt2 circle to x = (fromIntegral . round $ 100 * x) / 100 fun (x, y) = (to x, to y) yo = show . fun start = "Intersection: Circle " ++ yo pt3 ++ " " ++ show (to r) ++ " and " end = yo pt1 ++ " " ++ yo pt2 ++ ": " line = start ++ "Line " ++ end ++ show xs segment = start ++ "Segment " ++ end ++ show ys
segmentCircleIntersection
:: (Double, Double) -> (Double, Double) -> ((Double, Double), Double) -> [(Double, Double)]
segmentCircleIntersection pt1 pt2 circle =
filter (go p1 p2) $ lineCircleIntersection pt1 pt2 circle where [p1, p2] | pt1 < pt2 = [pt1, pt2] | otherwise = [pt2, pt1] go (x, y) (u, v) (i, j) | x == u = y <= j && j <= v | otherwise = x <= i && i <= u
lineCircleIntersection
:: (Double, Double) -> (Double, Double) -> ((Double, Double), Double) -> [(Double, Double)]
lineCircleIntersection (a1, b1) (a2, b2) ((a3, b3), r) = go delta
where (x1, x2) = (a1 - a3, a2 - a3) (y1, y2) = (b1 - b3, b2 - b3) (dx, dy) = (x2 - x1, y2 - y1) drdr = dx * dx + dy * dy d = x1 * y2 - x2 * y1 delta = r * r * drdr - d * d sqrtDelta = sqrt delta (sgnDy, absDy) = (sgn dy, abs dy) u1 = (d * dy + sgnDy * dx * sqrtDelta) / drdr u2 = (d * dy - sgnDy * dx * sqrtDelta) / drdr v1 = (-d * dx + absDy * sqrtDelta) / drdr v2 = (-d * dx - absDy * sqrtDelta) / drdr go x | 0 > x = [] | 0 == x = [(u1 + a3, v1 + b3)] | otherwise = [(u1 + a3, v1 + b3), (u2 + a3, v2 + b3)]
sgn :: Double -> Double sgn x
| 0 > x = -1 | otherwise = 1</lang>
- Output:
Intersection: Circle (3.0,-5.0) 3.0 and Line (-10.0,11.0) (10.0,-9.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (-10.0,11.0) (10.0,-9.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Line (-10.0,11.0) (-11.0,12.0): [(3.0,-2.0),(6.0,-5.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (-10.0,11.0) (-11.0,12.0): [] Intersection: Circle (3.0,-5.0) 3.0 and Line (3.0,-2.0) (7.0,-2.0): [(3.0,-2.0)] Intersection: Circle (3.0,-5.0) 3.0 and Segment (3.0,-2.0) (7.0,-2.0): [(3.0,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Line (3.0,-2.0) (7.0,-2.0): [(3.46,-2.0),(-3.46,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Segment (3.0,-2.0) (7.0,-2.0): [(3.46,-2.0)] Intersection: Circle (0.0,0.0) 4.0 and Line (0.0,-3.0) (0.0,6.0): [(0.0,4.0),(0.0,-4.0)] Intersection: Circle (0.0,0.0) 4.0 and Segment (0.0,-3.0) (0.0,6.0): [(0.0,4.0)] Intersection: Circle (4.0,2.0) 5.0 and Line (6.0,3.0) (10.0,7.0): [(8.0,5.0),(1.0,-2.0)] Intersection: Circle (4.0,2.0) 5.0 and Segment (6.0,3.0) (10.0,7.0): [(8.0,5.0)] Intersection: Circle (4.0,2.0) 5.0 and Line (7.0,4.0) (11.0,18.0): [(7.46,5.61),(5.03,-2.89)] Intersection: Circle (4.0,2.0) 5.0 and Segment (7.0,4.0) (11.0,18.0): [(7.46,5.61)]
Perl
<lang perl>use strict; use warnings; use feature 'say'; use List::Util 'sum';
sub find_intersection {
my($P1, $P2, $center, $radius) = @_; my @d = ($$P2[0] - $$P1[0], $$P2[1] - $$P1[1]); my @f = ($$P1[0] - $$center[0], $$P1[1] - $$center[1]); my $a = sum map { $_**2 } @d; my $b = 2 * ($f[0]*$d[0] + $f[1]*$d[1]); my $c = sum(map { $_**2 } @f) - $radius**2; my $D = $b**2 - 4*$a*$c;
return unless $D >= 0; my ($t1, $t2) = ( (-$b - sqrt $D) / (2*$a), (-$b + sqrt $D) / (2*$a) ); return unless $t1 >= 0 and $t1 <= 1 or $t2 >= 0 and $t2 <= 1;
my ($dx, $dy) = ($$P2[0] - $$P1[0], $$P2[1] - $$P1[1]); return ([$dx*$t1 + $$P1[0], $dy*$t1 + $$P1[1]], [$dx*$t2 + $$P1[0], $dy*$t2 + $$P1[1]])
}
my @data = (
[ [-10, 11], [ 10, -9], [3, -5], 3 ], [ [-10, 11], [-11, 12], [3, -5], 3 ], [ [ 3, -2], [ 7, -2], [3, -5], 3 ], [ [ 3, -2], [ 7, -2], [0, 0], 4 ], [ [ 0, -3], [ 0, 6], [0, 0], 4 ], [ [ 6, 3], [ 10, 7], [4, 2], 5 ], [ [ 7, 4], [ 11, 18], [4, 2], 5 ],
);
sub rnd { map { sprintf('%.2f', $_) =~ s/\.00//r } @_ }
for my $d (@data) {
my @solution = find_intersection @$d[0] , @$d[1] , @$d[2], @$d[3]; say 'For input: ' . join ', ', (map { '('. join(',', @$_) .')' } @$d[0,1,2]), @$d[3]; say 'Solutions: ' . (@solution > 1 ? join ', ', map { '('. join(',', rnd @$_) .')' } @solution : 'None'); say ;
}</lang>
- Output:
For input: (-10,11), (10,-9), (3,-5), 3 Solutions: (3,-2), (6,-5) For input: (-10,11), (-11,12), (3,-5), 3 Solutions: None For input: (3,-2), (7,-2), (3,-5), 3 Solutions: (3,-2), (3,-2) For input: (3,-2), (7,-2), (0,0), 4 Solutions: (-3.46,-2), (3.46,-2) For input: (0,-3), (0,6), (0,0), 4 Solutions: (0,-4), (0,4) For input: (6,3), (10,7), (4,2), 5 Solutions: (1,-2), (8,5) For input: (7,4), (11,18), (4,2), 5 Solutions: (5.03,-2.89), (7.46,5.61)
Phix
<lang Phix>constant epsilon = 1e-14 -- say atom cx, cy, r, x1, y1, x2, y2
function sq(atom x) return x*x end function
function within(atom x, y)
-- -- checks whether a point is within a segment -- ie: <-------d1-------> -- <--d2--><---d3---> -- within, d2+d3 ~= d1 -- x1,y1^ ^x,y ^x2,y2 -- vs: -- <-d2-> -- <-----------d3---------> -- not "", d2+d3 > d1 -- ^x,y - and obviously ditto when x,y is (say) out here^ -- -- (obviously only works when x,y is on the same line as x1,y1 to x2,y2) -- atom d1 := sqrt(sq(x2-x1) + sq(y2-y1)), -- distance between end-points d2 := sqrt(sq(x -x1) + sq(y -y1)), -- distance from point to one end d3 := sqrt(sq(x2-x ) + sq(y2-y )), -- distance from point to other end delta := (d2 + d3) - d1 return abs(delta) < epsilon -- true if delta is less than a small tolerance
end function
function pf(atom x,y)
return sprintf("(%g,%g)",{x,y})
end function
function intersects(bool bSegment) -- -- Returns the intersection points (if any) of a circle, center (cx,cy) with radius r, -- and line containing the points (x1,y1) and (x2,y2) being either infinite or limited -- to the segment drawn between those points. --
sequence res = {} atom A = y2 - y1, sqA = sq(A), B = x1 - x2, sqB = sq(B), C = x2*y1 - x1*y2, sqC = sq(C), sqr = r*r-cx*cx-cy*cy, a := sqA + sqB, b, c bool bDivA = false if abs(B)<epsilon then -- B is zero or close to it b = 2 * (B*C + A*B*cx - sqA*cy) c = sqC + 2*A*C*cx - sqA*sqr bDivA = true -- (and later divide by A instead!) else b = 2 * (A*C + A*B*cy - sqB*cx) c = sqC + 2*B*C*cy - sqB*sqr end if atom d := b*b - 4*a*c -- discriminant if d>=0 then -- (-ve means line & circle do not intersect) d = sqrt(d) atom ux,uy, vx,vy if bDivA then {uy,vy} = sq_div(sq_sub({+d,-d},b),2*a) {ux,vx} = sq_div(sq_sub(sq_mul(-B,{uy,vy}),C),A) else {ux,vx} = sq_div(sq_sub({+d,-d},b),2*a) {uy,vy} = sq_div(sq_sub(sq_mul(-A,{ux,vx}),C),B) end if if not bSegment or within(ux,uy) then res = append(res,pf(ux,uy)) end if if d!=0 and (not bSegment or within(vx,vy)) then res = append(res,pf(vx,vy)) end if end if return res
end function
-- cx cy r x1 y1 x2 y2 bSegment constant tests = {{3,-5,3,{{-10,11, 10,-9,false},
{-10,11,-11,12,true}, { 3,-2, 7,-2,false}}}, {0, 0,4,{{ 0,-3, 0, 6,false}, { 0,-3, 0, 6,true}}}, {4, 2,5,{{ 6, 3, 10, 7,false}, { 7, 4, 11, 8,true}}}}
for t=1 to length(tests) do
{cx, cy, r, sequence lines} = tests[t] string circle = sprintf("Circle at %s radius %d",{pf(cx,cy),r}) for l=1 to length(lines) do {x1, y1, x2, y2, bool bSegment} = lines[l] sequence res = intersects(bSegment) string ls = iff(bSegment?"segment":" line"), at = iff(length(res)?"intersect at "&join(res," and ") :"do not intersect") printf(1,"%s and %s %s to %s %s.\n",{circle,ls,pf(x1,y1),pf(x2,y2),at}) circle = repeat(' ',length(circle)) end for
end for</lang>
- Output:
Circle at (3,-5) radius 3 and line (-10,11) to (10,-9) intersect at (6,-5) and (3,-2). and segment (-10,11) to (-11,12) do not intersect. and line (3,-2) to (7,-2) intersect at (3,-2). Circle at (0,0) radius 4 and line (0,-3) to (0,6) intersect at (0,4) and (0,-4). and segment (0,-3) to (0,6) intersect at (0,4). Circle at (4,2) radius 5 and line (6,3) to (10,7) intersect at (8,5) and (1,-2). and segment (7,4) to (11,8) intersect at (8,5).
Raku
(formerly Perl 6) Extend solution space to 3D. Reference: this SO question and answers <lang perl6>sub LineCircularOBJintersection(@P1, @P2, @Centre, \Radius) {
my @d = @P2 »-« @P1 ; # d my @f = @P1 »-« @Centre ; # c
my \a = [+] @d»²; # d dot d my \b = 2 * ([+] @f »*« @d); # 2 * f dot d my \c = ([+] @f»²) - Radius²; # f dot f - r² my \Δ = b²-(4*a*c); # discriminant
if (Δ < 0) { return []; } else { my (\t1,\t2) = (-b - Δ.sqrt)/(2*a), (-b + Δ.sqrt)/(2*a); if 0 ≤ t1|t2 ≤ 1 { return @P1 »+« ( @P2 »-« @P1 ) »*» t1, @P1 »+« ( @P2 »-« @P1 ) »*» t2 } else { return [] } }
}
my \DATA = [
[ <-10 11>, < 10 -9>, <3 -5>, 3 ], [ <-10 11>, <-11 12>, <3 -5>, 3 ], [ < 3 -2>, < 7 -2>, <3 -5>, 3 ], [ < 3 -2>, < 7 -2>, <0 0>, 4 ], [ < 0 -3>, < 0 6>, <0 0>, 4 ], [ < 6 3>, < 10 7>, <4 2>, 5 ], [ < 7 4>, < 11 18>, <4 2>, 5 ], [ <5 2 −2.26 >, <0.77 2 4>, <1 4 0>, 4 ]
];
for DATA {
my @solution = LineCircularOBJintersection $_[0] , $_[1] , $_[2], $_[3]; say "For data set: ", $_; say "Solution(s) is/are: ", @solution.Bool ?? @solution !! "None";
}</lang>
- Output:
For data set: [(-10 11) (10 -9) (3 -5) 3] Solution(s) is/are: [(3 -2) (6 -5)] For data set: [(-10 11) (-11 12) (3 -5) 3] Solution(s) is/are: None For data set: [(3 -2) (7 -2) (3 -5) 3] Solution(s) is/are: [(3 -2) (3 -2)] For data set: [(3 -2) (7 -2) (0 0) 4] Solution(s) is/are: [(-3.4641016151377544 -2) (3.4641016151377544 -2)] For data set: [(0 -3) (0 6) (0 0) 4] Solution(s) is/are: [(0 -4) (0 4)] For data set: [(6 3) (10 7) (4 2) 5] Solution(s) is/are: [(1 -2) (8 5)] For data set: [(7 4) (11 18) (4 2) 5] Solution(s) is/are: [(5.030680985703315 -2.892616550038399) (7.459885052032535 5.60959768211387)] For data set: [(5 2 −2.26) (0.77 2 4) (1 4 0) 4] Solution(s) is/are: [(4.2615520237084015 2 -1.1671668246843006) (1.13386504516801 2 3.461514141193441)]
REXX
The formulae used for this REXX version were taken from the MathWorld webpage: circle line intersection. <lang rexx>/*REXX program calculates where (or if) a line intersects (or tengents) a cirle. */ /*───────────────────────────────────── line= x1,y1 x2,y2; circle is at 0,0, radius=r*/ parse arg x1 y1 x2 y2 cx cy r . /*obtain optional arguments from the CL*/ if x1== | x1=="," then x1= 0 /*Not specified? Then use the default.*/ if y1== | y1=="," then y1= -3 /* " " " " " " */ if x2== | x2=="," then x2= 0 /* " " " " " " */ if y2== | y2=="," then y2= 6 /* " " " " " " */ if cx== | cx=="," then cx= 0 /* " " " " " " */ if cy== | cy=="," then cy= 0 /* " " " " " " */ if r == | r =="," then r = 4 /* " " " " " " */ x_1= x1; x1= x1 + cx x_2= x2; x2= x2 + cx y_1= y1; y1= y1 + cy y_2= y2; y2= y2 + cy dx= x2 - x1; dy= y2 - y1 dr= sqrt(dx**2 + dy**2)
D= x1 * y2 - x2 * y1
ix1= ( D * dy + sgn(dy) * dx * sqrt(r**2 * dr**2 - D**2) ) / dr**2 ix2= ( D * dy - sgn(dy) * dx * sqrt(r**2 * dr**2 - D**2) ) / dr**2 iy1= (-D * dx + abs(dy) * sqrt(r**2 * dr**2 - D**2) ) / dr**2 iy2= (-D * dx - abs(dy) * sqrt(r**2 * dr**2 - D**2) ) / dr**2 incidence= (r**2 * dr**2 - D**2) / 1 say 'incidence=' incidence @potla= 'points on the line are: ' if incidence<0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line doesn't intersect the circle with radius: " r end
if incidence=0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line is tangent to circle with radius: " r end
if incidence>0 then do
say @potla ' ('||x_1","y_1') and ('||x_2","y_2') are: ' ix1","iy1 say "The line is secant to circle with radius: " r end
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ sgn: procedure; if arg(1)<0 then return -1; return 1 /*──────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2 do j=0 while h>9; m.j= h; h= h%2 + 1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g= (g+x/g) *.5; end /*k*/ numeric digits d; return g/1</lang>
- output when using the default inputs:
incidence= 1296 points on the line are: (0,-3) and (0,6) are: 0,4 The line is secant to circle with radius: 4
zkl
<lang zkl>const EPS=1e-14; // a close-ness to zero
// p1,p2 are (x,y), circle is ( (x,y),r )
fcn intersectLineCircle(p1,p2, circle, segment=False) // assume line {
cx,cy := circle[0].apply("toFloat"); r := circle[1].toFloat(); x1,y1 := p1.apply("toFloat"); x2,y2 := p2.apply("toFloat"); A,B,C,a := (y2 - y1), (x1 - x2), (x2*y1 - x1*y2), (A*A + B*B); b,c,bnz := 0,0,True; if(B.closeTo(0,EPS)){ // B is zero or close to it b = 2.0 * (B*C + A*B*cx - A*A*cy); c = C*C + 2.0*A*C*cx - A*A*(r*r - cx*cx - cy*cy); bnz = False }else{ b = 2.0*( A*C + A*B*cy - B*B*cx ); c = C*C + 2.0*B*C*cy - B*B*( r*r - cx*cx - cy*cy ); } d := b*b - 4.0*a*c; // discriminant if(d<0.0){ // no real solution? zero --> one solution if (d>-0.005) d=0.0; // close enough to zero else return(T); // no intersection } d=d.sqrt();
reg ux,uy, vx,vy; if(bnz){ ux,vx = (-b + d) / (2*a), (-b - d) / (2*a); uy,vy = -(A*ux + C) / B, -(A*vx + C) / B; }else{ uy,vy = (-b + d) / (2*a), (-b - d) / (2*a); ux,vx = -(B*uy + C) / A, -(B*vy + C) / A; }
if(segment){ within:='wrap(x,y){ // is (x,y) within segment p1 p2?
d1:=( (x2 - x1).pow(2) + (y2 - y1).pow(2) ).sqrt(); d2:=( (x - x1).pow(2) + (y - y1).pow(2) ).sqrt(); d3:=( (x2 - x) .pow(2) + (y2 - y) .pow(2) ).sqrt(); (d1 - d2 - d3).closeTo(0,EPS);
};
i1,i2 := within(ux,uy), within(vx,vy); if(d==0) return(if(i1) T(ux,uy) else T); return(T( i1 and T(ux,uy), i2 and T(vx,vy) ).filter()) }
if(d==0) return( T( T(ux,uy) ) ); return( T(ux,uy), T(vx,vy) )
}</lang> <lang zkl>circle:=T( T(3,-5),3 ); p1,p2 := T(-10,11), T( 10,-9); println("Circle @ ",circle); lcpp(p1,p2,circle); p2:=T(-11,12); lcpp(p1,p2,circle,True); p1,p2 := T(3,-2), T(7,-2); lcpp(p1,p2,circle);
circle:=T( T(0,0),4 ); p1,p2 := T(0,-3), T(0,6); println("\nCircle @ ",circle); lcpp(p1,p2,circle); lcpp(p1,p2,circle,True);
circle:=T( T(4,2),5 ); p1,p2 := T(6,3), T(10,7); println("\nCircle @ ",circle); lcpp(p1,p2,circle); p1,p2 := T(7,4), T(11,8); lcpp(p1,p2,circle,True);
fcn lcpp(p1,p2,circle,segment=False){
println(" %s %s -- %s intersects at %s" .fmt(segment and "Segment" or "Line ", p1,p2,intersectLineCircle(p1,p2, circle,segment)));
}</lang>
- Output:
Circle @ L(L(3,-5),3) Line L(-10,11) -- L(10,-9) intersects at L(L(6,-5),L(3,-2)) Segment L(-10,11) -- L(-11,12) intersects at L() Line L(3,-2) -- L(7,-2) intersects at L(L(3,-2)) Circle @ L(L(0,0),4) Line L(0,-3) -- L(0,6) intersects at L(L(0,4),L(0,-4)) Segment L(0,-3) -- L(0,6) intersects at L(L(0,4)) Circle @ L(L(4,2),5) Line L(6,3) -- L(10,7) intersects at L(L(8,5),L(1,-2)) Segment L(7,4) -- L(11,8) intersects at L(L(8,5))