  function z = DistE2Pts(u,v,a,b,h,k,phi)
% u and v are column n-vectors
% z = d(u(1),v(1),a,b,h,k,phi) +...+ d(u(n),v(n),a,b,h,k,phi)

% Translate and rotate the data so that the problem is equivalent to
% looking at the distance to E(a,b,0,0,0)...
u = u-h;
v = v-k;
W = [cos(phi) -sin(phi); sin(phi) cos(phi)]'*[u';v'];
u = W(1,:)';
v = W(2,:)';
z = 0;
for i=1:length(u)
    if u(i)>=0 && v(i)>=0
        % First quadrant...
        L = 0; R = pi/2;
    elseif u(i)<0 && v(i)>=0
        % Second quadrant...
        L = pi/2; R = pi;
    elseif u(i)<=0 && v(i)<=0
        % Third quadrant...
        L = pi; R = 3*pi/2;
    else
        % Fourth quadrant...
        L = 3*pi/2;R = 2*pi;
    end
    [tStar,dStar] = fminbnd(@(t) d(t,u(i),v(i),a,b),L,R);
    z = z+dStar;
end

  function z = d(t,alpha,beta,a,b)
% Objective function for fminbnd
z = (alpha - a*cos(t))^2 + (beta - b*sin(t))^2;
  
% -2 for optimizing over wrong section of ellipse or entire ellipse
% -2 for incorrect transformations
% -1 for doing more work than necessary (computing inverse tangents, doing
%       extra transformations, etc.)