%
% draws an ellipse with ax^2+bxy+cy^2+dx+ey+f=0
% overlaid with original points (rx,ry)
%
function drawellip(a,rx,ry)

   %  convert to standard form: ((x-cx)/r1)^2 + ((y-cy)/r2)^2 = 1
   % rotated by theta
   v = solveellipse(a);

   % draw ellipse with N points   
   N = 100;
   dx = (2*v(1))/N;
   theta = v(5);
   R = [ [ cos(theta) -sin(theta)]', [sin(theta) cos(theta)]'];
   for i = 1:N+1
        x = -v(1) + (i-1)*dx;
        y = v(2)*sqrt(1 - x*x/(v(1)*v(1)));
        d1 = R*[x y]';
        d2 = R*[x -y]';
        Xt(i) = d1(1) + v(3);
        Yt(i) = d1(2) + v(4);
        Xb(i) = d2(1) + v(3);
        Yb(i) = d2(2) + v(4);
   end

figure(1)
hold off
plot(rx,ry,'r')
hold on
plot(Xt,Yt,'g')
plot(Xb,Yb,'b')

