% reflections.m, reflections in a plane in 3D
% run this program several times, random data.
clear

v1=randn(3,1); v2=randn(3,1); A=[v1 v2];

% Construct the reflection matrix, M, and check for key properties
P=A*inv((A'*A))*A'; M=2*P-eye(3);
norm(M-M','fro'), norm(M*M-eye(3),'fro')

% Create an ON-basis for the plane and generate plane-plot
[q,r]=qr(A,0);
step=0.2; width=3; index=1;
for l=-width:step:width
	p1=l*q(:,1)-width*q(:,2);
	p2=l*q(:,1)+width*q(:,2);
	p3=(l+step/2)*q(:,1)+width*q(:,2);
	p4=(l+step/2)*q(:,1)-width*q(:,2);
	planeplot(:,index)=p1;
	planeplot(:,index+1)=p2;
	planeplot(:,index+2)=p3;
	planeplot(:,index+3)=p4;
	index=index+4;
end
ax=[ -width width -width width -width width];

% Generate some random vectors and project them
m=50;
dum=randn(3,m);
dummer=M*dum;

% plot some examples
figure(1), clf, axis(ax), axis equal, view(3), hold on
plot3([0 v1(1)], [0 v1(2)], [0 v1(3)], 'b')
plot3([0 v2(1)], [0 v2(2)], [0 v2(3)], 'b')
for l=1:5
	plot3(dum(1,l),dum(2,l),dum(3,l),'rp')
	plot3(dummer(1,l),dummer(2,l),dummer(3,l),'gp')
	plot3([dum(1,l) dummer(1,l)], [dum(2,l) dummer(2,l)], [dum(3,l) dummer(3,l)], 'y')
end
plot3(planeplot(1,:), planeplot(2,:), planeplot(3,:), 'k')
title('Five points and their reflections')

