function XYproj = ProjectPointsOntoConicByWEPeq3Schur(XY,A)
%
%  The WEP method
%  version: solving the cubic equation by MATLAB function "roots"
%
%   Projecting a given set of points onto a quadratic curve (conic)
%   and computing the distances from the points to the conic

%   The conic is defined by quadratic equation 
%      Ax^2 + 2Bxy + Cy^2 + 2Dx + 2Ey + F = 0

%   Input:  A = (A,B,C,D,E,F)' is the column vector of the parameters of the conic
%           XY(n,2) is the array of coordinates of n points x(i)=XY(i,1), y(i)=XY(i,2)

%   Output:  XYproj is the array of coordinates of projections
%
%        Nikolai Chernov,  August 2012

n = size(XY,1);       %  n is the number of points to be projected
XYproj = zeros(n,2);  %  prepare an array of projections
xy = zeros(4,2);

C0 = [A(1) A(2) A(4);
      A(2) A(3) A(5);
      A(4) A(5) A(6)];   %    3x3 matrix of conic's parameters
 
 for i=1:n
     u = XY(i,1);  v = XY(i,2);   %  coordinates of the point to be projected
     
                  %  Next compute D, the auxiliary matrix 
                  %  it comes from the orthogonality condition
                  
     D = [2*A(2) A(3)-A(1) A(1)*v-A(2)*u+A(5);
         A(3)-A(1) -2*A(2) A(2)*v-A(3)*u-A(4);
         A(1)*v-A(2)*u+A(5) A(2)*v-A(3)*u-A(4) 2*A(4)*v-2*A(5)*u];
     C = C0;

     cc = det(C);
     dd = det(D);
     if (abs(cc) > abs(dd))
         G = C;  C = D;  D = G;
         cd = cc;  cc = dd;  dd = cd;
     end
     cd = det(C+D);
     dc = det(C-D);
     a = -dd;   b = (cd+dc)/2-cc;   c = (dc-cd)/2+dd;   d = cc;
     E = roots([a,b,c,d]);    %  solving cubic equation
     
     jreal=0;  jreals=0;
     for j=1:3
         if isreal(E(j))
             jreal=j;
             jreals = jreals+1;
         end
     end
     if jreals==0
         disp(' No real e-values...');
         return;
     end
     if jreals<=2
         x = E(jreal);
     end
     if jreals==3
         x1 = E(1);
         x2 = E(2);
         x3 = E(3);
         x12 = abs(x1-x2);
         x13 = abs(x1-x3);
         x23 = abs(x2-x3);
         if ((x12<=x13)&&(x12<=x23))
             x = x3;
         end;
         if ((x13<=x12)&&(x13<=x23))
             x = x2;
         end;
         if ((x23<=x12)&&(x23<=x13))
             x = x1;
         end;
     end

     G = C-x*D;    %  parameters of the degenerate conic

     if det(G(1:2,1:2))>realmin
         det(G(1:2,1:2))
         disp(' Positive 2x2 determinant...');
         return;
     end

     [Q1,Q2] = GetLinesFromDegenerateConic(G);

     %   equation of the first  line is Q1(1)*x + Q1(2)*y + Q1(3) = 0
     %   equation of the second line is Q2(1)*x + Q2(2)*y + Q2(3) = 0

     %   Next find all the points of intersection of the above lines
     %   with the original conic (defined by the matrix C)

     %   The first line gives us up to two points

     k = 0;

     if abs(Q1(1)) > max(abs(Q1(2)),realmin)
         p = Q1(2)/Q1(1);  q = Q1(3)/Q1(1);
         a = A(1)*p*p-2*A(2)*p+A(3);
         b = A(2)*q-A(1)*p*q+A(4)*p-A(5);
         c = A(1)*q*q-2*A(4)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,2) = c/2/b;
             xy(k,1) = -p*xy(k,2) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,2) = (b + sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 else
                     k = k + 1;
                     xy(k,2) = (b - sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 end
             end
         end
     end

     if abs(Q1(2)) >= max(abs(Q1(1)),realmin)
         p = Q1(1)/Q1(2);  q = Q1(3)/Q1(2);
         a = A(3)*p*p-2*A(2)*p+A(1);
         b = A(2)*q-A(3)*p*q+A(5)*p-A(4);
         c = A(3)*q*q-2*A(5)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,1) = c/2/b;
             xy(k,2) = -p*xy(k,1) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,1) = (b + sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 else
                     k = k + 1;
                     xy(k,1) = (b - sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 end
             end
         end
     end

     %   The second line gives us up to two points

     if abs(Q2(1)) > max(abs(Q2(2)),realmin)
         p = Q2(2)/Q2(1);  q = Q2(3)/Q2(1);
         a = A(1)*p*p-2*A(2)*p+A(3);
         b = A(2)*q-A(1)*p*q+A(4)*p-A(5);
         c = A(1)*q*q-2*A(4)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,2) = c/2/b;
             xy(k,1) = -p*xy(k,2) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,2) = (b + sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 else
                     k = k + 1;
                     xy(k,2) = (b - sqrt(disc))/a;
                     xy(k,1) = -p*xy(k,2) - q;
                     k = k + 1;
                     xy(k,2) = c/a/xy(k-1,2);
                     xy(k,1) = -p*xy(k,2) - q;
                 end
             end
         end
     end

     if abs(Q2(2)) >= max(abs(Q2(1)),realmin)
         p = Q2(1)/Q2(2);  q = Q2(3)/Q2(2);
         a = A(3)*p*p-2*A(2)*p+A(1);
         b = A(2)*q-A(3)*p*q+A(5)*p-A(4);
         c = A(3)*q*q-2*A(5)*q+A(6);
         if a == 0
             k = k + 1;
             xy(k,1) = c/2/b;
             xy(k,2) = -p*xy(k,1) - q;
         else
             disc = b*b-a*c;
             if disc >= 0
                 if b > 0
                     k = k + 1;
                     xy(k,1) = (b + sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 else
                     k = k + 1;
                     xy(k,1) = (b - sqrt(disc))/a;
                     xy(k,2) = -p*xy(k,1) - q;
                     k = k + 1;
                     xy(k,1) = c/a/xy(k-1,1);
                     xy(k,2) = -p*xy(k,1) - q;
                 end
             end
         end
     end

     %       Select the intersection point closest to the point (u,v)
     if k > 0
         jmin = 1;
         dmin = realmax;
         for j=1:k
             if norm(xy(j,:) - XY(i,:))<dmin
                 dmin = norm(xy(j,:) - XY(i,:));
                 jmin = j;
             end
         end
         XYproj(i,:) = xy(jmin,:);   %    record the found projection
     else
         XYproj(i,:)=XY(i,:);
     end

     if ((~isfinite(XYproj(i,1)))||(~isfinite(XYproj(i,2))))
         XYproj(i,:) = XY(i,:);
     end
 end
 
end   %  ProjectPointsOntoConicByWEPeq3Schur