function [R, varargout] = euler_123_SCARF(alpha, beta, gamma)
% R = euler(alpha, beta, gamma)
% Rotation matrix R, "123" type.
% Output rotation matrix R is defined as "left-handed positive sign
% convention XYZ rotation" in http://en.wikipedia.org/wiki/Euler_angles
% This is equivalent to eq. (2) of Olsen et al, The Swarm Satellite Constellation
% Application and Research Facility (SCARF) and Swarm data products,
% EPS 65, 2013
%
% R = euler(alpha, beta, gamma);
% [R, R_a, R_b, R_g] = euler(alpha, beta, gamma);
% R_a = d R / d alpha, etc.
% Input: Euler angles alpha, beta and gamma (in radians)
%
% Nils Olsen, DTU Space, January 2014.

s_a = sin(alpha); s_b = sin(beta); s_g = sin(gamma);
c_a = cos(alpha); c_b = cos(beta); c_g = cos(gamma);

R(:,1,1) = c_b.*c_g;
R(:,2,1) = c_b.*s_g;
R(:,3,1) = -s_b;
R(:,1,2) = s_a.*s_b.*c_g-c_a.*s_g;
R(:,2,2) = s_a.*s_b.*s_g+c_a.*c_g;
R(:,3,2) = s_a.*c_b;
R(:,1,3) = c_a.*s_b.*c_g+s_a.*s_g;
R(:,2,3) = c_a.*s_b.*s_g-s_a.*c_g;
R(:,3,3) = c_a.*c_b;

% R = [c_b.*c_g       s_a.*s_b.*c_g-c_a.*s_g  c_a.*s_b.*c_g+s_a.*s_g;
%      c_b.*s_g       s_a.*s_b.*s_g+c_a.*c_g  c_a.*s_b.*s_g-s_a.*c_g  ;
%      -s_b           s_a.*c_b                c_a.*c_b   ];

if size(R,1) == 1
    R = squeeze(R);
end

if (nargout > 1)
    o = zeros(size(alpha));
    % derivative of R vrt. alpha
    R_a(:,1,1) = o;
    R_a(:,2,1) = o;
    R_a(:,3,1) = o;
    R_a(:,1,2) = c_a.*s_b.*c_g+s_a.*s_g;
    R_a(:,2,2) = c_a.*s_b.*s_g-s_a.*c_g;
    R_a(:,3,2) = +c_a.*c_b;
    R_a(:,1,3) = -s_a.*s_b.*c_g+c_a.*s_g;
    R_a(:,2,3) = -s_a.*s_b.*s_g-c_a.*c_g;
    R_a(:,3,3) = -s_a.*c_b;
%     R_a = [0                        0                        0;
%         c_a.*s_b.*c_g+s_a.*s_g   c_a.*s_b.*s_g-s_a.*c_g   +c_a.*c_b;
%         -s_a.*s_b.*c_g+c_a.*s_g  -s_a.*s_b.*s_g-c_a.*c_g   -s_a.*c_b  ]';
    
    % derivative of R vrt. beta
    R_b(:,1,1) = -s_b.*c_g;
    R_b(:,2,1) = -s_b.*s_g;
    R_b(:,3,1) = -c_b;
    R_b(:,1,2) = s_a.*c_b.*c_g;
    R_b(:,2,2) = s_a.*c_b.*s_g;
    R_b(:,3,2) = -s_a.*s_b;
    R_b(:,1,3) = c_a.*c_b.*c_g;
    R_b(:,2,3) = c_a.*c_b.*s_g;
    R_b(:,3,3) = -c_a.*s_b;
%     R_b = [-s_b.*c_g             -s_b.*s_g           -c_b;
%         s_a.*c_b.*c_g            s_a.*c_b.*s_g     -s_a.*s_b;
%         c_a.*c_b.*c_g            c_a.*c_b.*s_g     -c_a.*s_b   ]';
    
    % derivative of R vrt. gamma
    R_g(:,1,1) = -c_b.*s_g;
    R_g(:,2,1) =  +c_b.*c_g;
    R_g(:,3,1) = o;
    R_g(:,1,2) = -s_a.*s_b.*s_g-c_a.*c_g;
    R_g(:,2,2) = +s_a.*s_b.*c_g-c_a.*s_g;
    R_g(:,3,2) = o;
    R_g(:,1,3) = -c_a.*s_b.*s_g+s_a.*c_g;
    R_g(:,2,3) = c_a.*s_b.*c_g+s_a.*s_g;
    R_g(:,3,3) = o;
%     R_g = [-c_b.*s_g                  +c_b.*c_g                 0;
%         -s_a.*s_b.*s_g-c_a.*c_g    +s_a.*s_b.*c_g-c_a.*s_g   0;
%         -c_a.*s_b.*s_g+s_a.*c_g    c_a.*s_b.*c_g+s_a.*s_g    0   ]';

if size(R_a,1) == 1
    R_a = squeeze(R_a);
    R_b = squeeze(R_b);
    R_g = squeeze(R_g);
end

    varargout(1) = {R_a};
    varargout(2) = {R_b};
    varargout(3) = {R_g};
end

