OK, I understand that x/y swap. Let’s see what’s next in my tests.
Following Bourke’s solution, we had this to compute x and y of one of the desired points on the circle:
local p2_bourke_x = p2.x + (h*(p1.y-p0.y)/d)
local p2_bourke_y = p2.y + (h*(p1.x-p0.x)/d)
Note that the x calculation uses the y coordinates and the y calculation uses the x coordinates. Why is that, one wonders? Well, I figured it out, I’m not sure just how but ti goes like this:
If we were to interpolate between p0 and p1 in the obvious way, we would get a point along the p0-p1 line. But we want a point along the line perpendicular to that line. And given that a line’s slope is a/b, the slope of lines perpendicular to that one have slope -b/a. I was 90 percent of the way there and then searched to be sure and yes, the perpendicular line have slope -1/m if the line’s slope is m.
Slope is, as I learned it, “rise over run” or the y component over the x component. So the slope of the normal is “-run over rise”, and that’s why we swap x and y to move along the normal from P2, the crossing point.
There are, of course two points to be found, not just one. We’re not there yet.
Along the way I realized that my thinking was muddled, you can’t just add <a,h> to p0 to get the answer: that works only if the slope between the centers is zero. It would, I think, be possible to rotate that vector suitably and add it in but even then I’m not sure. So I’ll modify the current test to remove my scheme and then move on with Bourke’s.
Here again, we find that by deriving the math and going back to sources, I understand why the code works, which I did not and would not if I just tried to decode it from the source.
Here’s the test, revised:
_:test("calculate p3 per Bourke", 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
local p3_x = p2.x + (h*(p1.y-p0.y)/d) -- note x and y reversed
local p3_y = p2.y + (h*(p1.x-p0.x)/d)
local p3 = vector(p3_x, p3_y)
_:expect(p3.x).is(7.5)
_:expect(p3.y).is(6.61, 0.01)
end)
And it passes. Let’s do another step to get both values, with a separate function to do the job, since we need that sooner or later:
function circle_intersections(p0, r0, p1, r1)
local d = p1:dist(p0)
local a = calculate_a(p0, r0, p1, r1)
local h = math.sqrt(r0*r0 - a*a)
local dif = p1-p0
local p2 = p0 + dif * a / d
local step = vector(dif.y, dif.x) * h / d
return {p2+step, p2-step}
end
_:test("function returns both values", function()
local p0 = vector(0,0)
local r0 = 10
local p1 = vector(15,0)
local r1 = 10
local results = circle_intersections(p0, r0, p1, r1)
for i, result in results do
_:expect(result.x).is(7.5)
_:expect(math.abs(result.y)).is(6.61, 0.01)
end
end)
I just took the absolute value of the y coordinate in this case because it’s either 6.61 or -6.61. I did check that first.
So we are another small step toward what we really want.
We should probably deal with the special cases, which include
In fact, since this code works for us, let’s do what will probably be best for us: we intend only to use the two-point return. We’ll return an empty collection in the other cases.
_:test("function returns empty on errors", function()
local p0 = vector(0,0)
local r0 = 10
local p1 = vector(15,0)
local r1 = 5
local results = circle_intersections(p0, r0, p1, r1)
_:expect(#results).is(0)
end)
I added a quick check to the function:
function circle_intersections(p0, r0, p1, r1)
local d = p1:dist(p0)
if d >= r0 + r1 then
return {}
end
local a = calculate_a(p0, r0, p1, r1)
local h = math.sqrt(r0*r0 - a*a)
local dif = p1-p0
local p2 = p0 + dif * a / d
local step = vector(dif.y, dif.x) * h / d
return {p2+step, p2-step}
end
You might wonder why I didn’t raise an error if the conditions are not met. The answer is simple: I don’t want our code raising errors: you get crashed trains in that case. So I return a result that can be checked instead. We might decide on a different behavior … perhaps just returning a point along the line toward the given circle, reaching as far as we can. Now that I’ve thought of it, that might be better. But until this code gets some exercise, I can’t be sure.
You’ll also notice that I didn’t calculate the points x and y separately, as Bourke’s algorithm did. Instead, I computed the slope-inverted vector dif and used that. I think that is a righteous change that clarifies the code. I also think it needs additional testing because I didn’t negate one of the parameters when I created dif. Let’s do another example.
Also the next morning, let’s see what’s next. I think I’d like to check whether these new functions are actually returning results that meet the distance criteria. After all, that’s what we want, and it’s tedious at best to hand-calculate values.
Yes, perhaps I should have thought of this sooner. Yes, I’m happy to have thought of it now.
I’ll modify an existing test:
_:test("function at 45", function()
local p0 = vector(0,0)
local r0 = 10
local xy = math.sqrt(15*15/2)
local p1 = vector(xy,xy)
local r1 = 10
local results = circle_intersections(p0, r0, p1, r1)
verify_results(results, p0, r0, p1, r1)
end)
With that call in place, writing the function seems straightforward:
function verify_results(results, p0, r0, p1, r1)
_:expect(#results).is(2)
for i, result in results do
_:expect(result:dist(p0)).is(r0, 0.001)
_:expect(result:dist(p1)).is(r1, 0.001)
end
end
Disappointingly, the test fails. I think I’ll go back to Bourke’s formula for the results, since I’ve messed that up at least once already.
With his formula in place, I get this error:
Actual: 14.114378277661475, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 0.885621722338526, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 0.885621722338526, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 14.114378277661475, Expected: 10 +/- 0.001. Test: 'function at 45'.
With mine in place:
Actual: 14.114378277661475, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 0.885621722338526, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 0.885621722338526, Expected: 10 +/- 0.001. Test: 'function at 45'.
Actual: 14.114378277661475, Expected: 10 +/- 0.001. Test: 'function at 45'.
OK, whatever the defect is, it’s not in my formula - or it’s the same error in both. I’ll leave his in play while I try to see what has happened.
All I can think to do is to print intermediate results, which I’ll then try to verify graphically or via hand calculations.
I notice that the results coming out have equal x and y, and clearly they should not. I am reminded that the slope of a normal line is the negative of the reciprocal. I wonder if Bourke said +/- for the one and -/+ for the other and I didn’t notice.
Yes he did!
Fix that. Tests run. Whew!
Looking back at the code that the LLM gave me, it is in fact correct. My manual transcription was not. Some might argue that if I had just accepted the result, maybe checked it a bit, everything would have been fine. Even if I had just accepted it, everything would have been fine – at least I don’t see an error in the LLM code yet.
Well, if you want to take that chance, I can’t stop you. But I have an understanding of the code I have, and a sense of how I can modify and improve it, and I would not have that sense had I just copied the code.
And I know myself, and if I had decided to go with the LLM code I would not have tested it so extensively, because post-hoc testing is So Boring, and even if I did study and try to verify it, I wouldn’t understand it to the extent that I do now.
That said, I am pretty sure now that I have invested more time into doing this with code of my own than I would had I just gone with the LLM code. If speed without understanding is the game, this time the LLM won.
Me, I’m going to eat a bagel. Safe paths!