%ComputeTrnasformParams from defined vector v
%script to have a look at three dimentional transformations of a vector to one 
%of the axes. This is the transformation required to monitor the direction vector 
%in a monte Carlo simulation.
%There are 5 presented. The first is copied from some MC demonstration code.
%the next is the theory with repeats this patern.
%the 5th one is an alternate to the first two.
%note 1,2,4 rotate to the z axis - the axis of prpagation in MC simulation
%if you are moving in a different direction use either 3 or 4.

v(1,:) = [1,-2,3];
v(2,:) = [1,2,3];
v(3,:) = [-1,2,3];
v(4,:) = [-1,-2,3];

for i = 1:size(v,1)

a = v(i,1)/norm(v(i,:));
b = v(i,2)/norm(v(i,:));
c = v(i,3)/norm(v(i,:));
d = sqrt(a^2+b^2);
e = sqrt(b^2+c^2);
f = sqrt(a^2+c^2);

directioncosines = acos(v/norm(v));

theta  =  acos(c);  %angle of vetor to z axis
phi = atan2(b,a);
if b<0
   phi = 2*pi+phi;
end

disp(['a:',num2str(a),' b:',num2str(b),' c:',num2str(c),' d:',num2str(d),...
      '  Phi:',num2str(phi*180/pi),'  Theta:',num2str(theta)]);

%Transform 1
% code rotates about z axis by phi degrees then about y axis by theta degrees
codeT = [ cos(theta)*cos(phi), -sin(phi), sin(theta)*cos(phi);
          sin(phi)*cos(theta),  cos(phi) , sin(theta)*sin(phi);
         -sin(theta)         ,  0         , cos(theta) ];
codeTinv = inv(codeT);

%Transform 2
% Theory using the coordinates of the code: (moves to z axis)
%Rz rotates about z axis -phi degrees on the xy plane
Rz2 = [cos(phi),-sin(phi),0;sin(phi),cos(phi),0;0,0,1];
%Rz2 = [sin(phi),-cos(phi),0;cos(phi),sin(phi),0;0,0,1];
%Ry rotates -theta about y axis
Ry2 = [cos(theta),0,sin(theta);0,1,0;-sin(theta),0,cos(theta)];
theoryT2=Rz2*Ry2;

%Transform 3
% Theory using the coordinates of the code: (moves to x axis)
%Rz rotates about z axis -phi degrees on the xy plane
Rz3 = [cos(phi),-sin(phi),0;sin(phi),cos(phi),0;0,0,1];
%Ry rotates (pi/2-theta) about y axis
Ry3 = [sin(theta),0,-cos(theta);0,1,0;cos(theta),0,sin(theta)];
theoryT3=Rz3*Ry3;

%Transform 4
%theory using a different axis orientation (moves to y axis).
%Rx rotates thetax about x axis
Rx4 = [1,0,0;0,b/e,-c/e;0,c/e,b/e];
%Rz rotates about z axis thetaz degrees on the xy plane
Rz4 = [e,a,0;-a,e,0;0,0,1];
theoryT4=Rx4*Rz4;

%Transform 5
%theory using a different axis orientation (moves to z axis).
%Rx rotates thetax2 about x axis to the xz plane
Rx5 = [1,0,0;0,c/e,b/e;0,-b/e,c/e];
%Rz rotates about y axis thetay2 degrees on the xz plane
Ry5 = [e,0,a;0,1,0;-a,0,e];
theoryT5=Rx5*Ry5;

v(i,:)*Rz2
v(i,:)*Rz2*Ry2


end
%usage:
%disp('theory does:');
%out=v
%out = v*Rz2
%out = v*Rz2*Ry2

