[ax,ay,az,LONG,LAT,DEC]=rotation(x,y,z,A,D,W,sense,flatt)

function [ax,ay,az,LONG,LAT,DEC]=rotation(x,y,z,A,D,W,sense,flatt)

% function [ax,ay,az,LONG,LAT,DEC]=rotation(x,y,z,A,D,W,sense,flatt)
%
% rotational parameters of a planet
% algorithms by Montenbruck/Pfleger: "Astronomy with the Personal Computer"
%
% 05.05.04 M.Penzkofer

% constants
rho = pi/180;


% check input parameters
if (isequal(sense,'DIRECT') == 0 & isequal(sense,'RETROGRADE') == 0)
   disp('Enter DIRECT or RETROGRADE for parameter SENSE!');
   return
end

% compute unit vectors E(*,1)
E = gaussvec(90.0+A,90.0-D,W);

% copy of unit vector in direction of the rotation axis
ax = E(1,3);  ay = E(2,3);  az = E(3,3);

% planetocentric longitude and latitude
SX = - ( E(1,1)*x + E(2,1)*y + E(3,1)*z );
SY = - ( E(1,2)*x + E(2,2)*y + E(3,2)*z );
SZ = - ( E(1,3)*x + E(2,3)*y + E(3,3)*z );
T  = SZ / sqrt(SX*SX+SY*SY);
DEC  = atan(T) / rho;
LONG = atn2(SY,SX);

% planetographic longitude and latitude
if (isequal(sense,'DIRECT') == 1)
   LONG=-LONG;
end
if (LONG < 0.0)
   LONG=LONG+360.0;
end
LAT = atan( T / ((1.0-flatt)*(1.0-flatt)) ) / rho;