Circles of given radius through two points
Given two points on a plane and a radius, usually two circles of given radius can be drawn through the points.
- Exceptions
- r==0.0 should be treated as never describing circles (except in the case where the points are coincident).
- If the points are coincident then an infinite number of circles with the point on their circumference can be drawn, unless r==0.0 as well which then collapses the circles to a point.
- If the points form a diameter then return two identical circles or return a single circle, according to which is the most natural mechanism for the implementation language.
- If the points are too far apart then no circles can be drawn.
- Task detail
- Write a function/subroutine/method/... that takes two points and a radius and returns the two circles through those points, or some indication of special cases where two, possibly equal, circles cannot be returned.
- Show here the output for the following inputs:
p1 p2 r 0.1234, 0.9876 0.8765, 0.2345 2.0 0.0000, 2.0000 0.0000, 0.0000 1.0 0.1234, 0.9876 0.1234, 0.9876 2.0 0.1234, 0.9876 0.8765, 0.2345 0.5 0.1234, 0.9876 0.1234, 0.9876 0.0
- Ref
- Finding the Center of a Circle from 2 Points and Radius from Math forum @ Drexel
D
<lang d>import std.stdio, std.typecons, std.math;
class ValueException : Exception {
this(string msg_) { super(msg_); }
}
struct V2 { double x, y; } struct Circle { double x, y, r; }
/**Following explanation at: http://mathforum.org/library/drmath/view/53027.html
- /
Tuple!(Circle, Circle) circlesFromTwoPointsAndRadius(in V2 p1, in V2 p2, in double r) pure in {
assert(r >= 0, "radius can't be negative");
} body {
enum nBits = 40;
if (r.abs < (1.0 / (2.0 ^^ nBits)))
throw new ValueException("radius of zero");
if (feqrel(cast()p1.x, cast()p2.x) >= nBits &&
feqrel(cast()p1.y, cast()p2.y) >= nBits)
throw new ValueException("coincident points give" ~
" infinite number of Circles");
// Delta x, delta y between points. immutable dx = p2.x - p1.x; immutable dy = p2.y - p1.y;
// Dist between points.
immutable q = sqrt(dx ^^ 2 + dy ^^ 2);
if (q > 2.0 * r)
throw new ValueException("separation of points > diameter");
// Halfway point. immutable x3 = (p1.x + p2.x) / 2; immutable y3 = (p1.y + p2.y) / 2;
// Distance along the mirror line. immutable d = sqrt(r ^^ 2 - (q / 2) ^^ 2);
immutable c1 = Circle(x3 - d * dy / q, y3 + d * dx / q, r.abs); immutable c2 = Circle(x3 + d * dy / q, y3 - d * dx / q, r.abs); return typeof(return)(c1, c2);
}
void main() {
foreach (t; [tuple(V2(0.1234, 0.9876), V2(0.8765, 0.2345), 2.0),
tuple(V2(0.0000, 2.0000), V2(0.0000, 0.0000), 1.0),
tuple(V2(0.1234, 0.9876), V2(0.1234, 0.9876), 2.0),
tuple(V2(0.1234, 0.9876), V2(0.8765, 0.2345), 0.5),
tuple(V2(0.1234, 0.9876), V2(0.1234, 0.9876), 0.0)]) {
writefln("Through points:\n %s,\n %s\n and radius %f\n" ~
"You can construct the following circles:", t[]);
try {
writefln(" %s\n %s\n",
circlesFromTwoPointsAndRadius(t[])[]);
} catch (ValueException v)
writefln(" ERROR: %s\n", v.msg);
}
}</lang>
- Output:
Through points: V2(0.1234, 0.9876), V2(0.8765, 0.2345) and radius 2.000000 You can construct the following circles: Circle(1.86311, 1.97421, 2) Circle(-0.863212, -0.752112, 2) Through points: V2(0, 2), V2(0, 0) and radius 1.000000 You can construct the following circles: Circle(0, 1, 1) Circle(0, 1, 1) Through points: V2(0.1234, 0.9876), V2(0.1234, 0.9876) and radius 2.000000 You can construct the following circles: ERROR: coincident points give infinite number of Circles Through points: V2(0.1234, 0.9876), V2(0.8765, 0.2345) and radius 0.500000 You can construct the following circles: ERROR: separation of points > diameter Through points: V2(0.1234, 0.9876), V2(0.1234, 0.9876) and radius 0.000000 You can construct the following circles: ERROR: radius of zero
J
2D computations are often easier using the complex plane. <lang J> [INPUT =: _5]\0.1234, 0.9876 0.8765, 0.2345 2.0 0.0000, 2.0000 0.0000, 0.0000 1.0 0.1234, 0.9876 0.1234, 0.9876 2.0 0.1234, 0.9876 0.8765, 0.2345 0.5 0.1234, 0.9876 0.1234, 0.9876 0.0
average =: +/ % #
circles =: verb define"1
'P0 P1 R' =. (j./"1)_2[\y NB. Use complex plane
C =. P0 average@:, P1
BAD =: ":@:+. C
SEPARATION =. P0 |@- P1
if. 0 = SEPARATION do.
if. 0 = R do. 'Degenerate point at ' , BAD
else. 'Any center at a distance ' , (":R) , ' from ' , BAD , ' works.'
end.
elseif. SEPARATION (> +:) R do. 'No solutions.'
elseif. SEPARATION (= +:) R do. 'Duplicate solutions with center at ' , BAD
elseif. 1 do.
ORTHOGONAL_DISTANCE =. R * 1 o. _2 o. R %~ | C - P0
UNIT =: P1 *@:- P0
OFFSETS =: ORTHOGONAL_DISTANCE * UNIT * j. _1 1
C +.@:+ OFFSETS
end.
)
('x0 y0 x1 y1 r' ; 'center'),(;circles)"1 INPUT
┌───────────────────────────────┬────────────────────────────────────────────────────┐ │x0 y0 x1 y1 r │center │ ├───────────────────────────────┼────────────────────────────────────────────────────┤ │0.1234 0.9876 0.8765 0.2345 2 │_0.863212 _0.752112 │ │ │ 1.86311 1.97421 │ ├───────────────────────────────┼────────────────────────────────────────────────────┤ │0 2 0 0 1 │Duplicate solutions with center at 0 1 │ ├───────────────────────────────┼────────────────────────────────────────────────────┤ │0.1234 0.9876 0.1234 0.9876 2 │Any center at a distance 2 from 0.1234 0.9876 works.│ ├───────────────────────────────┼────────────────────────────────────────────────────┤ │0.1234 0.9876 0.8765 0.2345 0.5│No solutions. │ ├───────────────────────────────┼────────────────────────────────────────────────────┤ │0.1234 0.9876 0.1234 0.9876 0 │Degenerate point at 0.1234 0.9876 │ └───────────────────────────────┴────────────────────────────────────────────────────┘ </lang>
PL/I
<lang PL/I>twoci: Proc Options(main);
Dcl 1 *(5),
2 m1x Dec Float Init(0.1234, 0,0.1234,0.1234,0.1234),
2 m1y Dec Float Init(0.9876, 2,0.9876,0.9876,0.9876),
2 m2x Dec Float Init(0.8765, 0,0.1234,0.8765,0.1234),
2 m2y Dec Float Init(0.2345, 0,0.9876,0.2345,0.9876),
2 r Dec Float Init( 2, 1, 2,0.5 , 0);
Dcl i Bin Fixed(31);
Put Edit(' x1 y1 x2 y2 r '||
' cir1x cir1y cir2x cir2y')(Skip,a);
Put Edit(' ====== ====== ====== ====== = '||
' ====== ====== ====== ======')(Skip,a);
Do i=1 To 5;
Put Edit(m1x(i),m1y(i),m2x(i),m2y(i),r(i))
(Skip,4(f(7,4)),f(3));
Put Edit(twocircles(m1x(i),m1y(i),m2x(i),m2y(i),r(i)))(a);
End;
twoCircles: proc(m1x,m1y,m2x,m2y,r) Returns(Char(50) Var);
Dcl (m1x,m1y,m2x,m2y,r) Dec Float;
Dcl (cx,cy,bx,by,pb,x,y,x1,y1) Dec Float;
Dcl res Char(50) Var;
If r=0 then return(' radius of zero gives no circles.');
x=(m2x-m1x)/2;
y=(m2y-m1y)/2;
bx=m1x+x;
by=m1y+y;
pb=sqrt(x**2+y**2);
cx=(m2x-m1x)/2;
cy=(m2y-m1y)/2;
bx=m1x+x;
by=m1y+y;
pb=sqrt(x**2+y**2)
if pb=0 then return(' coincident points give infinite circles');
if pb>r then return(' points are too far apart for the given radius');
cb=sqrt(r**2-pb**2);
x1=y*cb/pb;
y1=x*cb/pb
Put String(res) Edit((bx-x1),(by+y1),(bx+x1),(by-y1))(4(f(8,4)));
Return(res);
End;
End;</lang>
Output:
x1 y1 x2 y2 r cir1x cir1y cir2x cir2y ====== ====== ====== ====== = ====== ====== ====== ====== 0.1234 0.9876 0.8765 0.2345 2 1.8631 1.9742 -0.8632 -0.7521 0.0000 2.0000 0.0000 0.0000 1 0.0000 1.0000 0.0000 1.0000 0.1234 0.9876 0.1234 0.9876 2 coincident points give infinite circles 0.1234 0.9876 0.8765 0.2345 1 points are too far apart for the given radius 0.1234 0.9876 0.1234 0.9876 0 radius of zero gives no circles.
Python
The function raises the ValueError exception for the special cases and usess try - except to catch these and extract the exception detail.
<lang python>from collections import namedtuple from math import sqrt
Pt = namedtuple('Pt', 'x, y') Circle = Cir = namedtuple('Circle', 'x, y, r')
def circles_from_p1p2r(p1, p2, r):
'Following explanation at http://mathforum.org/library/drmath/view/53027.html' if r == 0.0: raise ValueError('radius of zero') (x1, y1), (x2, y2) = p1, p2 if p1 == p2: raise ValueError('coincident points gives infinite number of Circles') # delta x, delta y between points dx, dy = x2 - x1, y2 - y1 # dist between points q = sqrt(dx**2 + dy**2) if q > 2.0*r: raise ValueError('separation of points > diameter') # halfway point x3, y3 = (x1+x2)/2, (y1+y2)/2 # distance along the mirror line d = sqrt(r**2-(q/2)**2) # One answer c1 = Cir(x = x3 - d*dy/q, y = y3 + d*dx/q, r = abs(r)) # The other answer c2 = Cir(x = x3 + d*dy/q, y = y3 - d*dx/q, r = abs(r)) return c1, c2
if __name__ == '__main__':
for p1, p2, r in [(Pt(0.1234, 0.9876), Pt(0.8765, 0.2345), 2.0),
(Pt(0.0000, 2.0000), Pt(0.0000, 0.0000), 1.0),
(Pt(0.1234, 0.9876), Pt(0.1234, 0.9876), 2.0),
(Pt(0.1234, 0.9876), Pt(0.8765, 0.2345), 0.5),
(Pt(0.1234, 0.9876), Pt(0.1234, 0.9876), 0.0)]:
print('Through points:\n %r,\n %r\n and radius %f\nYou can construct the following circles:'
% (p1, p2, r))
try:
print(' %r\n %r\n' % circles_from_p1p2r(p1, p2, r))
except ValueError as v:
print(' ERROR: %s\n' % (v.args[0],))</lang>
- Output:
Through points: Pt(x=0.1234, y=0.9876), Pt(x=0.8765, y=0.2345) and radius 2.000000 You can construct the following circles: Circle(x=1.8631118016581893, y=1.974211801658189, r=2.0) Circle(x=-0.8632118016581896, y=-0.7521118016581892, r=2.0) Through points: Pt(x=0.0, y=2.0), Pt(x=0.0, y=0.0) and radius 1.000000 You can construct the following circles: Circle(x=0.0, y=1.0, r=1.0) Circle(x=0.0, y=1.0, r=1.0) Through points: Pt(x=0.1234, y=0.9876), Pt(x=0.1234, y=0.9876) and radius 2.000000 You can construct the following circles: ERROR: coincident points gives infinite number of Circles Through points: Pt(x=0.1234, y=0.9876), Pt(x=0.8765, y=0.2345) and radius 0.500000 You can construct the following circles: ERROR: separation of points > diameter Through points: Pt(x=0.1234, y=0.9876), Pt(x=0.1234, y=0.9876) and radius 0.000000 You can construct the following circles: ERROR: radius of zero
Racket
Using library `plot/utils` for simple vector operations.
<lang racket>
- lang racket
(require plot/utils)
(define (circle-centers p1 p2 r)
(when (zero? r) (err "zero radius.")) (when (equal? p1 p2) (err "the points coinside.")) ; the midle point (define m (v/ (v+ p1 p2) 2)) ; the vector connecting given points (define d (v/ (v- p1 p2) 2)) ; the distance between the center of the circle and the middle point (define ξ (- (sqr r) (vmag^2 d))) (when (negative? ξ) (err "given radius is less then the distance between points.")) ; the unit vector orthogonal to the delta (define n (vnormalize (orth d))) ; the shift along the direction orthogonal to the delta (define x (v* n (sqrt ξ))) (values (v+ m x) (v- m x)))
- error message
(define (err m) (error "Impossible to build a circle:" m))
- returns a vector which is orthogonal to the geven one
(define orth (match-lambda [(vector x y) (vector y (- x))])) </lang>
Testing
> (circle-centers #(0.1234 0.9876) #(0.8765 0.2345) 2.0) '#(1.8631118016581893 1.974211801658189) '#(-0.8632118016581896 -0.7521118016581892) > (circle-centers #(0.0000 2.0000) #(0.0000 0.0000) 1.0) '#(0.0 1.0) '#(0.0 1.0) > (circle-centers #(0.1234 0.9876) #(0.1234 0.9876) 2.0) . . Impossible to find a circle: "the points coinside." > (circle-centers #(0.1234 0.9876) #(0.8765 0.2345) 0.5) . . Impossible to find a circle: "given radius is less then the distance between points." > (circle-centers #(0.1234 0.9876) #(0.1234 0.9876) 0.0) . . Impossible to find a circle: "zero radius."
Drawing circles:
<lang racket> (require 2htdp/image)
(define/match (point v)
[{(vector x y)} (λ (s) (place-image (circle 2 "solid" "black") x y s))])
(define/match (circ v r)
[{(vector x y) r} (λ (s) (place-image (circle r "outline" "red") x y s))])
(define p1 #(40 50)) (define p2 #(60 30)) (define r 20) (define-values (x1 x2) (circle-centers p1 p2 r))
((compose (point p1) (point p2) (circ x1 r) (circ x2 r))
(empty-scene 100 100))
</lang>
REXX
<lang rexx>/*REXX pgm finds 2 circles with a specific radius given two (X,Y) points*/ @. = @.1= 0.1234 0.9876 0.8765 0.2345 2 @.2= 0.0000 2.0000 0.0000 0.0000 1 @.3= 0.1234 0.9876 0.1234 0.9876 2 @.4= 0.1234 0.9876 0.8765 0.2345 0.5 @.5= 0.1234 0.9876 0.1234 0.9876 0 say ' x1 y1 x2 y2 radius cir1x cir1y cir2x cir2y' say ' ──────── ──────── ──────── ──────── ────── ──────── ──────── ──────── ────────'
do j=1 while @.j\== /*process all given points&radius*/
do k=1 for 4; w.k=f(word(@.j,k)); end /*k*/ /*format num.*/
say w.1 w.2 w.3 w.4 center(word(@.j,5),9) "───► " twoCircles(@.j)
end /*j*/
exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────F subroutine────────────────────────*/ f: return right(format(arg(1),,4),9) /*format a number with 4 dec dig.*/ /*──────────────────────────────────SQRT subroutine─────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits();numeric digits 11
g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end
numeric digits d; return g/1
.sqrtGuess: if x<0 then call sqrtErr; numeric form; m.=11; p=d+d%4+2
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2
/*──────────────────────────────────twoCircles subroutine───────────────*/ twoCircles: procedure; parse arg px py qx qy r . if r=0 then return 'radius of zero gives no circles.' x=(qx-px)/2; y=(qy-py)/2; bx=px+x; by=py+y; pb=sqrt(x**2+y**2) if pb=0 then return 'coincident points give infinite circles' if pb>r then return 'points are too far apart for the given radius' cb=sqrt(r**2-pb**2); x1=y*cb/pb; y1=x*cb/pb return f(bx-x1) f(by+y1) f(bx+x1) f(by-y1)</lang> output when using the default input(s):
x1 y1 x2 y2 radius cir1x cir1y cir2x cir2y ──────── ──────── ──────── ──────── ────── ──────── ──────── ──────── ──────── 0.1234 0.9876 0.8765 0.2345 2 ───► 1.8631 1.9742 -0.8632 -0.7521 0.0000 2.0000 0.0000 0.0000 1 ───► 0.0000 1.0000 0.0000 1.0000 0.1234 0.9876 0.1234 0.9876 2 ───► coincident points give infinite circles 0.1234 0.9876 0.8765 0.2345 0.5 ───► points are too far apart for the given radius 0.1234 0.9876 0.1234 0.9876 0 ───► radius of zero gives no circles.
Tcl
<lang tcl>proc findCircles {p1 p2 r} {
lassign $p1 x1 y1
lassign $p2 x2 y2
# Special case: coincident & zero size
if {$x1 == $x2 && $y1 == $y2 && $r == 0.0} {
return [list [list $x1 $y1 0.0]]
}
if {$r <= 0.0} {
error "radius must be positive for sane results"
}
if {$x1 == $x2 && $y1 == $y2} {
error "no sane solution: points are coincident"
}
# Calculate distance apart and separation vector
set dx [expr {$x2 - $x1}]
set dy [expr {$y2 - $y1}]
set q [expr {hypot($dx, $dy)}]
if {$q > 2*$r} {
error "no solution: points are further apart than required diameter"
}
# Calculate midpoint
set x3 [expr {($x1+$x2)/2.0}]
set y3 [expr {($y1+$y2)/2.0}]
# Fractional distance along the mirror line
set f [expr {($r**2 - ($q/2.0)**2)**0.5 / $q}]
# The two answers
set c1 [list [expr {$x3 - $f*$dy}] [expr {$y3 + $f*$dx}] $r]
set c2 [list [expr {$x3 + $f*$dy}] [expr {$y3 - $f*$dx}] $r]
return [list $c1 $c2]
}</lang> Demonstrating: <lang tcl>foreach {p1 p2 r} {
{0.1234 0.9876} {0.8765 0.2345} 2.0
{0.0000 2.0000} {0.0000 0.0000} 1.0
{0.1234 0.9876} {0.1234 0.9876} 2.0
{0.1234 0.9876} {0.8765 0.2345} 0.5
{0.1234 0.9876} {0.1234 0.9876} 0.0
} {
puts "p1:([join $p1 {, }]) p2:([join $p2 {, }]) r:$r =>"
if {[catch {
foreach c [findCircles $p1 $p2 $r] { puts "\tCircle:([join $c {, }])" }
} msg]} {
puts "\tERROR: $msg"
}
}</lang>
- Output:
p1:(0.1234, 0.9876) p2:(0.8765, 0.2345) r:2.0 => Circle:(1.863111801658189, 1.974211801658189, 2.0) Circle:(-0.8632118016581891, -0.752111801658189, 2.0) p1:(0.0000, 2.0000) p2:(0.0000, 0.0000) r:1.0 => Circle:(0.0, 1.0, 1.0) Circle:(0.0, 1.0, 1.0) p1:(0.1234, 0.9876) p2:(0.1234, 0.9876) r:2.0 => ERROR: no sane solution: points are coincident p1:(0.1234, 0.9876) p2:(0.8765, 0.2345) r:0.5 => ERROR: no solution: points are further apart than required diameter p1:(0.1234, 0.9876) p2:(0.1234, 0.9876) r:0.0 => Circle:(0.1234, 0.9876, 0.0)
XPL0
An easy way to solve this is to translate the coordinates so that one point is at the origin. Then rotate the coordinate frame so that the second point is on the X-axis. The circles' X coordinate is then half the distance to the second point. The circles' Y coordinates are easily seen as +/-sqrt(radius^2 - circleX^2). Now undo the rotation and translation. The method used here is a streamlining of these steps.
<lang XPL0>include c:\cxpl\codes;
proc Circles; real Data; \Show centers of circles, given points P & Q and radius real Px, Py, Qx, Qy, R, X, Y, X1, Y1, Bx, By, PB, CB; [Px:= Data(0); Py:= Data(1); Qx:= Data(2); Qy:= Data(3); R:= Data(4); if R = 0.0 then [Text(0, "Radius = zero gives no circles^M^J"); return]; X:= (Qx-Px)/2.0; Y:= (Qy-Py)/2.0; Bx:= Px+X; By:= Py+Y; PB:= sqrt(X*X + Y*Y); if PB = 0.0 then [Text(0, "Coincident points give infinite circles^M^J"); return]; if PB > R then [Text(0, "Points are too far apart for radius^M^J"); return]; CB:= sqrt(R*R - PB*PB); X1:= Y*CB/PB; Y1:= X*CB/PB; RlOut(0, Bx-X1); ChOut(0, ^,); RlOut(0, By+Y1); ChOut(0, 9\tab\); RlOut(0, Bx+X1); ChOut(0, ^,); RlOut(0, By-Y1); CrLf(0); ];
real Tbl; int I; [Tbl:=[[0.1234, 0.9876, 0.8765, 0.2345, 2.0],
[0.0000, 2.0000, 0.0000, 0.0000, 1.0],
[0.1234, 0.9876, 0.1234, 0.9876, 2.0],
[0.1234, 0.9876, 0.8765, 0.2345, 0.5],
[0.1234, 0.9876, 0.1234, 0.9876, 0.0]];
for I:= 0 to 4 do Circles(Tbl(I)); ]</lang>
- Output:
1.86311, 1.97421 -0.86321, -0.75211
0.00000, 1.00000 0.00000, 1.00000
Coincident points give infinite circles
Points are too far apart for radius
Radius = zero gives no circles