function quads()
% Calculates alternative conformations
% of a system of conjoined articulated
% quadrilaterals
%   Written by Vangelis Coutsias (C), 
%   Mathematics Department, UNM
%   
% Based on Coutsias et al., 
%     JQuantumChem 106, 176-189 (2006) 
close all

 %non-degenerate 
 [a,b,c,d] = specs(2)

%------------------------------------
%              c1
%   |------------------|
%   |   c2    (1)      |
% d1|---------|   c3   |b1
%   |d2     b2|--------|
%   |   (2)   |d3(3) b3|
%   |     a2  |t2  a3  |t3
%   |---------|--------|---
%   |t1         a1       
% t1, t2, t3 halftangents of exterior angles 
%------------------------------------
% poly coeffs, Bricard notation
% R. Bricard, JMathPurAp 3, 113-150(1987)
A(1) = (a(1)-b(1)-d(1)+c(1))*(a(1)-b(1)-d(1)-c(1));
B(1) = (a(1)-b(1)+d(1)+c(1))*(a(1)-b(1)+d(1)-c(1));
C(1) = -4*b(1)*d(1);
D(1) = (a(1)+b(1)-d(1)+c(1))*(a(1)+b(1)-d(1)-c(1));
E(1) = (a(1)+b(1)+d(1)+c(1))*(a(1)+b(1)+d(1)-c(1));

A(2) = (a(2)-b(2)-d(2)+c(2))*(a(2)-b(2)-d(2)-c(2));
B(2) = (a(2)+b(2)-d(2)+c(2))*(a(2)+b(2)-d(2)-c(2));
C(2) = -4*b(2)*d(2);
D(2) = (a(2)-b(2)+d(2)+c(2))*(a(2)-b(2)+d(2)-c(2));
E(2) = (a(2)+b(2)+d(2)+c(2))*(a(2)+b(2)+d(2)-c(2));

A(3) = (a(3)-b(3)+d(3)+c(3))*(a(3)-b(3)+d(3)-c(3));
B(3) = (a(3)+b(3)+d(3)+c(3))*(a(3)+b(3)+d(3)-c(3));
C(3) = -4*b(3)*d(3);
D(3) = (a(3)-b(3)-d(3)+c(3))*(a(3)-b(3)-d(3)-c(3));
E(3) = (a(3)+b(3)-d(3)+c(3))*(a(3)+b(3)-d(3)-c(3));
%-------------------------------------------------------
% A1*t3^2*t1^2 + B1*t3^2 + 2*C1*t3*t1 + D1*t1^2 + E1 = 0
% A2*t1^2*t2^2 + B2*t1^2 + 2*C2*t1*t2 + D2*t2^2 + E2 = 0
% A3*t2^2*t3^2 + B3*t2^2 + 2*C3*t2*t3 + D3*t3^2 + E3 = 0
%-------------------------------------------------------
P = zeros(3,3,3);

P(1:3,1,1) =   E(1:3);
P(1:3,1,3) =   D(1:3);
P(1:3,2,2) = 2*C(1:3);
P(1:3,3,1) =   B(1:3);
P(1:3,3,3) =   A(1:3);

%-----------------------------------------
% perturbation to assess pseudospectrum
w = 0.00000*rand(3,3,3);
P = P + w;
%-----------------------------------------
% Dixon resultant matrices
for k = 1:3
    L(k,1:3) = P(1,1:3,k);
    M(k,1:3) = P(2,k,1:3);
    N(k,1:3) = P(3,k,1:3);
end
for i = 1:3
    A(i,1:3) = M(i,2)*N(1,1:3) - M(i,1)*N(2,1:3);
    B(i,1:3) = M(i,3)*N(1,1:3) - M(i,1)*N(3,1:3);
    C(i,1:3) = M(i,3)*N(2,1:3) - M(i,2)*N(3,1:3);
    D(i,1:3) = L(i,1:3);
end
% The Dixon resultant is the determinant of the matrix polynomial
% D(t3) := R{1} + R{2} * t3 + R{3} * t3^2
for i = 1:3
    R{i} = [0      A(1,i) A(2,i) A(3,i) 0      B(1,i) B(2,i) B(3,i);...
            A(1,i) A(2,i) A(3,i) 0      B(1,i) B(2,i) B(3,i) 0     ;...
            0      B(1,i) B(2,i) B(3,i) 0      C(1,i) C(2,i) C(3,i);...
            B(1,i) B(2,i) B(3,i) 0      C(1,i) C(2,i) C(3,i) 0     ;...
            0      0      0      0      0      D(1,i) D(2,i) D(3,i);...
            0      0      0      0      D(1,i) D(2,i) D(3,i) 0     ;...
            0      D(1,i) D(2,i) D(3,i) 0      0      0      0     ;...
            D(1,i) D(2,i) D(3,i) 0      0      0      0      0     ];
end
%--------------------------------------
% generalized eigenproblem
ZZ = zeros(8,8);
II =   eye(8,8);
AA = [ ZZ II; -R{1} -R{2}]
BB = [ II ZZ;  ZZ    R{3}]
% eigenvalues  are t3 
% while t2, t3 are found in components 2 and 5 
% of the corresponding eigenvector
%                      *            *   
% eigenvectors are [1 t1 t1^2 t1^3 t2 t1*t2 t1^2*t2 t1^3*t2]'
[VV, DD] = eig(AA,BB);
%---------------------------------------
% Reconstruction of TriQuad system
% identify purely real eigenvalues
% only consider positive real parts
u = 0*rand(16,2); j = 0;
for i = 1:16
    dd(i) = DD(i,i);
    x(i)=real(dd(i))+u(i,1);
    y(i)=imag(dd(i))+u(i,2);
    %if (y(i) == 0)
   %    if (x(i) > 0 & x(i) ~= Inf)
   %    if (x(i) ~= Inf |x(i) ~= -Inf | y(i) ~= Inf| y(i) ~=-Inf)
            j = j+1;
            t3(j) = dd(i)%x(i);
            t1(j) = VV(2,i)/VV(1,i);
            t2(j) = VV(5,i)/VV(1,i);
            %    end
        %end
    %if (y(i) > 0)
   %    if (x(i) > 0 & x(i) ~= Inf)
   %    if (x(i) ~= Inf |x(i) ~= -Inf | y(i) ~= Inf| y(i) ~=-Inf)
   %        j = j+1;
   %        t3(j) = dd(i);%x(i);
   %        t1(j) = VV(2,i)/VV(1,i);
   %        t2(j) = VV(5,i)/VV(1,i);
   %         j = j+1;
   %        t3(j) = dd(i+1);%x(i);
   %        t1(j) = VV(2,i+1)/VV(1,i+1);
   %        t2(j) = VV(5,i+1)/VV(1,i+1);
   %    end
   %end
end
cos1 = (1-t1.^2)./(1+t1.^2);
sin1 = 2*t1./(1+t1.^2);
cos2 = (1-t2.^2)./(1+t2.^2);
sin2 = 2*t2./(1+t2.^2);
cos3 = (1-t3.^2)./(1+t3.^2);
sin3 = 2*t3./(1+t3.^2);
for i = 1:j
    Pa =           [0 0];
    Pb =      a(1)*[1 0];
    Pe =      a(2)*[1 0];
    Pg =      d(2)*[-cos1(i) sin1(i)];
    Pd =      d(1)*[-cos1(i) sin1(i)];
    Ph = Pe + d(3)*[ cos2(i) sin2(i)];
    Pf = Pe + b(2)*[ cos2(i) sin2(i)];
    Pl = Pb + b(3)*[ cos3(i) sin3(i)];
    Pc = Pb + b(1)*[ cos3(i) sin3(i)];
    if (y(i) == 0)
        figure
        hold on
        xx1 = [Pa(1),Pb(1),Pc(1),Pd(1),Pa(1)];
        yy1 = [Pa(2),Pb(2),Pc(2),Pd(2),Pa(2)];
        xx2 = [Pa(1),Pe(1),Pf(1),Pg(1),Pa(1)];
        yy2 = [Pa(2),Pe(2),Pf(2),Pg(2),Pa(2)];
        xx3 = [Pe(1),Pb(1),Pl(1),Ph(1),Pe(1)];
        yy3 = [Pe(2),Pb(2),Pl(2),Ph(2),Pe(2)];
        plot(xx1,yy1,'b',xx2,yy2,'r--',xx3,yy3,'g-.')
        hold off
    end
    % Verify constraints: s5, s6, s4 must be correct
    i
    error1(i) = c(1) - norm(Pc - Pd);
    error2(i) = c(2) - norm(Pf - Pg);
    error3(i) = c(3) - norm(Pl - Ph);
    fprintf(' closure errors \n');
    fprintf('err1= %g  err2= %g  err3= %g\n',...
             error1(i),error2(i),error3(i));
end
j
if (j > 0)
    fprintf(' angles t1, t2, t3 \n');
    i=1;i
fprintf('%g ',t1(1:j));
    i=2;i
fprintf('%g ',t2(1:j));
    i=3;i
fprintf('%g ',t3(1:j));
% compute Jacobian and its rank to test
% flexibility  J(i,j) = d(ci^2)/dtj
J = zeros(1:3,1:3,1:16);
J(1,1,1:16) = A(1)*t1.*t3.^2 + D(1)*t1 + C(1)*t3;
J(2,1,1:16) = A(2)*t1.*t2.^2 + B(2)*t1 + C(2)*t2;
J(3,2,1:16) = A(3)*t2.*t3.^2 + B(3)*t2 + C(3)*t3;
J(2,2,1:16) = A(2)*t1.^2.*t2 + D(2)*t2 + C(2)*t1;
J(3,3,1:16) = A(3)*t2.^2.*t3 + D(3)*t3 + C(3)*t2;
J(1,3,1:16) = A(1)*t1.^2.*t3 + B(1)*t3 + C(1)*t1;
J(:,:,4)
Jacob(1:16) = (J(1,1,1:16).*J(2,2,1:16).*J(3,3,1:16)+...
               J(1,3,1:16).*J(2,1,1:16).*J(3,2,1:16))
Jacob(1:16) = Jacob(1:16)./((1+t1(1:16).^2).*(1+t2(1:16).^2).*(1+t3(1:16).^2)).^2;
         fprintf('\n---------------\n');
fprintf(' Jacobian \n');
for i = 1:j
fprintf('%g   %g   %g\n',i,real(Jacob(i)), imag(Jacob(i)));
end
end
fprintf('\n---------------\n');
fprintf(' eigenvalues \n');
for i = 1:16
fprintf('%g   %g   %g\n',i,real(dd(i)),imag(dd(i)));
end
figure
plot(x,y,'ro')


function [a,b,c,d] = specs(cases)
if (cases == 1) 
 %non-degenerate 
    a(2)  = 1+13/16;%  19;
    a(3)  =   14/16;%   9.2;
    a(1)  = a(2) + a(3);
    d(1) = 1+14/16;%  18.1;
    b(2) = 1+ 2/16;%   4.1; 
    b(3) = 1+ 3/16;%   9.2;
    c(1) = 1+13/16;%  18.1;
    c(2) = 1+ 5/16;%   4.1;
    c(3) =   14/16;%  16;
    d(2) = 1+ 6/16;%   2.8;
    d(3) =   13/16;%  16;
    b(1) = 2+ 2/16;%  28.2;
end
%Doubly degenerate system
if (cases == 2)
    a(2)  =   19;
    a(3)  =    9.2;
    a(1)  = a(2) + a(3);
    d(1) =   18.1;
    b(2) =   19  ; 
    b(3) =    9.2;
    c(1) =   18.1;
    c(2) =    4.1;
    c(3) =   16  ;
    d(2) =    2.8;
    d(3) =  -16  ;
    b(1) =   28.2;
end
%Triply degenerate system
if (cases == 3)
    a(2)  =   19;
    a(3)  =    9.2;
    a(1)  = a(2) + a(3);
    d(1) =   18.1;
    b(2) =   19  ; 
    b(3) =    9.2;
    c(1) =   18.1;
    c(2) =    4.1;
    c(3) =   16  ;
    d(2) =    4.1;
    d(3) =  -16  ;
    b(1) =   28.2;
end
% Flexible: quad degeneracy
% (system of parallelograms)
if (cases == 4)
    a(2) = 10;%  19;
    a(3) =  3;%   9.2;
    a(1) =  a(2) + a(3);
    d(1) = 7+.1;%  18.1;
    b(2) = 4;%   4.1; 
    b(3) = 2;%   9.2;
    c(1) = 13;%  18.1;
    c(2) = 10;%   4.1;
    c(3) = 3;%  16;
    d(2) = 4;%   2.8;
    d(3) = 2;%  16;
    b(1) = 7+.15;%  28.2;
end