function [d,u,v] = InsideDist(x,y,theta);

n = length(x); x = [x;x(1)]; y = [y;y(1)];

% Set up the n linear systems so that we do not have to pivot.
% Note that all the systems have the same first column.

c = cos(theta); s = sin(theta);
if abs(c)>=abs(s);
% No pivoting required in any system...
A12 = x(1:n)-x(2:n+1);
A22 = y(1:n)-y(2:n+1);
b1 = x(1:n);
b2 = y(1:n);
m = s/c;
A11 = c;
else
% Pivoting required in every system, so set up interchanged systems...
A22 = x(1:n)-x(2:n+1);
A12 = y(1:n)-y(2:n+1);
b2 = x(1:n);
b1 = y(1:n);
m = c/s;
A11 = s;
end
% Eliminate the first unknown from the second equation...
A22 = A22 - m*A12; b2 = b2 - m*b1;

% We now have a triangular system which can be back solved...
t2 = b2./A22;
t1 = (b1 - A12.*t2)./A11;

% Assemble u and v. Involves finding the indices where
% the solution components fall in the right range...
idx = find(0<=t1 & 0<=t2 & t2<=1);
u = cos(theta)*t1(idx); v = sin(theta)*t1(idx);

% Compute inside distances...
m = length(u);
if m==0
d = 0;
else
% Compute the distances of each intersection point to origin and sort..
dP = sort(sqrt(u.^2 + v.^2));
d = sum(dP(2:2:m)-dP(1:2:m));
end