function [long_out,lat_out]=Euler_Rotation(long_pole,lat_pole,rotang,long_in,lat_in);
% function to rotation points around an axis by more than one angle
% The rotation axis is described by its latitude and longitude of intersection
% with the sphere
%Usage:
%   [long_out,lat_out]=Euler_Rotation(long_pole,lat_pole,rotang,long_in,lat_in)
%
% where: long_pole and lat_pole define the rotation axis
%        rotangle is vector of angles to rotate
%        long_in and lat_in are the coordinates of points to rotate
%        long_out and lat_out are the coordinates of the rotated points



% internally we use a colatitude
  colat_pole=90-lat_pole;
  colat_in=90-lat_in;
  
% make sure colat and long are column vectors
  colat_in=colat_in(:);
  long_in=long_in(:);
  
  
 %convert lat's and long's into x-y-z coordinates:
 scp=sin(colat_in*pi/180);ccp=cos(colat_in*pi/180);
 slp=sin(long_in*pi/180);clp=cos(long_in*pi/180);
 xyz=[scp.*clp scp.*slp ccp];
 [np,nc]=size(xyz); 
  
 % calculate sin and cos of pole colat. and long. and the rotation angles:
 sc=sin(colat_pole*pi/180);cc=cos(colat_pole*pi/180);
 sl=sin(long_pole*pi/180);cl=cos(long_pole*pi/180);
 sr=sin(rotang*pi/180);cr=cos(rotang*pi/180);
 nr=length(rotang);
 
 % matrix for coordinate transformation that moves the rotation axis to the z-axis 
 R1=[cc 0 -sc;0 1 0;sc 0 cc]*[cl sl 0;-sl cl 0 ;0 0 1];
 
 long_out=zeros(np,nr);
 lat_out=long_out;
 
 % loop to do the rotation by different angles:
 for i=1:nr,
     %  matrix for rotation about z
     R2=[cr(i) -sr(i) 0;sr(i) cr(i) 0;0 0 1];
 
     %here it is the grand transformation:
     xyzp=(R1'*R2*R1*xyz')';
 
     % now convert xyz back to latitudes and longitudes (force to be real to
     % avoid problems with infinities at poles
 
     colat_out=real(acos(xyzp(:,3))); % still in radians
     long_out(:,i)=(180/pi*asin(xyzp(:,2)./sin(colat_out)));
     lat_out(:,i)=90-180/pi*colat_out;
 end
 
 
 