function [v,d]=mtcee(M,plotwho)

% [v,d]=mtcee(M,plotwho)
% matrix times circle equals ellipse, M is a 2x2 matrix
% plotwho is optional, if you enter 1 the standard basis vectors (the coordinate 
% axes) and their images are plotted, if you enter 2 the eigenvectors are plotted
% (the default is 1)
% The outputs are the eigenvectors, v (they are scaled by their corresponding
% eigenvalues), and eigenvalues, d 

% The idea is that a 2-long vector can be thought of as a point in the plane. In
% this program we start with a bunch of points on the unit circle and see where 
% they move to when we multiply them by M (this is called the "image" of the unit
% circle and this new set of points is in fact an ellipse!)

% technical point: the eigenvalues and eigenvectors can be complex valued. 
% In this case only the basis vectors will be plotted, not the eigenvectors.

% What can you learn from these pictures? If the matrix is symmetric the 
% eigenvectors are orthogonal (perpendicular) and form the axes of the ellipse. 
% The absolute value of the determinant is the area of the ellipse and if the 
% determinant is negative then the ellipse has been reflected as well (indicated
% by a red eigenvector corresponding to a negative eigenvalue). The determinant
% is the product of the eigenvalues. An orthogonal matrix only reflects and/or 
% rotates the circle, no stretching. Pure imaginary eigenvalues (no real part)
% will only rotate the circle and rescale it's radius by a constant.
% Try some decompostitions like LU and QR and see what the factors look like. 
% The "action" of the matrix M will be the "composition" of the actions of it's
% factors. Try M^2, M^3, M^4...

if (nargin==1) | ~(plotwho==2)
    plotwho=1;
end

n=64;
k=[0:n];
c=[cos(2*pi*k/n);sin(2*pi*k/n)];
e=M*c;
x=[1,0,-1,0;0,1,0,-1];
y=M*x;

[v,d]=eig(M);
d=diag(d)';
if isreal(d)
    v1=v(:,1)*d(1);
    v2=v(:,2)*d(2);
    v=[v1,v2];
else
    plotwho=1;
end

figure
hold on
plot(c(1,:),c(2,:))
plot(e(1,:),e(2,:),'g')
if (plotwho==2)
    c1='k'; c2='k';
    if d(1)<0
        c1='r';
    end
    if d(2)<0
        c2='r';
    end
    plot([0,v1(1)],[0,v1(2)],c1);
    plot([0,v2(1)],[0,v2(2)],c2);

else
    xcs=reshape([zeros(1,4);x(1,:)],1,8);
    x2cs=reshape([zeros(1,4);x(2,:)],1,8);
    plot(xcs(1:2),x2cs(1:2),'r')
    plot(xcs(3:4),x2cs(3:4),'k')

    ycs=reshape([zeros(1,4);y(1,:)],1,8);
    y2cs=reshape([zeros(1,4);y(2,:)],1,8);
    plot(ycs(1:2),y2cs(1:2),'r')
    plot(ycs(3:4),y2cs(3:4),'k')
end
axis equal
hold off