  function [A,B,C,D,E,F,q] = IsAnEllipse(x,y)
  % x and y are column m-vectors.
  % q is the largest integer so that the best conic fit of the data
  % (x(1),y(1)),...,(x(q),y(q)) is an ellipse. A, B, C, D, E, and F are
  % the coefficients of the ellipse. If no such q exists, then q is
  % assigned the value of zero. 
  
  n = length(x);
 
  % Initializations
  q = 5;
  
  % Solve the q = 5 problem via QR...
  q = 5;
  M = [ x(1:q).*y(1:q)  y(1:q).*y(1:q)  x(1:q) y(1:q) ones(q,1)];
  b = -(x(1:q).*x(1:q));
  [Q,R] = qr(M,0);
  b = Q'*b;
  z = R\b;
  
  
  % In all that follows, we think of z = [B;C;D;E;F].
  % We always have an ellipse based on q points.
  % zNext are the coefficients for the conic based on q+1 points
  
  [R,zNext,b] = UpdateRzb(R,b,x(6),y(6));
  NextConicOK = 4*zNext(2)-zNext(1)^2 > 0;
  q = 5;
  while NextConicOK
      % Accept the next conic
      z = zNext;
      q = q+1;
      if q==n
          % Gone thru all n points
          NextConicOK = 0;
      else
          % Obtain the next conic by updating R and b
          [R,zNext,b] = UpdateRzb(R,b,x(q+1),y(q+1));
          NextConicOK = 4*zNext(2)-zNext(1)^2 > 0;
      end
  end
  A = 1; B = z(1); C = z(2); D = z(3); E = z(4); F = z(5);
      
  
  function [Rnext,znext,bnext] = UpdateRzb(R,b,nextPx,nextPy)
 % R is 5-by-5 and upper triangular
 % b is 5-by-1
 % nextPx and nextPy are scalars
 % norm(Mtilde*znext- btilde,2) = min where
 
   Mtilde = [ nextPx*nextPy  nextPy^2  nextPx  nextPy 1 ; R];
   btilde = [ -nextPx^2;b];
 
 % bnext = Q(:,1:5)'*b where  [Q,R] = qr(Mtilde,0)
 % But we handle the QR factorization ourselves thru a sequence
 % of 2-by-2 orthogonal updates of Mtilde and btilde...
   for k=1:5
       %Zero Mtilde(k+1,k)
       [LittleQ,T] = qr(Mtilde(k:k+1,k));   %  "LazyGivens"
       Mtilde(k:k+1,k:5) = LittleQ'*Mtilde(k:k+1,k:5);
       btilde(k:k+1) = LittleQ*btilde(k:k+1);
   end
   Rnext = Mtilde(1:5,1:5);
   bnext = btilde(1:5);
   znext = Rnext\bnext;
   
   
