Centre and radius of a circle passing through 3 points in a plane
You are encouraged to solve this task according to the task description, using any language you may know.
Write a function which returns the centre and radius of a circle passing through three point in a plane. Demonstrate the function using the points (22.83,2.07) (14.39,30.24) and (33.65,17.31)
ALGOL 68
Follows the lines of the C++ code at geeksforgeeks.org.
BEGIN # find the centre and radius of a circle through 3 points #
# follows the lines of the C++ code at #
# https://www.geeksforgeeks.org/equation-of-circle-when-three-points-on-the-circle-are-given/ #
MODE POINT = STRUCT( REAL x, y );
MODE CIRCLE = STRUCT( POINT centre, REAL radius );
# returns the circle that passes through p1, p2 and p3 #
PROC find circle = ( POINT p1, p2, p3 )CIRCLE:
BEGIN
REAL x1 = x OF p1, y1 = y OF p1, x2 = x OF p2, y2 = y OF p2, x3 = x OF p3, y3 = y OF p3;
REAL x12 = x1 - x2
, x13 = x1 - x3
, y12 = y1 - y2
, y13 = y1 - y3
, y31 = y3 - y1
, y21 = y2 - y1
, x31 = x3 - x1
, x21 = x2 - x1
;
REAL sx13 = x1^2 - x3^2
, sy13 = y1^2 - y3^2
, sx21 = x2^2 - x1^2
, sy21 = y2^2 - y1^2
;
REAL f = ( ( ( sx13 + sy13 ) * x12 )
+ ( ( sx21 + sy21 ) * x13 )
)
/ ( 2 * ( ( y31 * x12 ) - ( y21 * x13 ) ) )
, g = ( ( ( sx13 + sy13 ) * y12 )
+ ( ( sx21 + sy21 ) * y13 )
)
/ ( 2 * ( ( x31 * y12 ) - ( x21 * y13 ) ) )
;
REAL c = - (x1^2) - (y1^2) - ( 2 * g * x1 ) - ( 2 * f * y1 );
( ( -g, -f ), sqrt( g^2 + f^2 - c ) )
END # find circle # ;
CIRCLE c = find circle( ( 22.83, 2.07 ), ( 14.39, 30.24 ), ( 33.65, 17.31 ) );
print( ( "Centre = ( ", fixed( x OF centre OF c, -10, 6 ) ) );
print( ( ", ", fixed( y OF centre OF c, -10, 6 ) ) );
print( ( " ), Radius = ", fixed( radius OF c, -10, 6 ), newline ) )
END- Output:
Centre = ( 18.978516, 16.265411 ), Radius = 14.708624
BASIC
BASIC256
call findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
end
subroutine findCircle(x1, y1, x2, y2, x3, y3)
x12 = x1 - x2
x13 = x1 - x3
y12 = y1 - y2
y13 = y1 - y3
y31 = y3 - y1
y21 = y2 - y1
x31 = x3 - x1
x21 = x2 - x1
sx13 = x1 * x1 - x3 * x3
sy13 = y1 * y1 - y3 * y3
sx21 = x2 * x2 - x1 * x1
sy21 = y2 * y2 - y1 * y1
f = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
g = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
c = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
h = -g
k = -f
r = sqr(h * h + k * k - c)
print "Centre is at ("; h; ", "; k; ")"
print "Radius is "; r
print
print "Check radius as the distance between the centre and the first point:"
print sqr((22.83 - h) ^ 2 + (2.07 - k) ^ 2)
end subroutine
- Output:
Centre is at (18.9785156601, 16.2654107977) Radius is 14.708624 Check radius as the distance between the centre and the first point: 14.708624
Chipmunk Basic
100 sub findcircle(x1,y1,x2,y2,x3,y3)
110 x12 = x1-x2
120 x13 = x1-x3
130 y12 = y1-y2
140 y13 = y1-y3
150 y31 = y3-y1
160 y21 = y2-y1
170 x31 = x3-x1
180 x21 = x2-x1
190 '
200 sx13 = x1*x1-x3*x3
210 sy13 = y1*y1-y3*y3
220 sx21 = x2*x2-x1*x1
230 sy21 = y2*y2-y1*y1
240 '
250 f = (sx13*x12+sy13*x12+sx21*x13+sy21*x13)/(y31*x12-y21*x13)/2
260 g = (sx13*y12+sy13*y12+sx21*y13+sy21*y13)/(x31*y12-x21*y13)/2
270 '
280 c = -x1*x1-y1*y1-2*g*x1-2*f*y1
290 h = -g
300 k = -f
310 r = sqr(h*h+k*k-c)
320 '
330 print "Centre is at ( ";h;", ";k;")"
340 print "Radius is ";r
350 print
360 print "Check radius as the distance between the centre and the first point:"
370 print sqr((22.83-h)^2+(2.07-k)^2)
380 end sub
390 cls
400 findcircle(22.83,2.07,14.39,30.24,33.65,17.31)
410 end
Gambas
Sub findCircle(x1 As Float, y1 As Float, x2 As Float, y2 As Float, x3 As Float, y3 As Float)
Dim x12 As Float = x1 - x2
Dim x13 As Float = x1 - x3
Dim y12 As Float = y1 - y2
Dim y13 As Float = y1 - y3
Dim y31 As Float = y3 - y1
Dim y21 As Float = y2 - y1
Dim x31 As Float = x3 - x1
Dim x21 As Float = x2 - x1
Dim sx13 As Float = x1 * x1 - x3 * x3
Dim sy13 As Float = y1 * y1 - y3 * y3
Dim sx21 As Float = x2 * x2 - x1 * x1
Dim sy21 As Float = y2 * y2 - y1 * y1
Dim f As Float = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
Dim g As Float = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
Dim c As Float = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
Dim h As Float = -g
Dim k As Float = -f
Dim r As Float = Sqr(h * h + k * k - c)
Print "Centre is at (" & h & ", " & k & ")"
Print "Radius is "; r
Print
Print "Check radius as the distance between the centre and the first point:"
Print Sqr((22.83 - h) ^ 2 + (2.07 - k) ^ 2)
End Sub
Public Sub Main()
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
End
- Output:
Centre is at (18.9785156601488, 16.2654107977159) Radius is 14,7086239783342 Check radius as the distance between the centre and the first point: 14,7086239783342
PureBasic
Procedure.d findCircle(x1, y1, x2, y2, x3, y3)
Define.d x12, x13, y12, y13, y31, y21, x31, x21, sx13, sy13, sx21, sy21
Define.d f, g, c, h, k, r
x12 = x1 - x2
x13 = x1 - x3
y12 = y1 - y2
y13 = y1 - y3
y31 = y3 - y1
y21 = y2 - y1
x31 = x3 - x1
x21 = x2 - x1
sx13 = x1 * x1 - x3 * x3
sy13 = y1 * y1 - y3 * y3
sx21 = x2 * x2 - x1 * x1
sy21 = y2 * y2 - y1 * y1
f = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
g = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
c = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
h = -g
k = -f
r = Sqr(h * h + k * k - c)
PrintN("Centre is at (" + Str(h) + ", " + Str(k) + ")")
PrintN("Radius is " + Str(r))
PrintN("")
PrintN("Check radius as the distance between the centre and the first point: " + Str(Sqr(Pow((22.83 - h), 2) + Pow((2.07 - k), 2))))
EndProcedure
OpenConsole()
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
PrintN(#CRLF$ + "--- terminado, pulsa RETURN---"): Input()
PrintN(#CRLF$ + "Press ENTER to exit"): Input()
CloseConsole()
Yabasic
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
END
SUB findCircle(x1, y1, x2, y2, x3, y3)
x12 = x1 - x2
x13 = x1 - x3
y12 = y1 - y2
y13 = y1 - y3
y31 = y3 - y1
y21 = y2 - y1
x31 = x3 - x1
x21 = x2 - x1
sx13 = x1 * x1 - x3 * x3
sy13 = y1 * y1 - y3 * y3
sx21 = x2 * x2 - x1 * x1
sy21 = y2 * y2 - y1 * y1
f = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
g = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
c = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
h = -g
k = -f
r = SQRT(h * h + k * k - c)
PRINT "Centre is at (", h, ", ", k, ")"
PRINT "Radius is ", r
PRINT
PRINT "Check radius as the distance between the centre and the first point:"
PRINT SQRT((22.83 - h) ^ 2 + (2.07 - k) ^ 2)
END SUB
- Output:
Centre is at (18.9785, 16.2654) Radius is 14.7086 Check radius as the distance between the centre and the first point: 14.7086
C++
This follows the lines of the C++ code here.
#include <iostream>
#include <math.h>
using namespace std;
void findCircle(float x1, float y1, float x2, float y2, float x3, float y3) {
float x12 = x1 - x2;
float x13 = x1 - x3;
float y12 = y1 - y2;
float y13 = y1 - y3;
float y31 = y3 - y1;
float y21 = y2 - y1;
float x31 = x3 - x1;
float x21 = x2 - x1;
float sx13 = pow(x1, 2) - pow(x3, 2);
float sy13 = pow(y1, 2) - pow(y3, 2);
float sx21 = pow(x2, 2) - pow(x1, 2);
float sy21 = pow(y2, 2) - pow(y1, 2);
float f = ((sx13) * (x12) + (sy13) * (x12) + (sx21) * (x13) + (sy21) * (x13))
/ (2 * ((y31) * (x12) - (y21) * (x13)));
float g = ((sx13) * (y12) + (sy13) * (y12) + (sx21) * (y13) + (sy21) * (y13))
/ (2 * ((x31) * (y12) - (x21) * (y13)));
float c = -pow(x1, 2) - pow(y1, 2) - 2 * g * x1 - 2 * f * y1;
float h = -g;
float k = -f;
float r = sqrt(h * h + k * k - c);
cout << "Centre is at (" << h << ", " << k << ")" << endl;
cout << "Radius is " << r;
}
int main() {
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31);
return 0;
}
Dart
import 'dart:math';
void findCircle(
double x1, double y1, double x2, double y2, double x3, double y3) {
double x12 = x1 - x2;
double x13 = x1 - x3;
double y12 = y1 - y2;
double y13 = y1 - y3;
double y31 = y3 - y1;
double y21 = y2 - y1;
double x31 = x3 - x1;
double x21 = x2 - x1;
double sx13 = x1 * x1 - x3 * x3;
double sy13 = y1 * y1 - y3 * y3;
double sx21 = x2 * x2 - x1 * x1;
double sy21 = y2 * y2 - y1 * y1;
double f = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) /
(2 * (y31 * x12 - y21 * x13));
double g = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) /
(2 * (x31 * y12 - x21 * y13));
double c = -pow(x1, 2) - pow(y1, 2) - 2 * g * x1 - 2 * f * y1;
double h = -g;
double k = -f;
double r = sqrt(h * h + k * k - c);
print("Centre is at ($h, $k)");
print("\nCheck radius as the distance between the centre and the first point: $r");
}
void main() {
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31);
}
- Output:
Centre is at (18.978515660148815, 16.265410797715866) Check radius as the distance between the centre and the first point: 14.70862397833418
EasyLang
proc circ x1 y1 x2 y2 x3 y3 . cx cy cr .
x12 = x1 - x2
x13 = x1 - x3
y12 = y1 - y2
y13 = y1 - y3
y31 = y3 - y1
y21 = y2 - y1
x31 = x3 - x1
x21 = x2 - x1
sx13 = x1 * x1 - x3 * x3
sy13 = y1 * y1 - y3 * y3
sx21 = x2 * x2 - x1 * x1
sy21 = y2 * y2 - y1 * y1
f = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
g = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
c = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
cx = -g
cy = -f
cr = sqrt (cx * cx + cy * cy - c)
.
circ 22.83 2.07 14.39 30.24 33.65 17.31 cx cy cr
print "Centre: (" & cx & ", " & cy & ") Radius: " & cr- Output:
Centre: (18.98, 16.27) Radius: 14.71
F#
// Centre and radius of a circle passing through 3 points in a plane. Nigel Galloway: February 20th., 2024
let c (a,b) (c,d) (e,f)=(0.5*((a*a+b*b)*(f-d)+(c*c+d*d)*(b-f)+(e*e+f*f)*(d-b))/(a*(f-d)+c*(b-f)+e*(d-b)),
0.5*((a*a+b*b)*(e-c)+(c*c+d*d)*(a-e)+(e*e+f*f)*(c-a))/(b*(e-c)+d*(a-e)+f*(c-a)))
let d n g = let n,g=fst n-fst g,snd n-snd g in sqrt(n*n+g*g)
let circ P1 P2 P3 = let c=c P1 P2 P3 in (c,d c P1)
let centre,radius=circ (22.83, 2.07) (14.39, 30.24) (33.65, 17.31)
printfn $"Centre = %A{centre}, radius = %f{radius}"
- Output:
Centre = (18.97851566, 16.2654108), radius = 14.708624
FreeBASIC
Sub findCircle(x1 As Double, y1 As Double, x2 As Double, y2 As Double, x3 As Double, y3 As Double)
Dim x12 As Double = x1 - x2
Dim x13 As Double = x1 - x3
Dim y12 As Double = y1 - y2
Dim y13 As Double = y1 - y3
Dim y31 As Double = y3 - y1
Dim y21 As Double = y2 - y1
Dim x31 As Double = x3 - x1
Dim x21 As Double = x2 - x1
Dim sx13 As Double = x1 * x1 - x3 * x3
Dim sy13 As Double = y1 * y1 - y3 * y3
Dim sx21 As Double = x2 * x2 - x1 * x1
Dim sy21 As Double = y2 * y2 - y1 * y1
Dim f As Double = (sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13) / (y31 * x12 - y21 * x13) / 2
Dim g As Double = (sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13) / (x31 * y12 - x21 * y13) / 2
Dim c As Double = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
Dim h As Double = -g
Dim k As Double = -f
Dim r As Double = Sqr(h * h + k * k - c)
Print "Centre is at (" & h & ", " & k & ")"
Print "Radius is"; r
Print
Print "Check radius as the distance between the centre and the first point:"
Print Sqr((22.83 - h)^2 + (2.07 - k)^2)
End Sub
findCircle(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
Sleep
- Output:
Centre is at (18.97851566014882, 16.26541079771587) Radius is 14.70862397833418 Check radius as the distance between the centre and the first point: 14.70862397833418
jq
Adapted from Julia
Works with jq, the C implementation of jq
Works with gojq, the Go implementation of jq
# Emit {x,y,r} corresponding to the circle through the three points
# specified as [x,y] pairs.
def findcircle($p1; $p2; $p3):
def assertEq($p; $q): if ($p - $q)|length < 1e-12 then . else "assertion failed: \($p) != \($q)" | error end;
def ss($a; $b) : ($a|.*.) + ($b|.*.);
$p1 as [$a,$b]
| $p2 as [$c,$d]
| $p3 as [$e,$f]
| ($a - $e) as $ae
| ($d - $b) as $db
| ($b - $f) as $bf
| ($e - $c) as $ec
| ($c - $a) as $ca
| ($f - $d) as $fd
| ss($a; $b) as $a2b2
| ss($c; $d) as $c2d2
| ss($e; $f) as $e2f2
| {x: (0.5 * ($a2b2 * $fd + $c2d2 * $bf + $e2f2 * $db) / ($a * $fd + $c * $bf + $e * $db)),
y: (0.5 * ($a2b2 * $ec + $c2d2 * $ae + $e2f2 * $ca) / ($b * $ec + $d * $ae + $f * $ca)) }
# any one of these should do / be nearly identical:
| [ss(.x-$a; .y-$b), ss(.x-$c; .y-$d), ss(.x-$e; .y-$f)] as $r123
| assertEq( $r123|max; $r123|min )
| .r = (($r123 | add) / 3 | sqrt) ;
findcircle( [22.83, 2.07]; [14.39, 30.24]; [33.65, 17.31])
| "Centre = \([.x, .y]), radius = \(.r)"- Output:
Centre = [18.978515660148815,16.26541079771587], radius = 14.708623978334177
Julia
function findcircle(p1, p2, p3)
a, b = p1
c, d = p2
e, f = p3
a2b2 = a * a + b * b
ae = a - e
db = d - b
c2d2 = c * c + d * d
bf = b - f
ec = e - c
e2f2 = e * e + f * f
ca = c - a
fd = f - d
cx = 0.5 * (a2b2 * fd + c2d2 * bf + e2f2 * db) / (a * fd + c * bf + e * db)
cy = 0.5 * (a2b2 * ec + c2d2 * ae + e2f2 * ca) / (b * ec + d * ae + f * ca)
# any one of these should do / be nearly identical:
r123 = [(cx-a)^2 + (cy-b)^2, (cx-c)^2 + (cy-d)^2, (cx-e)^2 + (cy-f)^2]
@assert maximum(r123) - minimum(r123) < 1e-12
r = sqrt(sum(r123) / length(r123))
return (cx, cy), r
end
ctr, r = findcircle((22.83, 2.07), (14.39, 30.24), (33.65, 17.31))
println("Centre = $ctr, radius = $r")
- Output:
Centre = (18.978515660148815, 16.26541079771587), radius = 14.708623978334177
Kotlin
import kotlin.math.sqrt
data class Point(val x: Double, val y: Double)
fun findCircle(p1: Point, p2: Point, p3: Point): Pair<Point, Double> {
fun sq(x: Double) = x * x
val centreX =
0.5 * (
(sq(p1.x) + sq(p1.y)) * (p3.y - p2.y) +
(sq(p2.x) + sq(p2.y)) * (p1.y - p3.y) +
(sq(p3.x) + sq(p3.y)) * (p2.y - p1.y)
) / (
p1.x * (p3.y - p2.y) +
p2.x * (p1.y - p3.y) +
p3.x * (p2.y - p1.y)
)
val centreY =
0.5 * (
(sq(p1.x) + sq(p1.y)) * (p3.x - p2.x) +
(sq(p2.x) + sq(p2.y)) * (p1.x - p3.x) +
(sq(p3.x) + sq(p3.y)) * (p2.x - p1.x)
) / (
p1.y * (p3.x - p2.x) +
p2.y * (p1.x - p3.x) +
p3.y * (p2.x - p1.x)
)
val centre = Point(centreX, centreY)
val radius = sqrt(sq(centreX - p1.x) + sq(centreY - p1.y))
return Pair(centre, radius)
}
fun main() {
findCircle(Point(22.83,2.07), Point(14.39,30.24), Point(33.65,17.31))
.let { (c, r) ->
println("Centre = $c")
println("Radius = $r")
}
}
- Output:
Centre = Point(x=18.978515660148815, y=16.26541079771587) Radius = 14.70862397833418
Phix
with javascript_semantics
function circle(sequence p1, p2, p3)
atom {a,b} = p1, {c,d} = p2, {e,f} = p3,
a2b2 = a*a+b*b, ae = a-e, db = d-b,
c2d2 = c*c+d*d, bf = b-f, ec = e-c,
e2f2 = e*e+f*f, ca = c-a, fd = f-d,
cx = 0.5*(a2b2*fd+c2d2*bf+e2f2*db)/(a*fd+c*bf+e*db),
cy = 0.5*(a2b2*ec+c2d2*ae+e2f2*ca)/(b*ec+d*ae+f*ca)
-- any one of these should do / be nearly identical:
sequence r123 = {power(cx-a,2)+power(cy-b,2),
power(cx-c,2)+power(cy-d,2),
power(cx-e,2)+power(cy-f,2)}
assert((max(r123)-min(r123))<1e-12)
atom r = sqrt(average(r123))
return {{cx,cy},r}
end function
atom {{cx,cy},r} = circle({22.83,2.07},{14.39,30.24},{33.65,17.31})
printf(1,"Centre = %.8f, %.7f, radius = %.6f\n",{cx,cy,r})
- Output:
Centre = 18.97851566, 16.2654108, radius = 14.708624
Raku
Don't bother defining all the intermediate variables.
sub circle( (\𝒳ᵢ, \𝒴ᵢ), (\𝒳ⱼ, \𝒴ⱼ), (\𝒳ₖ, \𝒴ₖ) ) {
:center(
my $𝒞ₓ = ((𝒳ᵢ² + 𝒴ᵢ²) × (𝒴ₖ - 𝒴ⱼ) + (𝒳ⱼ² + 𝒴ⱼ²) × (𝒴ᵢ - 𝒴ₖ) + (𝒳ₖ² + 𝒴ₖ²) × (𝒴ⱼ - 𝒴ᵢ)) /
(𝒳ᵢ × (𝒴ₖ - 𝒴ⱼ) + 𝒳ⱼ × (𝒴ᵢ - 𝒴ₖ) + 𝒳ₖ × (𝒴ⱼ - 𝒴ᵢ)) / 2,
my $𝒞ᵧ = ((𝒳ᵢ² + 𝒴ᵢ²) × (𝒳ₖ - 𝒳ⱼ) + (𝒳ⱼ² + 𝒴ⱼ²) × (𝒳ᵢ - 𝒳ₖ) + (𝒳ₖ² + 𝒴ₖ²) × (𝒳ⱼ - 𝒳ᵢ)) /
(𝒴ᵢ × (𝒳ₖ - 𝒳ⱼ) + 𝒴ⱼ × (𝒳ᵢ - 𝒳ₖ) + 𝒴ₖ × (𝒳ⱼ - 𝒳ᵢ)) / 2
),
radius => (($𝒞ₓ - 𝒳ᵢ)² + ($𝒞ᵧ - 𝒴ᵢ)²).sqrt
}
say circle (22.83,2.07), (14.39,30.24), (33.65,17.31);
- Output:
(center => (18.97851566 16.2654108) radius => 14.70862397833418)
You may Attempt This Online!
Swift
import Foundation
import Matrix
extension Matrix where Element: SignedNumeric {
func minor(row: Int, column: Int) -> Element {
var submatrix = self
_ = submatrix.remove(rowAt: row - 1)
_ = submatrix.remove(columnAt: column - 1)
return submatrix.determinant as Element
}
}
enum MatrixErrors: Error {
case notEnoughPoints, tooManyPoints, pointsOnALine, miscError
}
func circleFrom3Points(points: (Double,Double)... ) throws -> (Double,Double,Double){
var pointArray: [[Double]] = [[0,0,0,0]]
for p in points {
pointArray.append([pow(p.0, 2) + pow(p.1, 2), p.0, p.1, 1])
}
guard pointArray.count > 3 else { throw MatrixErrors.notEnoughPoints }
guard pointArray.count < 5 else { throw MatrixErrors.tooManyPoints }
var matrix = Matrix(elements:pointArray)
var m11 = matrix.minor(row: 1, column: 1)
guard m11 != 0 else { throw MatrixErrors.pointsOnALine }
var m12 = matrix.minor(row: 1, column: 2)
var m13 = matrix.minor(row: 1, column: 3)
let x = 0.5 * m12 / m11
let y = -0.5 * m13 / m11
let r = (pow(x - pointArray[1][1],2) + pow(y - pointArray[1][2],2)).squareRoot()
return (x,y,r)
}
do {
let (x,y,r) = try circleFrom3Points(points: (22.83,2.07), (14.39,30.24), (33.65,17.31))
print("x:\(x), y:\(y), r: \(r)")
} catch {
debugPrint(error)
exit(1)
}
x:18.978515660148812, y:16.265410797715873, r: 14.708623978334185
Wren
This follows the lines of the C++ code here.
var findCircle = Fn.new { |x1, y1, x2, y2, x3, y3|
var x12 = x1 - x2
var x13 = x1 - x3
var y12 = y1 - y2
var y13 = y1 - y3
var y31 = y3 - y1
var y21 = y2 - y1
var x31 = x3 - x1
var x21 = x2 - x1
var sx13 = x1 * x1 - x3 * x3
var sy13 = y1 * y1 - y3 * y3
var sx21 = x2 * x2 - x1 * x1
var sy21 = y2 * y2 - y1 * y1
var f = sx13 * x12 + sy13 * x12 + sx21 * x13 + sy21 * x13
f = f / (y31 * x12 - y21 * x13) / 2
var g = sx13 * y12 + sy13 * y12 + sx21 * y13 + sy21 * y13
g = g / (x31 * y12 - x21 * y13) / 2
var c = -x1 * x1 - y1 * y1 - 2 * g * x1 - 2 * f * y1
var h = -g
var k = -f
var r = (h * h + k * k - c).sqrt
return [h, k, r]
}
var hkr = findCircle.call(22.83, 2.07, 14.39, 30.24, 33.65, 17.31)
System.print("Centre is at %([hkr[0], hkr[1]])")
System.print("Radius is %(hkr[2])")
System.print("\nCheck radius as the distance between the centre and the first point:")
System.print(((22.83 - hkr[0]).pow(2) + (2.07 - hkr[1]).pow(2)).sqrt)
- Output:
Centre is at [18.978515660149, 16.265410797716] Radius is 14.708623978334 Check radius as the distance between the centre and the first point: 14.708623978334