(42S) Determine Circle From Three Given Points

[gerry_in_polo]

Calculates the center and radius of the circle passing through three given points. Example, enter X1=1, Y1=4, X2=-1, Y2=2, X3=4, Y3=-3. Result, center (Xc,Yc)=(2.5,0.5) and radius R=3.8079.

[Dieter]

(07-15-2018 08:43 AM)gerry_in_polo Wrote:  Calculates the center and radius of the circle passing through three given points. Example, enter X1=1, Y1=4, X2=-1, Y2=2, X3=4, Y3=-3. Result, center (Xc,Yc)=(2.5,0.5) and radius R=3.8079.

Thank you very much for your program. This problem, solving a circle through three given points, reminds me of the early days with my first programmable calculator, a 34C. I had a HP67/97 book which included auch a program, and I tried to adapt it for the 34C.

Essentially the problem boils down to solving a linear equation system with three unknowns. In matrix notation it can be written this way:

Code:

(x1  y1  -1)     (2*xc          )     (x1² + y1²)
(x2  y2  -1)  *  (2*yc          )  =  (x2² + y2²)
(x3  y3  -1)     (xc² + yc² - r²)     (x3² + y3²)

Since the 42s features matrix support, this approach can be implemented in a quite compact program. So here is my alternative solution:

Code:

00 { 184-Byte Prgm }
01>LBL "CIRCLE"
02 3
03 ENTER
04 NEWMAT
05 STO "C"
06 3
07 ENTER
08 1
09 NEWMAT
10 STO "B"
11 1.003
12 STO 00
13>LBL 00
14 INDEX "C"
15 RCL 00
16 "X"
17 AIP
18 |-"^Y"
19 AIP
20 |-"?"
21 1
22 STOIJ
23 PROMPT
24 X<>Y
25 STOEL
26 J+
27 X^2
28 X<>Y
29 STOEL
30 J+
31 X^2
32 +
33 -1
34 STOEL
35 INDEX "B"
36 RCL 00
37 1
38 STOIJ
39 R^
40 STOEL
41 ISG 00
42 GTO 00
43 RCL "B"
44 RCL "C"
45 ÷
46 STO "X"
47 INDEX "X"
48 RCLEL
49 2
50 ÷
51 STO "XC"
52 X^2
53 I+
54 RCLEL
55 2
56 ÷
57 STO "YC"
58 X^2
59 +
60 I+
61 RCLEL
62 -
63 SQRT
64 STO "R"
65 CLV "B"
66 CLV "C"
67 CLV "X"
68 "Xc=\LF   "
69 ARCL "XC"
70 AVIEW
71 STOP
72 "Yc=\LF   "
73 ARCL "YC"
74 AVIEW
75 STOP
76 "R=\LF   "
77 ARCL "R"
78 AVIEW
79 END

Note: "|-" means "append", "^" is "↑" and "\LF" is a line feed.

Your example:

Code:

XEQ "CIRCLE"
X1↑Y^1?
           1  [ENTER]  4 [R/S]
X2↑Y^2?
          -1  [ENTER]  2 [R/S]
X3↑Y^3?
           4  [ENTER] -3 [R/S]
Xc=
   2,5
              [R/S]
Yc=
   0,5
              [R/S]
R=
   3,807786...

Since I am not very familiar with 42s matrix programming the program can probably be improved with more efficient code. So if someone can provide a better version: go ahead.

Dieter

[gerry_in_polo]

The matrix approach does not readily illustrate the basic principle of determining the circle which encompasses three points: the intersection of the two perpendicular inspectors of the two chords linking the first and second and then the second and third given points as the circle's center and subsequent computation of the radius.

[gerry_in_polo]

Another observation: your program is 79 lines long, mine is 87, not really a significant difference although a different approach is always worth knowing.

[Dieter]

(07-16-2018 01:19 AM)gerry_in_polo Wrote:  Another observation: your program is 79 lines long, mine is 87, not really a significant difference ...

Your program can even be shorter by using RCL-arithmetic.
On the other hand my version has a quite elaborate output routine. Line 65 ff. can be replaced with a simple

Code:

65 VIEW "XC"
66 STOP
67 VIEW "YC"
68 STOP
69 VIEW "R"

Then take a look at the byte count. With the above simplification it's down to 149 Bytes.

(07-16-2018 01:19 AM)gerry_in_polo Wrote:  ...although a different approach is always worth knowing.

Yes, that's exactly why I posted it. This is just another way of solving the same problem. Isn't this what forums are all about?
But let me explain why I chose this method:

The approach here is the general circle equation: (X–Xc)² + (Y–Yc)² = R². The program solves the resulting equation system for the three given points. There is one significant advantage of this method: it doesn't matter if two Y-values are the same. Try (0|2), (4|2) and (2|0) and the program correctly returns the center at (2|2) and a radius of 2. This will not work with the other approach where a division by (Y2–Y1) and (Y3–Y2) takes place.

Dieter

[Thomas Klemm]

The HP-42S allows to use complex numbers to solve the problem.

Formula
The three complex number \(a\), \(b\) and \(c\) represent the vertices of the triangle.
Define the following:
\[ \begin{eqnarray}
p=c+a \\
P=c-a \\
\\
q=a+b \\
Q=a-b
\end{eqnarray} \]
Then for the center \(z\) of the circle the following holds true:
(The \(\cdot\) is for the dot-product.)
\[ \begin{eqnarray}
P\cdot(z-\frac{p}{2})=0 \\
Q\cdot(z-\frac{q}{2})=0
\end{eqnarray} \]
This means that the vector from the center to the midpoint is perpendicular to the sides of the triangle.

Given:
\[ \begin{eqnarray}
u=P\cdot p \\
v=Q\cdot q
\end{eqnarray} \]

We can rewrite these equations:
\[ \begin{eqnarray}
P\cdot z=P\cdot\frac{p}{2}=\frac{u}{2} \\
Q\cdot z=Q\cdot\frac{p}{2}=\frac{v}{2}
\end{eqnarray} \]

Or then:
\[ \begin{eqnarray}
P\cdot 2z=u \\
Q\cdot 2z=v
\end{eqnarray} \]

Using:
\[ \begin{eqnarray}
P=P_x+iP_y \\
Q=Q_x+iP_y \\
\\
z=x+iy \\
w=v+iu
\end{eqnarray} \]

This equation can be written using matrix notation:
\[ \begin{bmatrix}
P_x & P_y \\
Q_x & Q_y
\end{bmatrix}

\begin{bmatrix}
2x \\
2y
\end{bmatrix}
=
\begin{bmatrix}
u \\
v
\end{bmatrix} \]

If we transpose this matrix we can rewrite the equation as:
\[ \begin{eqnarray}
g=Q_x+iP_x \\
h=Q_y+iP_y \\
\\
2x\cdot g+2y\cdot h=w
\end{eqnarray} \]

This equation can be cross multiplied both by \(h\) and \(g\):
(Reminder: \(h\times h=g\times g=0\).)
\[ \begin{eqnarray}
2x\cdot h\times g=h\times w \\
2y\cdot h\times g=w\times g=-g\times w
\end{eqnarray} \]

The determinant is:
\[
D=h\times g
\]

We end up with:
\[ \begin{eqnarray}
x=\frac{h\times w}{2D} \\
y=-\frac{g\times w}{2D}
\end{eqnarray} \]

Program

Code:

00 { 69-Byte Prgm }  ; X Y Z T
01 LBL "CENTER"      ; c b a
02 ENTER             ; c c b a
03 R↑                ; a c c b
04 STO+ ST Z         ; a c c+a=p b
05 -                 ; c-a=P p b b
06 LASTX             ; a P p b
07 R↓                ; P p b a
08 DOT               ; u b a a
09 STO 00            ; u b a a
10 X<> ST L          ; P b a a
11 R↓                ; b a a P
12 STO+ ST Z         ; b a a+b=q P
13 -                 ; a-b=Q q P P
14 DOT               ; v P P P
15 STO 01            ; v P P P
16 X<> ST L          ; Q P P P
17 COMPLEX           ; Qy Qx P P
18 R↑                ; P Qy Qx P
19 COMPLEX           ; Py Px Qy Qx
20 X<>Y              ; Px Py Qy Qx
21 R↓                ; Py Qy Qx Px
22 COMPLEX           ; h Qx Px Px
23 X<> ST Z          ; Px Qx h Px
24 COMPLEX           ; g h Px Px
25 RCL 01            ; v g h Px
26 RCL 00            ; u v g h
27 COMPLEX           ; w g h h
28 R↓                ; g h h w
29 CROSS             ; h×g=D h w w
30 STO÷ ST Z         ; D h w/D w
31 X<> ST L          ; g h w/D w
32 X<> ST Z          ; w/D g h w
33 CROSS             ; g×w/D=2x h w w
34 X<>Y              ; h 2x w w
35 LASTX             ; w/D h 2x w
36 CROSS             ; h×w/D=-2y 2x w w
37 +/-               ; 2y 2x w w
38 COMPLEX           ; 2z w w w
39 2                 ; 2 2z w w
40 ÷                 ; z w w w
41 END               ; z w w w

The radius of the circle isn't calculated. But that should be easy now that we know the center.

Example:
A = (1, 4)
B = (-1, 2)
C = (4, -3)

1 ENTER 4 COMPLEX
-1 ENTER 2 COMPLEX
4 ENTER -3 COMPLEX
XEQ "CENTER"
2.5 i0.5


Kind regards
Thomas

[Dieter]

(07-17-2018 01:21 PM)Thomas Klemm Wrote:  The HP-42S allows to use complex numbers to solve the problem.

Hey, that's a really great new approach! Hats off!

(07-17-2018 01:21 PM)Thomas Klemm Wrote:  The radius of the circle isn't calculated. But that should be easy now that we know the center.

The simplest way of calculating the radius may be storing one of the points (as a complex number) at the beginning (insert STO "P" between line 01 and 02), and add a final ENTER RCL–"P" ABS at the end. This will return the center coordinates in Y and the radius in X:


 1  ENTER  4  COMPLEX
-1  ENTER  2  COMPLEX
 4  ENTER -3  COMPLEX
XEQ "CENTER"
=>
2,5 i0,5
3,80788655...


Dieter