Centre and radius of a circle passing through 3 points in a plane: Difference between revisions

Line 441:
 
=={{header|jq}}==
'''Adapted from [[#WrenJulia|WrenJulia]]'''
 
'''Works with jq, the C implementation of jq'''
Line 449:
# Emit {x,y,r} corresponding to the circle through the three points
# specified as [x,y] pairs.
def findCirclefindcircle($p1; $p2; $p3):
 
def assertEq($p; $q): if ($p - $q)|length < 1e-12 then . else "assertion failed: \($p) != \($q)" | error end;
def s($a;$b): $a*$a + $b*$b;
def t($a;$b): $a*$a - $b*$b;
 
def ss($a; $p1b) as: [($x1,a|.*.) + ($y1]b|.*.);
| $p2 as [$x2, $y2]
| $p3p1 as [$x3a, $y3b]
| $p2 as [$x2c, $y2d]
| $p3 as [$e,$f]
 
| ($x1a - $x2e) as $x12ae
| ($x1d - $x3b) as $x13db
| ($y1b - $y2f) as $y12bf
| ($y1e - $y3c) as $y13ec
| ($y3c - $y1a) as $y31ca
| ($y2f - $y1d) as $y21fd
| ($x3 - $x1) as $x31
| ($x2 - $x1) as $x21
 
| tss($x1a; $x3b) as $tx13a2b2
| tss($y1c; $y3d) as $ty13c2d2
| tss($x2e; $x1f) as $tx21e2f2
| t($y2; $y1) as $ty21
 
| { yx: (0.5 * ($tx13a2b2 * $x12fd + $ty13c2d2 * $bf + $e2f2 * $db) / ($a * $x12fd + $tx21c * $x13bf + $ty21e * $x13db)),
xy: (0.5 * ($tx13a2b2 * $y12ec + $ty13c2d2 * $ae + $e2f2 * $ca) / ($b * $y12ec + $tx21d * $y13ae + $ty21f * $y13ca)) }
# any one of these should do / be nearly identical:
| .y = .y / ($y21 * $x13 - $y31 * $x12) / 2
| [ss(.x =-$a; .x /y-$b), ss(.x-$x21 *c; .y-$y13d), ss(.x- $x31 *e; .y-$y12f)] /as 2$r123
| assertEq( $r123|max; $r123|min )
| (2 * (.x * $x1 + .y * $y1) - s($x1; $y1)) as $c
| .r = (($r123 s(.x;| .yadd) -/ $c)3 | sqrt) ;
 
demofindcircle( [22.83, 2.07]; [14.39, 30.24]; [33.65, 17.31])
def demo($p1; $p2; $p3):
| "Centre is= at [\([.x, .y]), radius = \(.yr)]",
def sq: .*.;
def d($a; $b): (($a[0] - $b[0])|sq) + (($a[1] - $b[1])|sq) | sqrt;
findCircle($p1; $p2; $p3)
| "Centre is at [\(.x), \(.y)]",
"Radius is \(.r)",
 
"\nCheck radius as the distance between the centre and the first point: "
+ "\( d($p1; [.x,.y]) )"
;
 
demo([22.83, 2.07]; [14.39, 30.24]; [33.65, 17.31])
</syntaxhighlight>
{{output}}
<pre>
Centre is at= [18.978515660148815, 16.26541079771586626541079771587], radius = 14.708623978334177
Radius is 14.708623978334177
 
Check radius as the distance between the centre and the first point: 14.708623978334177
</pre>
 
2,458

edits