clf, subplot(1,2,1), plot(0,0), hold on, axis([-2,2,-2,2]), axis manual
for i=1:4
    [x(i),y(i)]=ginput(1)
    if mod(i,2) 
         plot(x(i),y(i),'*') 
    else 
         plot(x(i),y(i),'r*') 
    end
end 

plot([x(1),x(3)],[y(1),y(3)],'LineWidth',2)
plot([x(2),x(4)],[y(2),y(4)],'r','LineWidth',2)

C=[x(1) y(1) 0    0;
   0    0    x(1) y(1);
   x(3) y(3) 0    0;
   0    0    x(3) y(3)]

b=[x(2) y(2) x(4) y(4)]'

s=inv(C)*b
A=reshape(s,2,2)'
%collect lines by threads
K = 5;
X = linspace(x(1), x(3), K);
Y = linspace(y(1), y(3), K);
R = A*[X;Y];
for i=1:K
   plot([X(i),R(1,i)],[Y(i),R(2,i)],'k','LineWidth',1)
end

test1to2=A*[x(1);y(1)]-[x(2);y(2)]
test3to4=A*[x(3);y(3)]-[x(4);y(4)]

%A = randn(2,2); %random matrix, Gaussian sampling
subplot(1,2,2)
axis equal
I = eye(2,2);
N = 500;
alpha = 0:2*pi/N:2*pi; %cut 2pi into N pieces

U = [cos(alpha); sin(alpha)]; % unit circle
resultA = A*U; %transforming unit circle by matrix A
resultI = I*U; %ditto for I

xA = resultA(1,:); yA = resultA(2,:); %splits results into components
xI = resultI(1,:); yI= resultI(2,:);

%hold off;
plot(xA,yA,'r','LineWidth', 2);
hold on;
plot(xI,yI,'b','LineWidth', 2);

M = 50;
alpha = 0:2*pi/M:2*pi; %cut 2pi into M pieces
U = [cos(alpha); sin(alpha)];
resultA = A*U;
for i=1:M %create pairs of points on circle and A*circle
    plot([U(1,i) resultA(1,i)],[U(2,i) resultA(2,i)], 'k','LineWidth',1);
end