(41) Intersection points between circles - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP-41C Software Library (/forum-11.html) +--- Thread: (41) Intersection points between circles (/thread-15954.html) (41) Intersection points between circles - rawi - 11-25-2020 03:07 PM 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 RE: Intersection points between circles - Albert Chan - 11-25-2020 06:16 PM 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 RE: Intersection points between circles - rawi - 11-25-2020 06:48 PM 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. RE: Intersection points between circles - Albert Chan - 11-26-2020 12:22 AM 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 RE: Intersection points between circles - Albert Chan - 11-30-2020 02:19 PM We can also get x directly, from Law of Cosine. Let R2 = angle corresponded to side r2 (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}$$ RE: Intersection points between circles - rawi - 11-30-2020 05:19 PM 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 and analogue with y1 and y2. We set sg1 = 1 if x2>=x1 and sg1 = -1 if x2= 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 RE: Intersection points between circles - SlideRule - 11-30-2020 08:51 PM perhaps Computer Numerical Control, pages 37 & 38 helps? BEST! SlideRule RE: Intersection points between circles - Albert Chan - 11-30-2020 09:03 PM (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 RE: Intersection points between circles - Albert Chan - 12-01-2020 12:17 AM (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= 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 RE: Intersection points between circles - rawi - 12-01-2020 09:25 AM 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 RE: Intersection points between circles - Albert Chan - 12-01-2020 11:59 AM (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 RE: Intersection points between circles - rawi - 12-01-2020 02:04 PM 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. We set sg = 1 if x2>=x1 and sg1 = -1 if x2Y 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` RE: Intersection points between circles - Albert Chan - 12-10-2020 01:34 PM (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)