function f = EuToRod

%initial orientation
InitOrient = EuToG(108*pi/180, 44*pi/180, 58*pi/180);  

% Data file
while(1)
   fileIn = input('Open data file (*.txt): ','s');
   fidIn = fopen(fileIn, 'rt');
   if fidIn ~= -1		%if filename exists, the loop ends. Otherwise input is required once more.
      break;
   end
    disp('Warning: the file does not exist in the indicated directory! Please input again.');
end

% Output file
fileOut = input('Output file: ','s');
fidOut = fopen(fileOut, 'wt');

% Input orientation data
Data = fscanf(fidIn,'%g', [6,inf]);
Data = Data';

[rowCount, colCount] = size(Data);
fclose(fidIn);

% New matrix with new orientation in Rodrigues space
RodOrients = Data; 

% Symmetry operation matrix 
global e; 
e(:,:,1)=[1 0 0; 0 1 0; 0 0 1]; e(:,:,2)=[1 0 0;0 -1 0;0 0 -1]; e(:,:,3)=[1 0 0;0 0 -1;0 1 0]; e(:,:,4)=[1 0 0;0 0 1; 0 -1 0];
e(:,:,5)=[-1 0 0;0 1 0;0 0 -1]; e(:,:,6)=[-1 0 0;0 -1 0;0 0 1]; e(:,:,7)=[-1 0 0;0 0 -1;0 -1 0]; e(:,:,8)=[-1 0 0;0 0 1;0 1 0];
e(:,:,9)=[0 1 0;-1 0 0;0 0 1]; e(:,:,10)=[0 1 0;0 0 -1;-1 0 0]; e(:,:,11)=[0 1 0;1 0 0;0 0 -1]; e(:,:,12)=[0 1 0;0 0 1;1 0 0];
e(:,:,13)=[0 -1 0;1 0 0;0 0 1]; e(:,:,14)=[0 -1 0;0 0 -1;1 0 0]; e(:,:,15)=[0 -1 0;-1 0 0;0 0 -1]; e(:,:,16)=[0 -1 0;0 0 1;-1 0 0];
e(:,:,17)=[0 0 1;0 1 0;-1 0 0]; e(:,:,18)=[0 0 1;1 0 0;0 1 0]; e(:,:,19)=[0 0 1;0 -1 0;1 0 0]; e(:,:,20)=[0 0 1;-1 0 0;0 -1 0];
e(:,:,21)=[0 0 -1;0 1 0;1 0 0]; e(:,:,22)=[0 0 -1;-1 0 0;0 1 0]; e(:,:,23)=[0 0 -1;0 -1 0;-1 0 0]; e(:,:,24)=[0 0 -1;1 0 0;0 -1 0];

for n = 1 : rowCount
    [CosMis, V] = Mis(EuToG(Data(n, 3), Data(n, 4), Data(n, 5)), InitOrient);
    % sqrt((1 - CosMis) / (1 + CosMis)): tan(theta/2)^2
    RodOrients(n, [3, 4, 5]) = sqrt((1 - CosMis) / (1 + CosMis)) * V;
end

fprintf(fidOut,'%d\t%d\t%#.4g\t%#.4g\t%#.4g\t%d\n',RodOrients'); 
fclose(fidOut);

function g = EuToG(phi1,phi,phi2)
% Transform a orientation from [phi1, phi, phi2] (in radius) to a matrix

%u=g(1,1);r=g(1,2);h=g(1,3);
%v=g(2,1);s=g(2,2);k=g(2,3);
%w=g(3,1);t=g(3,2);l=g(3,3);
g = [cos(phi1)*cos(phi2)-sin(phi1)*sin(phi2)*cos(phi), sin(phi1)*cos(phi2)+cos(phi1)*sin(phi2)*cos(phi), sin(phi2)*sin(phi);
    -cos(phi1)*sin(phi2)-sin(phi1)*cos(phi2)*cos(phi), -sin(phi1)*sin(phi2)+cos(phi1)*cos(phi2)*cos(phi),cos(phi2)*sin(phi);
    sin(phi1)*sin(phi), -cos(phi1)*sin(phi), cos(phi)];

function [CosMis, V] = Mis(A, B)
% Calculate the misorientation angle and rotation axis of A relative to B
% CosMis: cosine of the misorientation angle
% V: 1 * 3 vector, the rotation axis

global e;
% MaxTrig corresponding to the value of g1(1,1)+g1(2,2)+g1(3,3) when the
% misorientation angle is minimum
MaxTrig = 1.91; % when misorientation is 63 deg
g = A * B';
for i=1:24
	g1 = e(:,:,i) * g;
   % The bigger temp is, the smaller mis angle (0 ~ pi)
   temp = g1(1,1) + g1(2,2) + g1(3,3);
   if temp > MaxTrig
      MaxTrig = temp;
      record = i;
   end
end

g1 = e(:, :, record) * g;
CosMis = (MaxTrig - 1) / 2;
V = [g1(2,3) - g1(3,2), g1(3,1) - g1(1,3), g1(1,2) - g1(2,1)];
VNorm = norm(V); % sqrt(v1^2 + v2^2 + v3^2)
if VNorm ~= 0
   V = V / VNorm;
end