  function P2()
% Checks out DistE2Pts
close all
% Generate and display the data set...
n = 100;
% The ellipse ((x-h0)/a)^2 + ((y-k0)/b)^2 = 1 rotated phi0 radians
% counterclockwise...
a0 = 5; 
b0 = 3; 
phi0 =30*(pi/180); 
h0 = 2; 
k0 = -3;
% Add this much noise to each point...
delta = .2;
[x,y] = MakeEllipseData(a0,b0,h0,k0,phi0,n,delta);
% Display the point set...
plot(x,y,'.k')

% Generate and display a reasonable fitting ellipse...
h = mean(x);
k = mean(y);
distVec = sqrt((x-h).^2+(y-k).^2);
[a,idx] = max(distVec);
b = min(distVec);
phi = atan2(y(idx)-k,x(idx)-h);
hold on
PlotEllipse(a,b,h,k,phi)
m = max(a,b);
axis([h-1.1*m h+1.1*m k-1.1*m k+1.1*m])
axis equal
shg
% Assess the quality of the fit...
z = DistE2Pts(x,y,a,b,h,k,phi);
title(sprintf('DistE2Pts = %6.3f',z),'fontsize',18)


  function [x,y] = MakeEllipseData(a,b,h,k,phi,n,delta)
% Generates N points on E(a,b,h,k,phi) and adds O(delta) noise
tVals = linspace(0,2*pi,n)';
x1 = h+a*cos(tVals)+delta*(-1+2*rand(n,1));
y1 = k+b*sin(tVals)+delta*(-1+2*rand(n,1));
x = cos(phi)*x1-sin(phi)*y1;
y = sin(phi)*x1+cos(phi)*y1;

  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')


       



