function [x_i,y_i] = line_intersect_b(A1,B1,C1,A2,B2,C2)
% function [x_i,y_i] = line_intersect_b(A1,B1,C1,A2,B2,C2)
% This function checks to see if two lines,
% both in the plane z = 0, intersect.
% If they do, it returns the coordinates of the
% point of intersection and plots the point.
% If the lines don't intersect, it returns values of
% NaN ("Not a Number") for x and y, and no plot
% is produced
% Input parameters
% A1, B1, C1 = coefficients of equation for line 1
% (A1)(x) + (B1)(y) = C1
% A2, B2, C2 = coefficients of equation for line 2
% (A2)(x) + (B2)(y) = C2
% The equations of the lines are
% A1.*x + B1.*y = C1
% A2.*x + B2.*y = C2
% or, in matrix form
% |A1 B1| |x| = |C1|
% |A2 B2| |y| = |C2|
% Example
% [x_i,y_i] = line_intersect_b(1,2,3,4,5,6)
% Name: Steve Martel
% Date: 9/15/04

 
% Check to see if lines both lines parallel the x- or y-axes.
if A1 == 0 & A2 == 0
    % Both lines parallel the x-axis and do not intersect
    x_i = NaN;
    y_i = NaN;
    return
elseif B1 == 0 & B2 == 0
    % Both lines parallel the y-axis and do not intersect
    x_i = NaN;
    y_i = NaN;
    return
end

% Check to see if either line parallels the x- or y-axes.
if A1 == 0 & B2 == 0
    % Line 1 parallels the x-axis and line 2 parallels the y-axis
    x_i = C2./A2;
    y_i = C1./B1;
    return
elseif B1 == 0 & A2 == 0 
    % Line 1 parallels the y-axis and line 2 parallels the x-axis
    x_i = C1./A1;
    y_i = C2./B2;
    return
end

% Check to see if just one line parallels the x- or y-axes.
if A1 == 0
    % Just line 1 parallels the x-axis
    y_i = C1./B1;
    x_i = (C1 - B1.*y_int)./A1;
    return
elseif B1 == 0
    % Just line 1 parallels the y-axis
    x_i = C1./A1;
    y_i = (C1 - B1.*x_int)./B1;
    return
elseif A2 == 0
    % Just line 2 parallels the x-axis
    y_i = C2./B2;
    x_i = (C2 - B2.*y_int)./A2;
    return
elseif B2 == 0
    % Just line 2 parallels the y-axis
    x_i = C2./A2;
    y_i = (C2 - B2.*x_int)./B2;
    return
end

% If the code gets to here, neither line parallels the x- axis or y-axis
% and the slopes aren't zero or infinity
% Find slopes of lines
m1 = A1./B1;
m2 = A2./B2;
% Check to see if the slopes are parallel
if m1 == m2
    % The lines are parallel to the y-axis and do not intersect
    x_i = NaN;
    y_i = NaN;
    return
end

% The lines are not parallel and therefore intersect
% Find the y-intercepts of each line
b1 = C1./B1;
b2 = C2./B2;

% Find the coordinates of the intersection
x_i = (b2 - b1)./(m1 - m2);
y_i = m1.*x_i + b1;

% Plot the intersection, extending the lines a bit from that point
% Prepare two subplots, the first showing the x- and y- intercepts
% of the lines, the second showing their intersection
x_int_1 = -b1./m1;
x_int_2 = -b2./m2;

figure(1)
clf
subplot(1,2,1)
plot([x_int_1 0],[0 b1],'b',[x_int_2 0],[0 b2],'r')
title (['Lines intersect at x = ',num2str(x_i),' y = ',num2str(y_i)])
xlabel('x')
ylabel('y')
    
subplot(1,2,2)
plot(x_i,y_i,'+')
title (['Lines intersect at x = ',num2str(x_i),' y = ',num2str(y_i)])
xlabel('x')
ylabel('y')
