Intersection of circle chords

The following code contains 3 circles, the lines from each pair of circle's intersections (with circcirc fxn), the slope/y-intercepts for each of the 3 lines, and the intersections themselves (with plotted lines to confirm visually). I cannot figure out why the bjk intercept is approximately .158305 off from what it is supposed to be, and that is what I require help with. What could be causing this?

clc
clear all
% circle centers/radii
xCenterk = 5;
yCenterk = 2;
xCenteri = 5;
yCenteri = 5;
xCenterj = 3;
yCenterj = 3;
x=0:.001:100;
y=0:.001:100;
theta = 0 : 0.001 : 2*pi;
radiusk = 1.5;
radiusi = 2;
radiusj = 1.5;

%generate 3 circles with parameters above.
xk = radiusk * cos(theta) + xCenterk;
yk = radiusk * sin(theta) + yCenterk;
xi = radiusi * cos(theta) + xCenteri;
yi = radiusi * sin(theta) + yCenteri;
xj = radiusj * cos(theta) + xCenterj;
yj = radiusj * sin(theta) + yCenterj;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%plotting each circle pair's intersections (or chord)%%%%%%%%%%%%%%%%%%%%

[xoutij,youtij] = circcirc(xCenteri,yCenteri,radiusi,xCenterj,yCenterj,radiusj);
mapshow(xoutij,youtij,'DisplayType','point','Marker','o')
plot(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]))

[xoutik,youtik] = circcirc(xCenteri,yCenteri,radiusi,xCenterk,yCenterk,radiusk);
mapshow(xoutik,youtik,'DisplayType','point','Marker','o')
plot(line([xoutik(1,1) xoutik(1,2)],[youtik(1,1) youtik(1,2)]))

[xoutjk,youtjk] = circcirc(xCenterj,yCenterj,radiusj,xCenterk,yCenterk,radiusk);
mapshow(xoutjk,youtjk,'DisplayType','point','Marker','o')
plot(line([xoutjk(1,1) xoutjk(1,2)],[youtjk(1,1) youtjk(1,2)]))
hold on;

% each lines slope

slopeij = @(line) (youtij(1,2) - youtij(1,1))/(xoutij(1,2) - xoutij(1,1));
mij = slopeij(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]));

slopeik = @(line) (youtik(1,2) - youtik(1,1))/(xoutik(1,2) - xoutik(1,1));
mik = slopeik(line([xoutik(1,1) xoutik(1,2)],[youtik(1,1) youtik(1,2)]));

slopejk = @(line) (youtjk(1,2) - youtjk(1,1))/(xoutjk(1,2) - xoutjk(1,1));
mjk = slopejk(line([xoutjk(1) xoutjk(2)],[youtjk(1) youtjk(2)]));

% each lines y-intercept (round is used as a makeshift catch for slope = 0)

if round(100*mij)==0;
    bij = youtij(1,1);
else
    interceptij = @(line,m) xoutij(1,2) - m*xoutij(1,1);
    bij = interceptij(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]),mij);
end

if round(100*mik)==0;
    bik = youtik(1,1);
else
    interceptik = @(line,m) xoutik(1,2) - m*xoutik(1,1);
    bik = interceptik(line([xoutik(1,1) xoutik(1,2)],[youtik(1,1) youtik(1,2)]),mik);
end

if round(100*mjk)==0;
    bjk = youtjk(1,1);
else
    interceptjk = @(line,m) xoutjk(1,2) - m*xoutjk(1,1);
    bjk = interceptjk(line([xoutjk(1,1) xoutjk(1,2)],[youtjk(1,1) youtjk(1,2)]),mjk);
end

%%%Plotting intersection points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1intersect = (bik-bij)/(mij-mik);
y1intersect = mij*x1intersect + bij;
x2intersect = (bjk-bij)/(mij-mjk);
y2intersect = mij*x2intersect+bij;
x3intersect = (bjk-bik)/(mik-mjk);
y3intersect = mjk*x3intersect+bjk;
plot(x1intersect,y1intersect,'marker','o');
plot(x2intersect,y2intersect,'marker','o');
plot(x3intersect,y3intersect,'marker','o');
axis equal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

aij=0:.01:10;
b3=mij*aij+bij;
plot(aij,b3,'m');
axis equal

aik=0:.01:10;
b3=mik*aik+bik;
plot(aik,b3,'m');
axis equal

ajk=0:.01:10;
b3=mjk*ajk+bjk;
plot(ajk,b3,'m');
axis equal

plot(xk, yk,'b',xi,yi,'r',xj,yj,'c');
axis equal;
hold on;

[Sean]

3 line intersections

I just checked on Wolfram Alpha, and the 3 lines intersect properly, but not on MatLab--so I must be saying something wrong?

[Sean]

Something like

yij=(r1^2-x1^2-y1^2-r2^2+x2^2+y2^2-2*(x2-x1)*x)/(2*(y2-y1));
yik=(r1^2-x1^2-y1^2-r3^2+x3^2+y3^2-2*(x3-x1)*x)/(2*(y3-y1));
yjk=(r3^2-x3^2-y3^2-r2^2+x2^2+y2^2-2*(x2-x3)*x)/(2*(y2-y3));

where 1=i, 2=j, 3=k works instead of all of this code.

Answer 1

[Geoff Hayes]

Sean - rather than pasting 80+ lines of code in your question, it would be better to just attach it the m file to your question (use the paperclip button to do so).

A couple of comments on your code. Use of

 plot(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]))

is incorrect as using the line command creates the line on the figure, and the plot(line(...)) plots the handle returned by line at (1,h) (where h is the handle). So you can remove the plot command and just use

 line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]);

for all plotting of the chords.

It is unclear why an anonymous function is used to calculate the slope for each of these chords/lines

 slopeij = @(line) (youtij(1,2) - youtij(1,1))/(xoutij(1,2) - xoutij(1,1));
 mij = slopeij(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]));

Use of line as the input to the slopeij function should be avoided since line is a MATLAB built-in function (and the variable line is not used used at all!). For each of the slope calculations, you can just replace the above with

 mij = (youtij(1,2) - youtij(1,1))/(xoutij(1,2) - xoutij(1,1);

As for calculating the y-intercept, you may also want to check for when the slope is infinite (as well as just zero) to handle those cases where the line of intersection (between two circles) is vertical. Use

 isinf(mij)

to determine this.

Again, use of anonymous functions is not needed for the y-intercepts nor should you be using line as an input variable which is not even called. You can simply use the local variables that have been defined for the slope and the points on the line. So instead of

 interceptij = @(line,m) xoutij(1,2) - m*xoutij(1,1);
 bij = interceptij(line([xoutij(1,1) xoutij(1,2)],[youtij(1,1) youtij(1,2)]),mij);

with

 bij =  xoutij(1,2) - mij*xoutij(1,1);

But even this line of code seems to be incorrect. The equation of a line is y=mx*b where m is the slope and b is the y-intercept. Solving for b leads to

 b = y-m*x

which is different from the bij above which use x coordinates from two different points, rather than using the x and y coordinates from a single point. The above should then be replaced with

 bij =  youtij(1,1) - mij*xoutij(1,1);

This should occur for all three y-intercepts.

Try the above and see what happens!