Now that I’ve done the math, let’s code the function. With tests, of course.
The key issue in finding the points on a given circle that are a given distance from a given point (they sure do give us a lot of stuff, don’t they), the key issue is finding the point on the line between the given point and the circle center where the normal to that line will intersect the circle at the desired points. Then since you have the lengths of the hypotenuses (hypoteneece?), given the lengths along the line, which is the base of a right triangle, you can get the height by the Pythagorean Theorem which we all know and love.
We’ll see all that unfold shortly, I hope. Yesterday, after briefly falling into the LLM trap, I derived the formula myself, as shown in yesterday’s article. Now let’s build some code, supported by tests.
In yesterday’s work, calling the distance from the given circle to the key point a, the formula we derived, with plenty of help from Paul Bourne and some kind person on Stack Overflow was:
a = (d^2 + r0^2 - r1^2)/ 2*d
where d = the distance between the two circle centers.
I’ll write that function and then figure out how to test it.
function _:featureCircleCircle()
_:describe("TDD Paul Bourke circle intersections", function()
function calculate_a(p0, r0, p1, r1)
local d = p0:dist(p1)
return (d*d + r0*r0 - r1*r1) / (2*d)
end
_:test("trivial case", function()
local p0 = vector(0,0)
local r0 = 1
local p1 = vector(10,0)
local r1 = 9
local a = calculate_a(p0, r0, p1, r1)
_:expect(a).is(1, 0.0001)
end)
end)
end
The two circles just touch, so the a value is 1. Like the test says, trivial case. Let’s do a harder one, but still easy. How about two circles, radius 1, spaced apart by 1. THe a value should be 1/2, I think.
_:test("close and equal radii case", function()
local p0 = vector(0,0)
local r0 = 1
local p1 = vector(1,0)
local r1 = 1
local a = calculate_a(p0, r0, p1, r1)
_:expect(a).is(0.5, 0.0001)
end)
I could do more extensive testing, and perhaps it would be wise to do so, but I want to move on to the actual answer.
We have a right triangle with base a and hypotenuse r0, so the height of the triangle can be calculated and tested like this:
_:test("calculate h", function()
local p0 = vector(0,0)
local r0 = 10
local p1 = vector(15,0)
local r1 = 10
local a = calculate_a(p0, r0, p1, r1)
_:expect(a).is(7.5, 0.001)
local h = math.sqrt(r0*r0-a*a)
_:expect(h).is(6.61, 0.01)
end)
That’s the example from yesterday. I had calculated the 6.61 by hand, as shown in one of yesterdays’ pictures.
Now we have that value, and we need to convert it to a point. The points of intersection, it seems to me, will be at the given circle’s center plus <a, h> or plus <a, -h>
Bourne does it another way. He calculates the point along the center-connecting line via lerp:
p2 = p0 + a*(p1-p0)/d
And then the two answer points as
p31 = p2 + h*(p1-p0)/d
p32 = p2 - h*(p1-p0)/d
Let’s test it both ways. First time I run the test, I think my answer is right and his wrong. In any case they do not agree. Perhaps I transcribed his solution incorrectly. I’ll remove that code and do again.
Same result. Here’s the test and results:
_:test("calculate p3 both ways", function()
local p0 = vector(0,0)
local r0 = 10
local p1 = vector(0,15)
local r1 = 10
local a = calculate_a(p0, r0, p1, r1)
local h = math.sqrt(r0*r0-a*a)
local d = p1:dist(p0)
local p2 = p0 + (p1-p0)*a / d
local p2_bourke = p2 + (p1-p0)*h / d
local p2_janet = p0 + vector(a, h)
print("p2 bourke", p2_bourke)
print("p2 janet", p2_janet)
_:expect(p2_janet).is(p2_bourke)
end)
Feature: TDD Paul Bourke circle intersections
p2 bourke <0, 14.114378277661476, 0>
p2 janet <7.5, 6.614378277661476, 0>
Actual: table[y]=6.614378277661476, Expected: 14.114378277661476. Test: 'calculate p3 both ways'.
Actual: table[x]=7.5, Expected: 0. Test: 'calculate p3 both ways'.
Same thing again. I think I’m right. So does the code stolen from the LLM.
Ah. Transcribed wrong He uses y in calculating x and vice versa. Redo test:
_:test("calculate p3 both ways", function()
local p0 = vector(0,0)
local r0 = 10
local p1 = vector(15,0)
local r1 = 10
local a = calculate_a(p0, r0, p1, r1)
local h = math.sqrt(r0*r0-a*a)
local d = p1:dist(p0)
local p2 = p0 + (p1-p0)*a / d
print("p2", p2)
local p2_bourke_x = p2.x + (h*(p1.y-p0.y)/d)
local p2_bourke_y = p2.y + (h*(p1.x-p0.x)/d)
local p2_bourke = vector(p2_bourke_x, p2_bourke_y)
local p2_janet = p0 + vector(a, h)
print("p2 bourke", p2_bourke)
print("p2 janet", p2_janet)
_:expect(p2_janet).is(p2_bourke)
end)
Now he agrees with me. So far so good. Good progress. I want to check my shorter scheme for something not on the x axis. For now, though it’s time for a break and perhaps a bit of Sunday breakfast.
Safe paths!