function out=gs(A)
%
% GS(A) where A is a 3x3 real matrix preferrably with 
% entries between -2 and 2. Self-demonstrating with a
% reasonable, typical example.
% GS graphically demonstrates the Gram-Schmidt orthogo-
% nalization algorithm.
%
% Original: August 1999
% Updated November 2003: Fatter lines for better visibility.
% All rights reserved. Matt Kawski. 
% http://math.asu.edu/~kawski
a=2
if nargin < 1
   A=[3/2 -3/2 0; 0 -2 -1/3;-1/2 0 2]'
   % A=a*rand(3,3);
end
clf
v1=A(:,1)
v2=A(:,2)
v3=A(:,3)
clf
axis(a*[-1 1 -1 1 -1 1])
hold on
plot3([0,v1(1)],[0,v1(2)],[0,v1(3)],'b-','LineWidth',3)
plot3([0,v2(1)],[0,v2(2)],[0,v2(3)],'b-','LineWidth',3)
plot3([0,v3(1)],[0,v3(2)],[0,v3(3)],'b-','LineWidth',3)
axis(a*[-1 1 -1 1 -1 1])
hold on
disp('Strike any key to continue')
pause
w1=v1
R11=sqrt(w1'*w1)
Q1=w1/R11
clf
axis(a*[-1 1 -1 1 -1 1])
hold on
plot3([0,v2(1)],[0,v2(2)],[0,v2(3)],'b-','LineWidth',3)
plot3([0,v3(1)],[0,v3(2)],[0,v3(3)],'b-','LineWidth',3)
plot3([0,Q1(1)],[0,Q1(2)],[0,Q1(3)],'r-','LineWidth',3)
if R11>1
   plot3([Q1(1),v1(1)],[Q1(2),v1(2)],[Q1(3),v1(3)],'b-','LineWidth',3)
end
disp('Strike any key to continue')
pause
R12=Q1'*v2
p1=R12*Q1
w2=v2-p1
R22=sqrt(w2'*w2)
Q2=w2/R22
plot3([p1(1),v2(1)],[p1(2),v2(2)],[p1(3),v2(3)],'g-','LineWidth',3)
if R12 > max(R11,1)
   plot3([max([Q1(1),v1(1)]),p1(1)],...
      [max([Q1(2),v1(2)]),p1(2)],...
      [max([Q1(3),v1(3)]),p1(3)],'g-','LineWidth',3)
end
if R12 < 0 
   plot3([0,p1(1)],[0,p1(2)],[0,p1(3)],'g-','LineWidth',3)
end
disp('Strike any key to continue')
pause
plot3([w2(1),v2(1)],[w2(2),v2(2)],[w2(3),v2(3)],'g-','LineWidth',3)
plot3([0,Q2(1)],[0,Q2(2)],[0,Q2(3)],'r-','LineWidth',3)
if R22>1
   plot3([Q2(1),w2(1)],[Q2(2),w2(2)],[Q2(3),w2(3)],'b-','LineWidth',3)
end
disp('Strike any key to continue')
pause

%tt=[-2:1:2];
%[S,T]=meshgrid(tt,tt);
%mesh(S*Q1(1)+T*Q2(1),S*Q1(2)+T*Q2(2),S*Q1(3)+T*Q2(3),ones(size(S)))
%hidden off

nn = 2
for k=[-nn:0.2:nn]
    plot3([k*Q1(1)-nn*Q2(1),k*Q1(1)+nn*Q2(1)],...
          [k*Q1(2)-nn*Q2(2),k*Q1(2)+nn*Q2(2)],...
          [k*Q1(3)-nn*Q2(3),k*Q1(3)+nn*Q2(3)],'c-')
    plot3([-nn*Q1(1)+k*Q2(1),nn*Q1(1)+k*Q2(1)],...
          [-nn*Q1(2)+k*Q2(2),nn*Q1(2)+k*Q2(2)],...
          [-nn*Q1(3)+k*Q2(3),nn*Q1(3)+k*Q2(3)],'c-')
end
plot3([0,Q1(1)],[0,Q1(2)],[0,Q1(3)],'r-','LineWidth',3)
plot3([0,Q2(1)],[0,Q2(2)],[0,Q2(3)],'r-','LineWidth',3)
disp('Strike any key to continue')
pause

R13=Q1'*v3
R23=Q2'*v3
p2=R13*Q1+R23*Q2
w3=v3-p2
R33=sqrt(w3'*w3)
Q3=w3/R33
plot3([p2(1),v3(1)],[p2(2),v3(2)],[p2(3),v3(3)],'g-','LineWidth',3)
plot3([0,p2(1)],[0,p2(2)],[0,p2(3)],'g-','LineWidth',3)
disp('Strike any key to continue')
pause
plot3([w3(1),v3(1)],[w3(2),v3(2)],[w3(3),v3(3)],'g-','LineWidth',3)
if R33>1
   plot3([Q3(1),w3(1)],[Q3(2),w3(2)],[Q3(3),w3(3)],'b-','LineWidth',3)
end
plot3([0,Q3(1)],[0,Q3(2)],[0,Q3(3)],'r-','LineWidth',3)
plot3([0,Q2(1)],[0,Q2(2)],[0,Q2(3)],'r-','LineWidth',3)
plot3([0,Q1(1)],[0,Q1(2)],[0,Q1(3)],'r-','LineWidth',3)
disp('Strike any key to continue')
pause
for k=[0:3:180], view(k,20), pause(1/2), end