function [x, y, z, xdot, ydot, zdot] = kepler2cart(varargin)
%KEPLER2CART            Converts Kepler elements to Cartesian coordinates
%
%   [x, y, z, xdot, ydot, zdot] = KEPLER2CART(a, e, i, o, O, M, mu) or 
%   equivalently, KEPLER2CART( [Keplerian elements], mu) will convert the 
%   given orbital elements into Cartesian coordinates. The input value [mu] 
%   is a scalar value -- the standard gravitational parameter of the
%   central body in question. This is the same as calling KEPLER2CART(a, e, 
%   i, o, O, theta, mu, 'M').
%
%   In case the Keplerian elements are given as a single Matrix, the 
%   Keplerian coordinates are taken columnwise from it. In the more general 
%   case of providing each element seperately, the elements may be matrices 
%   of any dimension, as long as the dimensions are equal for each input. 
%
%   Note that the Mean anomaly [M] is used for the default conversion. 
%   Calling KEPLER2CART(a, e, i, o, O, theta, mu, 'theta') (or its
%   matrix form) will use the true nomaly [theta] directly. Note that this 
%   is NOT the default mode of operation, and the last argument (the string 
%   'theta') must be included explicitly, in both the matrix form of
%   individual element form. 

%   Author: Rody P.S. Oldenhuis
%   Delft University of Technology
%   E-mail: oldenhuis@dds.nl
%   Last edited 14/Nov/2008

%%  initialize

    % default errortrap
    error(nargchk(2, 8, nargin));
    
    % parse input
    thorM = 'M';
    if (nargin >= 2) && (nargin <= 3)
        Kepler = varargin{1};
        mu     = varargin{2};        
        thorM  = 'M';
        a      = Kepler(:, 1);        e      = Kepler(:, 2);
        i      = Kepler(:, 3);        o      = Kepler(:, 4);
        O      = Kepler(:, 5);        Mth    = Kepler(:, 6);
    end
    if (nargin == 3)        
        thorM  = varargin{3};         
    end
    if (nargin >= 7)        
        a      = varargin{1};        e      = varargin{2};
        i      = varargin{3};        o      = varargin{4};
        O      = varargin{5};        Mth    = varargin{6};
        mu     = varargin{7};        thorM  = 'M';
    end
    if (nargin == 8)         
        thorM  = varargin{8};        
    end
    
    % select theta or M
    if isempty(thorM), thorM = 'M'; end
    th = strcmpi(thorM, 'theta');

    % make it work for [n]-dimensional input
    sizea = size(a);  
    a   = a(:);    e   = e(:);
    i   = i(:);    o   = o(:);
    O   = O(:);    Mth = Mth(:);
    

%%  start conversion    
    
    % compute constants required for transformation
    l1 =  cos(O).*cos(o) - sin(O).*sin(o).*cos(i);
    l2 = -cos(O).*sin(o) - sin(O).*cos(o).*cos(i);
    m1 =  sin(O).*cos(o) + cos(O).*sin(o).*cos(i);
    m2 = -sin(O).*sin(o) + cos(O).*cos(o).*cos(i);
    n1 =  sin(o).*sin(i);
    n2 =  cos(o).*sin(i);
    H  =  sqrt( mu.*a.*(1 - e.^2) );    
    if (th)
        theta = Mth;
        r     = a.*(1-e.^2) ./ (1 + e.*cos(theta));
    else
        [r, theta] = aeM2rtheta(a, e, Mth);
    end 
         
    % transform
    xi   = r.*cos(theta);
    eta  = r.*sin(theta);
    one  = [l1; m1; n1] .* [xi; xi; xi];
    two  = [l2; m2; n2] .* [eta; eta; eta];    
    rvec = one + two;
    H    = repmat(H, 3, 1);
    Vvec = mu./H .* [-l1.*sin(theta) + l2.*(e+cos(theta)) ;
                     -m1.*sin(theta) + m2.*(e+cos(theta)) ;
                     -n1.*sin(theta) + n2.*(e+cos(theta))];

    % assign values and reshape
    rows = size(rvec, 1) / 3;    
    x    = reshape(rvec(1:rows)       , sizea);
    y    = reshape(rvec(rows+1:2*rows), sizea);
    z    = reshape(rvec(2*rows+1:end) , sizea);
    xdot = reshape(Vvec(1:rows)       , sizea);
    ydot = reshape(Vvec(rows+1:2*rows), sizea);
    zdot = reshape(Vvec(2*rows+1:end) , sizea);
    
end