(41) Intersection points between circles

[rawi]

Hi,

if somebody else already did a program to find the intersection points of two circles I could not find it.

So here is mine:

Program CIP
Program to determine the intersection points of two circles. Circles are defined by coordinates of central point and radius. Output is coordinates of the intercept points of the two circles.
If there is no solution message is given. If both solutions are the same there is only 1 solution.

Use of registers:
Circle 1: Coordinates of midpoint: x: 01, y: 02
Radius: 03
Circle 2: Coordinates of midpoint: x: 04, y: 05
Radius: 06
07: Distance between the two midpoints.
08: Angle between x-axis and line between midpoints
09: Angle at midpoint of first circle of the triangle formed by midpoints and first intersection of cycles
10: For computation
Output:
Coordinates of intersection points: Point 1: x: 11, y: 12; Point 2: x: 13, y: 14

Flags used: 01, 02

Example:
Circle 1 has the central point x=13, y=4 and the radius 6.8
Circle 2 has the central point x=6, y=10 and the radius 4

Keystrokes:
XEQ Alpha CIP Alpha
-> CIRCLE1? 13 ENTER 4 ENTER 6.8 R/S
-> CIRCLE2? 6 ENTER 10 ENTER 4 R/S
-> 9.9987 (x-value of 1st intersection point) X<>Y 10.1018 (y-value of 1st intersection point) – R/S -
-> 6.5109 (x-value of 2nd intersection point) X<>Y 6.0328 (y-value of 2nd intersection point)

A revised version that is shorter and does not use Flag 02 is below. Thanks to Albert Chan who showed me ways to improve the program!

Code:

01 LBL „CIP“
02 CF 01
03 CF 02
04 “CIRCLE1?”
05 PROMPT        Input of midpoint and radius of circle 1
06 STO 03
07 RDN
08 STO 02
09 RDN
10 STO 01
11 “CIRCLE2?”
12 PROMPT        Input of midpoint and radius of circle 2
13 STO 06
13 RDN
15 STO 05
16 RCL 02
17 X>Y?
18 SF 02
19 -
20 X²
21 X<>Y
22 STO 04
23 RCL 01 
24 X>Y?
25 SF 01
26 -
27 X²
28 +
29 SQRT
30 STO 07        Distance between midpoints
31 RCL 05
32 RCL 02
33 -
34 X<>Y
35 /
36 ASIN
37 ABS
38 STO 08        Angle between line between midpoints and x-axis
39 RCL 07
40 X²
41 RCL 03
42 X²
43 +
44 RCL 06
45 X²
46 -
47 2
48 /
49 RCL 07
50 /
51 RCL 03
52 /
53 ACOS
54 STO 09
55 +
56 COS
57 LASTX
58 SIN
59 RCL 03
60 *
61 FS? 02
62 CHS
63 RCL 02
64 +
65 STO 12
66 X<>Y
67 RCL 03
68 *
69 FS? 01
70 CHS
71 RCL 01 
72 +
73 STO 11
74 STOP        X-register: x-coordinate, Y register: y-coordinate of first intersection point
75 RCL 08
76 RCL 09
77 -
78 COS
79 LASTX
80 SIN
81 RCL 03
82 *
83 FS? 02
84 CHS
85 RCL 02
86 +
87 STO 14
88 X<>Y
89 RCL 03
90 *
91 FS? 01
92 CHS
93 RCL 01
94 +
95 STO 13
96 GTO 02
97 LBL 01
98 “NO SOLUTION”
99 AON
100 STOP
101 AOFF
102 LBL 02
103 CF 01
104 CF 02
105 END

[Albert Chan]

Hi, rawi

We can get the 2 circles intersection points, without trig. functions.

lua> I = require'complex'.I
lua> p1, r1 = 13+4*I, 6.8
lua> p2, r2 = 6+10*I, 4.0
lua> d = (p2-p1):abs()
lua> d -- center-to-center distance
9.219544457292887

|r1-r2| ≤ d ≤ |r1+r2|, thus 2 circles do intersect somewhere.

My code, from Area of Triangle post

Code:

function area(a,b,c)                -- area of triangle: side a,b,c
    if a > b then a, b = b, a end
    if a > c then a, c = c, a end   -- a = shortest side, b-c guaranteed exact
    b, c = b-c, b*c
    b = (a+b)*(a-b)/4
    return sqrt(b*(c-b)), b/c       -- area, sin(A/2)^2
end

lua> h = area(r1,r2,d) * 2 / d -- height of triangle
lua> d1 = sqrt((r1+h)*(r1-h)) -- base of triangle
lua> d1, h
6.249766489755484      2.6796303520316775

If p1=(0,0), p2=(d,0), then intersected points at (d1, ±h)
Note: we setup with r1 ≥ r2, to force d1 ≥ d/2 > 0, see post #4

-- Complex multiply and addition for rotatation and translation, to get intersection points:

lua> v1 = (p2-p1) / d
lua> p1 + v1*(d1+h*I)
6.51094321225877+6.032767080968563*I
lua> p1 + v1*(d1-h*I)
9.99870384656476+10.101821154325554*I

---

Update: we can remove r1 ≥ r2 requirement, by combine this post and post #4
It might also be faster, since real part of v2 does not use square root.

Note: Lua math.sqrt(x) return -NaN if x < 0, is exactly what we wanted.
With "bad" triangle, area = -NaN ⇒ v2 = -NaN-NaN*I ⇒ 2 circles do not intersect.

lua> v2 = (1/2 + (r1+r2)*(r1-r2)/(2*d*d)) + 2*area(r1,r2,d) / (d*d) * I
lua> p1 + (p2-p1)*v2
6.510943212258768+6.032767080968563*I
lua> p1 + (p2-p1)*v2:conj()
9.99870384656476+10.101821154325554*I

Intersection points might be closer to mid-point of p1, p2. We may use that as refercence.

Code:

function intersection(p1, r1, p2, r2)
    local d = (p2-p1):abs()
    local v2 = ((r1+r2)*(r1-r2) + 4*area(r1,r2,d)*I) / (d*d)
    p1, p2 = (p1+p2)/2, (p2-p1)/2   -- mid-point as reference    
    return p1 + p2*v2, p1 + p2*v2:conj()
end

lua> intersection(p1, r1, p2, r2)
6.510943212258769+6.032767080968563*I      9.99870384656476+10.101821154325554*I

[rawi]

Hi Albert,

thank you very much. As a statistician up to now I had not to handle complex numbers. So - frankly spoken - I have some difficulties to understand your solution. But more and more I have the feeling that I should dig in complex numbers. There seem to be a lot of things that can be solved more elegant and easier than by solving it only with real numbers.

[Albert Chan]

Another way is to solve for x, then get h (half chord length)

XCAS> circle1 := x^2 + y^2 - r1^2 = 0
XCAS> circle2 := (x-d)^2 + y^2 - r2^2 = 0
XCAS> chord := simplify(circle1 - circle2); // equation for chord

\((2d)x + (r_2^2 - r_1^2 - d^2)=0\\⇒ x = \large{d\over2} + {r_1^2 - r_2^2 \over 2d}\)

XCAS> p1, r1 := 13 + 4*i, 6.8
XCAS> p2, r2 := 6 + 10*i, 4.0
XCAS> d := abs(p1-p2)                                 → √85 ≈ 9.21954445729
XCAS> d1 := d/2 + (r1+r2)*(r1-r2)/(2*d)      → 6.24976648976
XCAS> h := sqrt((r1+d1)*(r1-d1))                 → 2.67963035203

Rotate/translate d1±h*i, back to intersection coordinates.

XCAS> v1 := (p2-p1)/d;   → (-7+6*i)/√85 = sign(p2-p1)
XCAS> p1 + v1*(d1+h*i)   → 6.51094321226 + 6.03276708097*i
XCAS> p1 + v1*(d1-h*i)   → 9.99870384656 + 10.1018211543*i

[Albert Chan]

We can also get x directly, from Law of Cosine.
Let R₂ = angle corresponded to side r₂ (the blue side)

\(r_2^2 = d^2 + r_1^2 - 2\;d\;r_1 \cos(R_2) = d^2 + r_1^2 - 2\;d\;x
\\ ⇒ x = \large{d\over2} + {r_1^2 - r_2^2 \over 2d}\)

Let Δ = area of the triangle, vertices = (0,0), (x,h), (d,0)

\(Δ = \large{d\;h\over2}\)
\(⇒ h = \large{2Δ \over d}\)

[rawi]

Hi Albert,

thank you very much for the interesting discussion.
The solution that is used in my program indeed uses as well the Law of Cosine. The difference between your solution and mine is that I use directly the coordinates given. You first transfer it to an easier problem by setting one midpoint to (0;0) and the other to (d;0), but then you have it to transfer back. I am not yet sure which solution is easier.

Unfortunately I do not know how to insert a image that is not on a website. So I try to describe it without a picture.

Let
(x1,y1) be the coordinates of the first midpoint
(x2,y2) be the coordinates of the second midpoint
r1 be the radius of the first circle
r2 be the radius of the second circle

First the distance between midpoints is computed (steps 16-30): z=((x1-x2)^2+(y1-y2)^2)^0.5

Then we compute the angle between line between midpoints and paralell to x-axis (steps31-38):
alpha = arcsin((y2-y1)/z)

The triangle made by the midpoints and the first intersection have the sides z, r1 and r2. Therefore we can compute the angle between line between midpoints and line from midpoint of first circle to first intersection point by using the Law of Cosine (steps 39-54):

beta = arccos((z²+r1²-r2²)/(2*z*r1))

Now we must take into account, whether x2<x1 or x2>=x1 and analogue with y1 and y2.
We set sg1 = 1 if x2>=x1 and sg1 = -1 if x2<x1
and likewise sg2 = 1 fi y2 >= y1 and sg2 = -1 if y2 < y1.
This is done by using flag 01 and flag 02.
Then we can compute the coordinates of the two intersection points (sx1, sy1) (steps 55-73) and (sx2, sy2) (steps 75-95):

sx1 = x1 + r1*(cos(alpha+beta))*sg1
sy1 = y1 + r1*(sin(alpha+beta))*sg2
sx2 = x1 + r1*(cos(alpha-beta))*sg1
sy2 = y1 + r1*(sin(alpha-beta))*sg2

This is the solution given.
I hope this makes it a little more transparent.

Best

Raimund

[SlideRule]

perhaps Computer Numerical Control, pages 37 & 38 helps?

BEST!
SlideRule

[Albert Chan]

(11-30-2020 05:19 PM)rawi Wrote:  Then we compute the angle between line between midpoints and paralell to x-axis (steps31-38):
alpha = arcsin((y2-y1)/z)

What we wanted is the angle between vector ((x1,y1) to (x2,y2)), and the x-axis.
We should have: -pi ≤ alpha ≤ +pi

Using arcsin, we have -pi/2 ≤ alpha ≤ pi/2. sg1/sg2 correction will not fix this.
Correction should be:

      alpha = arcsin((y2-y1)/z); if x2 < x1 then alpha = pi - alpha end

Or, alpha = arccos((x2-x1)/z); if y2 < y1 then alpha = -alpha end

Code:

function intersection(p1, r1, p2, r2)
    p2 = p2 - p1                    -- p2 relative to p1
    local d, a = p2:abs(), p2:arg() -- polar form
    local b = acos((d*d+(r1+r2)*(r1-r2)) / (2*d*r1))
    return p1 + r1*((a+b)*I):exp(), p1 + r1*((a-b)*I):exp()
end

lua> I = require'complex'.I
lua> p1, r1 = 13+4*I, 6.8
lua> p2, r2 = 6+10*I, 4.0
lua>
lua> intersection(p1,r1, p2,r2)
6.510943212258769+6.032767080968563*I      9.998703846564759+10.101821154325552*I
lua> intersection(p2,r2, p1,r1)
9.99870384656476+10.101821154325554*I      6.51094321225877+6.032767080968563*I

Quote:The difference between your solution and mine is that I use directly the coordinates given. You first transfer it to an easier problem by setting one midpoint to (0;0) and the other to (d;0), but then you have it to transfer back. I am not yet sure which solution is easier.

Your code do the polar form side, mine the rectangular side. Both ways are equivalent.

circles intersections = p1 + r1 * exp((a±b)*I) = p1 + exp(a*I) * (r1*exp(±b*I))

exp(a*I) = (p2-p1) / d
r1*exp(±b*I) = d1 ± h*I

[Albert Chan]

(11-30-2020 05:19 PM)rawi Wrote:  alpha = arcsin((y2-y1)/z)
...
We set sg1 = 1 if x2>=x1 and sg1 = -1 if x2<x1
and likewise sg2 = 1 fi y2 >= y1 and sg2 = -1 if y2 < y1.

sx1 = x1 + r1*(cos(alpha+beta))*sg1
sy1 = y1 + r1*(sin(alpha+beta))*sg2
sx2 = x1 + r1*(cos(alpha-beta))*sg1
sy2 = y1 + r1*(sin(alpha-beta))*sg2

If x2 ≥ x1, we have the correct angle. No correction is needed.

If x2 < x1, calculated alpha should really be pi - alpha

cos((pi - alpha) ± beta) = cos(pi - (alpha ∓ beta)) = cos(alpha ∓ beta) * -1
sin((pi - alpha) ± beta) = sin(pi - (alpha ∓ beta)) = sin(alpha ∓ beta)

→ sy does not require correction.

Code:

function intersection(p1,r1, p2, r2)
    p2 = p2 - p1        -- p2 relative to p1
    local d = p2:abs()
    local a = asin(p2:imag()/d)
    local b = acos((d*d+(r1+r2)*(r1-r2)) / (2*d*r1))
    local sg1 = p2:real()<0 and -1 or 1
    return p1 + r1*(cos(a+b)*sg1 + sin(a+b)*I),
           p1 + r1*(cos(a-b)*sg1 + sin(a-b)*I)
end

lua> I = require'complex'.I
lua> p1, r1 = 13+4*I, 6.8
lua> p2, r2 = 6+10*I, 4.0
lua>
lua> intersection(p1,r1,p2,r2)
9.99870384656476+10.101821154325554*I      6.510943212258768+6.032767080968563*I
lua> intersection(p2,r2,p1,r1)
9.99870384656476+10.101821154325554*I      6.51094321225877+6.032767080968563*I

[rawi]

Hi Albert,

You wrote:

Quote:We should have: -pi ≤ alpha ≤ +pi

Using arcsin, we have -pi/2 ≤ alpha ≤ pi/2. sg1/sg2 correction will not fix this.
Correction should be:

      alpha = arcsin((y2-y1)/z); if x2 < x1 then alpha = pi - alpha end

Or, alpha = arccos((x2-x1)/z); if y2 < y1 then alpha = -alpha end

I forgot to mention that I take the absolute figure for computing alpha. I compute
alpha = abs(arcsin((y2-y1)/z)) (see line 37 of program listing). So I limit it to the interval from 0 to pi/2.

You wrote:

Quote:sy does not require correction.

The results tell me that the correction is needed.
Take the following example:
x1 = -1, y1= -3, r1 = 4; x2 = 6, y2 = 2, r2 = 5
If computed it with correction I get:
sx1 = 1,3106; sy1=0.2651; sx2= 2,8380; sy2 = -1.8732 which is correct
If I keep the circles, but use the former circle 2 as circle 1 and vice versa I get:
sx1 = 2,8380; sy1= -1.8732; sx2 = 1.3106; sy2 = 0.2651
So the numbers of the intersection points are exchanged, but they still are correct.

If I compute without correction of sy (I delete lines 83, 84, 61, 62 of the program) I get:
For x1=-1, y1=-3, r1 = 4, x2 = 6, y2 = 2, r2 = 5 the same results as before.
But if I change the number of circles, i.e.
x1=6, y1=2, r1=5, x2=-1, y2=-3, r2=4
I get
x1 = 2.8380, y1=5.8380, x2=1,3106, y2=3,7349 which is not correct for the y-coordinates.

Perhaps by limiting alpha from 0 to pi/2 and doing the correction both for sx and sy I get correct results? I tested the program with numerous examples and it always delivered correct results.

Best
Raimund

[Albert Chan]

(12-01-2020 09:25 AM)rawi Wrote:  I forgot to mention that I take the absolute figure for computing alpha. I compute
alpha = abs(arcsin((y2-y1)/z)) (see line 37 of program listing). So I limit it to the interval from 0 to pi/2.
...
If I compute without correction of sy (I delete lines 83, 84, 61, 62 of the program) I get:
For x1=-1, y1=-3, r1 = 4, x2 = 6, y2 = 2, r2 = 5 the same results as before.
But if I change the number of circles, i.e.
x1=6, y1=2, r1=5, x2=-1, y2=-3, r2=4
I get
x1 = 2.8380, y1=5.8380, x2=1,3106, y2=3,7349 which is not correct for the y-coordinates.

Try it again, but remove the ABS applied to alpha.
It is this extra ABS that made the program complicated, requiring sg2 correction to undo.

Quote:Perhaps by limiting alpha from 0 to pi/2 and doing the correction both for sx and sy I get correct results?

Yes, it does. Here is the proof that sg2 correction undo ABS.

If y2-y1 ≥ 0, ABS does nothing to alpha, and we had already shown it work in previous post.

If y2-y1 < 0, we have 2 possibilities:

x2-x1 ≥ 0: calculated alpha should really be -alpha:
cos(-alpha ± beta) = cos(alpha pi ∓ beta) * 1
sin(-alpha ± beta) = sin(alpha pi ∓ beta) * -1

x2-x1 < 0: calculated alpha should really be -(pi - alpha):
cos(-(pi - alpha) ± beta) = cos(pi - (alpha ± beta)) = cos(alpha ± beta) * -1
sin(-(pi - alpha) ± beta) = -sin(pi - (alpha ± beta)) = sin(alpha ± beta) * -1

[rawi]

Hi Albert,

I tried it and you are totally right. Thanks for pointing that out to me.

Here is the improved method description:

Let (x1, y1) be the coordinates of midpoint of circle 1 and (x2, y2) be the coordinates of midpoint of circle 2 and r1 resp. r2 be the radii.
Let z be the distance between the midpoints: z=((x1-x2)²+(y1-y2)²)^0.5
The angle alpha between the line between the midpoints and the parallel to the x-axis through (x1, y1) using the sinus-function is:

alpha = arcsin((y2-y1)/z)
Beta can be computed by using the Law of Cosine:

beta = arccos((z²+r1²-r2²)/(2*z*r1))
We must take into account, whether x2<x1 or x2>=x1.
We set sg = 1 if x2>=x1 and sg1 = -1 if x2<x1
Then the coordinates of the two intersection points are:

sx1 = x1 + r1*(cos(alpha+beta))*sg
sy1 = y1 + r1*(sin(alpha+beta))
sx2 = x1 + r1*(cos(alpha-beta))*sg
sy2 = y1 + r1*(sin(alpha-beta))

And here is the improved code, which is shorter and easier than the previous one:

Code:

01 LBL „CIP“
02 CF 01
03 “CIRCLE1?”
04 PROMPT        Input of midpoint and radius of circle 1
05 STO 03
06 RDN
07 STO 02
08 RDN
09 STO 01
10 “CIRCLE2?”
11 PROMPT        Input of midpoint and radius of circle 2
12 STO 06
13 RDN
14 STO 05
15 RCL 02
16 -
17 X²
18 X<>Y
19 STO 04
20 RCL 01 
21 X>Y?
22 SF 01
23 -
24 X²
25 +
26 SQRT
27 STO 07        Distance between midpoints
28 RCL 05
29 RCL 02
30 -
31 X<>Y
32 /
33 ASIN
34 STO 08        Angle between line between midpoints and x-axis
35 RCL 07
36 X²
37 RCL 03
38 X²
39 +
40 RCL 06
41 X²
42 -
43 2
44 /
45 RCL 07
46 /
47 RCL 03
48 /
49 ACOS
55 STO 09
51 +
52 COS
53 LASTX
54 SIN
55 RCL 03
56 *
57 RCL 02
58 +
59 STO 12
60 X<>Y
61 RCL 03
62 *
63 FS? 01
64 CHS
65 RCL 01 
66 +
67 STO 11
68 STOP        X-register: x-coordinate, Y register: y-coordinate of first intersection point
69 RCL 08
70 RCL 09
71 -
72 COS
73 LASTX
74 SIN
75 RCL 03
76 *
77 RCL 02
78 +
79 STO 14
80 X<>Y
81 RCL 03
82 *
83 FS? 01
84 CHS
85 RCL 01
86 +
87 STO 13
88 GTO 02
89 LBL 01
90 “NO SOLUTION”
91 AON
92 STOP
93 AOFF
94 LBL 02
95 CF 01
96 END

[Albert Chan]

(12-01-2020 02:04 PM)rawi Wrote:  sx1 = x1 + r1*(cos(alpha+beta))*sg
sy1 = y1 + r1*(sin(alpha+beta))
sx2 = x1 + r1*(cos(alpha-beta))*sg
sy2 = y1 + r1*(sin(alpha-beta))

To reduce errors, it may be better if we pick circle1 with the smaller radius, r1 ≤ r2.

With circle1 having smaller radius, beta could be huge. (up to pi)
We could use complement of beta instead:

gamma = arcsin((z²+r1²-r2²)/(2*z*r1)) = pi/2 - beta

With gamma, the coordinates are:

sx1 = x1 + r1 * sin(gamma-alpha) * sg
sy1 = y1 + r1 * cos(gamma-alpha)
sx2 = x1 + r1 * sin(gamma+alpha) * sg
sy2 = y1 − r1 * cos(gamma+alpha)

Prove:

If x ≥ 0, sg = 1
cos(alpha ± beta) = cos(beta ± alpha) = sin(gamma ∓ alpha) * 1
sin(alpha ± beta) = ±sin(beta ± alpha) = ±cos(gamma ∓ alpha)

If x <0, sg = -1, alpha should really be pi - alpha
cos(pi-alpha ∓ beta) = cos(alpha ± beta) * -1
sin(pi-alpha ∓ beta) = sin(alpha ± beta)