function Apr22()
% Fitting elliptical data
% Parameters of the "true" ellipse...
a = 5; b = 3; h = 1; k = 2; phi = pi/6;

close all
n = 20;
for delta = [0 .1 .15 .2 .3]
    for tau = [1 .7 .5 .4  .3  .2]*pi
        tMin = pi-tau;
        tMax = pi+tau;
        [x,y] = MakeEllipseData(a,b,h,k,phi,n,delta,tMin,tMax);
        % The following statement puts a data point "opposite" the cluster
        % [xn,yn] = MakeEllipseData(a,b,h,k,phi,1,delta,0,0); x(n)=xn;y(n)=yn;
        plot(x,y,'.r','markersize',20)
        hold on
        PlotEllipse(a,b,h,k,phi)
        DisplayFit(x,y)
        hold off
        title(sprintf('delta = %5.2f',delta),'fontsize',16)
        shg
        pause
    end
end

function DisplayFit(x,y)
v = axis;
xmin = v(1)-.5; xmax = v(2)+.5;
ymin = v(3)-.5; ymax = v(4)+.5;
plot([xmin xmax],[0 0],'--k',[0 0],[ymin ymax],'--k')
axis equal
[B,C,D,E,F] = ConicLS(x,y);
N = 100;
x = linspace(xmin,xmax,N);
y = linspace(ymin,ymax,N);
Z = zeros(N,N);
for i=1:N
    for j=1:N
        Z(i,j) = x(j)^2 + B*x(j)*y(i) + C*y(i)^2 + D*x(j) + E*y(i) + F;
    end
end
contour(x,y,Z,[0 0],'b','linewidth',2)



  function [B,C,D,E,F] = ConicLS(x,y)
% x and y are column n-vectors
% Determines real numbers B,C,D,E,F so that
%     sum_{i=1}^{n} (x(i)^2 + B*x(i)y(i) + C*y(i)^2 + Dx(i) + Ey(i) + F)^2
% is minimized.
n = length(x);
M = [x.*y  y.^2  x y ones(n,1)];
alfa = M\(-x.^2);
B = alfa(1); C = alfa(2); D = alfa(3); E = alfa(4); F = alfa(5);

  function [x,y] = MakeEllipseData(a,b,h,k,phi,n,delta,tMin,tMax)
% Generates n points on E(a,b,h,k,phi) and adds O(delta) noise
% The range of t values is defined by tmin and tmax, 0<=tmin<tmax<=2pi
tVals = linspace(tMin,tMax,n)';
x1 = a*cos(tVals);
y1 = b*sin(tVals);
x = h+cos(phi)*x1-sin(phi)*y1 + delta*(-1+2*rand(n,1));
y = k+sin(phi)*x1+cos(phi)*y1 + delta*(-1+2*rand(n,1));

  function PlotEllipse(a,b,h,k,phi)
% Assumes hold is on and displays E(a,b,h,k,phi)
tVals = linspace(0,2*pi,200);
W = [cos(phi) -sin(phi); sin(phi) cos(phi)]*[a*cos(tVals);b*sin(tVals)];
x = h+W(1,:);
y = k+W(2,:);
plot(x,y,'r')