function variants

% first we have to set up some global matrices

% Crystal 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];

% note how Matlab automatically sets up variables for you
%  although square brackets around the input are essential here !!!
angles_in = input('Enter three Bunge Euler angles, in degrees [phi1, PHI, phi2]:  ');
%  angles_in
%  if you want the program to echo your input, un-comment the line above
%  (remove the "%" )

r_angles_in = angles_in * pi / 180.;
%  again, Matlab automatically figures out that you want r_angles_in to be
%  an array (vector) just like angles_in
fprintf(' angles in radians = %.4g  %.4g  %.4g \n',r_angles_in);

disp('  now we feed the 3 angles to the subroutine "EuToG" that converts angles');
disp('  to a matrix; in this case we have to spell out each argument');
Orient_matrix = EuToG(r_angles_in(1),r_angles_in(2),r_angles_in(3));
disp('  now we echo the matrix just calculated');
Orient_matrix

disp('The next step is to multiply the orientation (matrix) by a symmetry operator');
disp(' except we put this in a loop so that we can obtain a series of results');

% set up an Output file
fileOut = input('Enter a name for the Output file: ','s');
fidOut = fopen(fileOut, 'wt');
    fprintf(fidOut,' loop_number  phi1   PHI   phi2  (degrees) \n');


for i = 1:24
    disp('  next we pre-multiply by a symmetry (crystal) matrix ');
    new_orient = e(:,:,i) * Orient_matrix;
    disp('  now we echo the new matrix just calculated');
    fprintf(' for loop number: %d  \n',i);
    new_orient
    disp('and convert it back to Euler angles');
% note how we set up the function call
%  the output variables are in the square brackets on the LHS
%  and the inputs to the fucntion Euler are in the parentheses
    [d1,d2,d3] = Euler(new_orient,'B');
    fprintf(' angles in degrees = %.4g  %.4g  %.4g \n',d1,d2,d3);
    fprintf(fidOut,' %d   %10.4f  %10.4f  %10.4f \n',i,d1,d2,d3);
   
end

fclose(fidOut);

end


function g = EuToG(phi1,phi,phi2)
% Transform a orientation from [phi1, phi, phi2] (in radians) 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)];

end

function [d1,d2,d3] = Euler(a,nomen)
%  the variables d1,d2 & d3 are the outputs of the function Euler
%  which is copied from the fortran subroutine of the same name
%  a, the matrix and nomen, defining the type of Euler angle, are inputs
%  NB  conversion from fortran required adaptation from having the matrix 
%  in transposed form (as used by Canova in LApp!!!)

% c                          Last revision 20nov90 UFK
% c      common a(3,3),grvol(1152),epsga(5),ist1,ist2,sqr3,sqrh,ident(3,3)
% c  SPECIAL VERSION WITHOUT COMMON BLOCK
% 	real a(3,3),d1,d2,d3,th,sth,cth,sph,cph,sps,cps,ps,ph,dth,dph,dps
% 	real pi,rad
%       char nomen
%       save kor
% 	PI=3.14159265
% 	rad=57.29578
% c
% c *** this subroutine calculates the euler angles associated with the
% c     transformation matrix a(i,j) if iopt=1 and viceversa if iopt=2
% c ***  Note that a is sample (rows) in terms of crystal (columns);
% c      -- opposite of standard g (e.g.Bunge) - this is Canova's
% c ***  Note that in this version, the Euler angles are defined symmetrically:
% c      so that interchanging phi and psi means transposing a.
% c      ("Kocks" nomen: defined going from +X to +Y in both COD and SOD)
% c ***  However, other angle conventions are translated, according to
% c  nomen="K" - Kocks (as internally) -- also sometimes "N"...
% c        "R" - Roe          (Psi=psi, Phi=180-phi)
% c        "B" - Bunge        (phi1=90+psi, Phi=Theta, phi2=90-phi)
% c        any other - Canova (Theta first, phiC=90+phi, omega=90-psi)
% c ***   Note: only in symm. notation does a point with all Euler angles
% c             between 0 and 90 deg appear in the +x,+y quadrant!
% c       If you want to see an individual point in this quadrant and are:
% c       using Roe,    the third angle must be between 90 and 180;
% c             Bunge,      first                                 ;
% c             Canova,     second                                .
% c ***  Input and output Euler angles d1,d2,d3 in degrees

% fprintf('nomen = %s \n',nomen);
% disp('echoing the input matrix ...');
% 
% a

   if abs(a(3,3)) < 0.999
       %  treatment depends on how close the 2nd angle is to zero
       %  if not too close, then use standard ATAN2 method
      th = acos(a(3,3));
      sth = sin(th);

	if abs(a(2,3)/sth) < 1.e-35 && abs(a(1,3)/sth) < 1.e-35
        disp('Euler: low values for 2,3 and 1,2');
		ph = pi/4.;
	else
      ph = atan2(a(2,3)/sth,a(1,3)/sth);
%      d3 = atan2(a(1,3)/sth,a(2,3)/sth);
    end

	if abs(a(3,2)/sth) < 1e-35 && abs(a(3,1)/sth) < 1.e-35
        disp('Euler: low values for 3,2 and 3,1');
		ps=pi/4.;
    else
      ps = atan2(a(3,2)/sth,a(3,1)/sth);
%      d1 = atan2(a(3,1)/sth,-1.*a(3,2)/sth);
    end
%   d2 = th;
%         fprintf('  angles in degrees = %10.4f  %10.4f  %10.4f \n',ps*180/pi,th*180/pi,ph*180/pi);
%         fprintf('  Bunge angles ? in degrees = %10.4f  %10.4f  %10.4f \n',d1*180/pi,d2*180/pi,d3*180/pi);
%  this line was a debugging line   
    
% ccccc   if it bombs out here, both a-components in one arg. are zero
% c       (this should not be possible, but has happened, probably fixed)
% c  ADR: i 01 - above is the fix!!
% c
   else

 if abs(a(2,1)/sth) < 1.e-35 && abs(a(1,1)/sth) < 1.e-35
		ps = pi/4.;
	else
		ps = 0.5*atan2(a(2,1),-a(1,1));
    end

      ph=-ps;
% c       The above still have the problem that they give too many
% c       equivalents of the same grain in DIOROUT and density file. Therefore:
%       if kerr == 1 && kor ~= ior 
%         disp('WARNING: orientation has Theta < 1 deg.: sometimes problems');
%         kor=ior;
%       end
   
   end
   
   dth = th * 180/pi;
   dph = ph * 180/pi;
   dps = ps * 180/pi;
   d1 = dps;
   d2 = dth;
    if nomen == 'K' || nomen == 'N'
        d3 = dph;
      elseif nomen == 'R'
        d3 = 180. - dph;
      elseif nomen == 'B'
        d1 = dps + 90.;
        d3 = 90. - dph;
      else
        d1 = dth;
        d2 = dph + 90.;
        d3 = 90. - dps;
      end

      if(d1 > 360.) 
          d1 = d1-360.;
      end
      if(d3 > 360.) 
          d3 = d3-360.;
      end
      if(d1 < 0.) 
          d1 = d1+360.;
      end
      if(d3 < 0.) 
          d3 = d3+360.;
      end

end