% Demonstration of mapping of sphere under a linear transformation
% using the SVD
% Input: 3-by-3 matrix a
% Output: image of unit sphere under a, and principle axes
%
hold off, clg
theta1=(0:2*pi/1000:2*pi);
mm=max(size(theta1));
circle1=[cos(theta1);sin(theta1)];
m=20;
theta2=(0:2*pi/m:2*pi);
[u,s,v]=svd(a);
s=diag(s)
nrm=max(s);
i=1;
for theta = theta2,
    circle2=v*[sin(theta)*ones(1,mm);cos(theta)*circle1];
    ellipse=a*circle2;
    plot3(ellipse(1,:),ellipse(2,:),ellipse(3,:),'.w');
    if i==1,
       i=2;
       axis(1.0*nrm*[-1,1,-1,1,-1,1]);
       axis('square')
       hold on
    end
end
plot3([0,s(1)*u(1,1)],[0,s(1)*u(2,1)],[0,s(1)*u(3,1)],'b-')
plot3([0,s(2)*u(1,2)],[0,s(2)*u(2,2)],[0,s(2)*u(3,2)],'r-')
plot3([0,s(3)*u(1,3)],[0,s(3)*u(2,3)],[0,s(3)*u(3,3)],'g-')