function [d,u,v] = InsideDist(x,y,theta)
% x and y are column vectors that define a polygon.
% theta is the angle associated with the ray
% If the ray does not intersect the polygon then d = 0 and
%    u and v are empty vectors.
% If the ray intersects the polygon then u and v are column
%      vectors with the property that (u(i),v(i)) are the 
%      intersection points arranged with increasing distance
%      from the origin.
n = length(x);
% Set up A and B arrays...
A = zeros(n,4);
B = [x y];
idx = [(2:n) 1];
A = [cos(theta)*ones(n,1) x-x(idx) sin(theta)*ones(n,1) y-y(idx)];
% Solve the 2x2 systems...
% Pivot if necessary. (Note that all the systems have the
% same first column)
if abs(cos(theta))<abs(sin(theta))
    % Swap the 2 equations...
    A(:,[1 3]) = A(:,[3 1]);
    A(:,[2 4]) = A(:,[4 2]);
    B(:,[1 2]) = B(:,[2 1]);
end
L21 = A(:,3)./A(:,1);
U22 = A(:,4) - L21.*A(:,2);
X = B;
X(:,2) = X(:,2) - L21.*X(:,1);
X(:,2) = X(:,2)./U22;
X(:,1) = (X(:,1) - A(:,2).*X(:,2))./A(:,1);
% Find all the "valid" intersection points...
idx = find(0<=X(:,1) & 0<=X(:,2) & X(:,2)<=1 );
u = X(idx,1)*cos(theta);
v = X(idx,1)*sin(theta);
% Compute the inside distance...
if length(u)>0
    dist = u.^2 + v.^2;
    [z,jdx] = sort(dist);
    u = u(jdx);
    v = v(jdx);
    M = length(u);
    d = sum(sqrt((u(1:2:M)-u(2:2:M)).^2 + (v(1:2:M)-v(2:2:M)).^2));
else
    d = 0;
end

% 1 point for setting up A and B 
% 1 point for pivot part of LU
% 1 point for rest of LU
% 1 point for getting u and v
% 1 point for d

% Up to -1 for lack of vectorization.

