% square.m
clf; hold on;
axis([0 10 0 10])
axis('equal')
p=100; q=100; r=100; s=100;

disp('arrange the sides of the square counterclockwise')
disp('enter points P_i belonging to side A')
disp('by clicking the mouse in the graphical window.')
disp('Finish the input by pressing the Return key')
[Px,Py] = ginput(p); plot(Px,Py,'o')
disp('enter points Q_i for side B ')
[Qx,Qy] = ginput(q); plot(Qx,Qy,'x')
disp('enter points R_i for side C ')
[Rx,Ry] = ginput(r); plot(Rx,Ry,'*')
disp('enter points S_i for side D ')
[Sx,Sy] = ginput(s); plot(Sx,Sy,'+')

zp = zeros(size(Px)); op =  ones(size(Px));
zq = zeros(size(Qx)); oq =  ones(size(Qx));
zr = zeros(size(Rx)); or =  ones(size(Rx));
zs = zeros(size(Sx)); os =  ones(size(Sx));

A = [ zp op zp  Px  Py
      zq zq oq  Qy -Qx
      or or zr  Rx  Ry
      os zs os  Sy -Sx]
[cc, n] = clsq(A,2)
% compute the vector of constants for the four lines
c(1) = cc(2); c(2) = cc(3); c(3) = c(1)+cc(1); 
c(4) = c(2)+cc(1);  c=c(:);
% compute the 4 corners of the square 
B  = [n  [-n(2) n(1)]']
X = -B* [c([1 3 3 1])'; c([2 2 4 4])']
X = [X X(:,1)]
plot(X(1,:), X(2,:))
hold off;

