aem2rtheta.m

It s verry good for you

function [r, th] = aeM2rtheta(a, e, M, varargin)
% AEM2RTHETA   converts orbital elements [a], [e] and [M] to [r] and [th].
%
%   [r, th] = AEM2RTHETA(a, e, M) converts the semi-major axis [a], the
%   eccentricity [e] and mean anomaly [M], to the polar coordinates [r] and
%   [th]. The mean anomaly [M] should be given in radians. M2RTHETA accepts 
%   matrix input of any size, as long as all matrices have the same number 
%   of elements.
%
%   AEM2RTHETA(a, e, M, tol) uses a maximum error-tolerance given by [tol].
%   The default is 1e-6. 
%
%   Algorithm: AEM2RTHETA works correct for both elliptic and hyperbolic 
%   orbits (e < 1 and e > 1, respectively). To calculate [theta] from [M], 
%   AEM2RTHETA uses a carefully selected first approximate root of Kepler's 
%   equation, followed by a Newton-Raphson iteration scheme if the first 
%   estimate is not within the limits set by [tol]. See [Seppo Mikkola, 
%   "A cubic approximation for Kepler's equation", 1987] for more details. 

%   Author: Rody P.S. Oldenhuis
%   Delft University of Technology
%   E-mail: [e-mail]
%   Last edited 29/Jan/2009

    % errortraps
    error(nargchk(3, 4, nargin));
    if (numel(M) ~= numel(a)) || (numel(M) ~= numel(e)) || (numel(e) ~= numel(a))
        error('Number of elements of all input arguments must be equal');
    end 
        
    % parse input
    tol = 1e-6;
    if (nargin >= 4)
        tol = varargin{1};
        if isempty(tol), tol = 1e-6; end        
    end    
    
    % transform matrices into vectors
    sza = size(a);
    M   = M(:);
    a   = a(:);
    e   = e(:);
    
    % initial values
    r  = zeros(numel(a), 1);
    th = r; 
    
    % seperate elliptical from hyperbolic orbits
    ell = (e < 1);
    hyp = (e > 1); 
    
    % Elliptic orbits
    if any(ell)  
        
        % extract appropriate values
        Me = M(ell); ee = e(ell);
        
        % put M in interval (-pi < M < +pi)
        if any(Me > pi) || any(Me < -pi)
            Me  = mod(Me, 2*pi);
            Mep = Me > pi;
            if any(Mep)
                Me(Mep) = Me(Mep) - 2*pi;
            end
        end
        
        % apprioximate root from Mikkola's paper             
        gamma = (4*ee + 1/2);
        alpha = (1 - ee)./gamma;
        beta  = Me./(2*gamma);
        sgnb  = sign(beta);
        z  = (beta + sgnb.*sqrt(beta.^2 + alpha.^3));
        z  = sign(z).*abs(z).^(1/3);           
        z(z==0) = 1e-100;  % prevent 1-over-0 warning
        s  = z - alpha./z;
        s  = s - 0.078*s.^5./(1 + ee);
        Ee = Me + ee.*(3*s - 4*s.^3);
        
        % perform at least one Newton-Raphson iteration, additional ones if 
        % required by [tol]       
        Ep = inf;
        while any( abs(Ee - Ep) > tol)
            Ep = Ee;           
            Ee = Ee + (Me + ee.*sin(Ee) - Ee) ./ (1 - ee.*cos(Ee));
        end    
        
        % correct for M = 0 case
        Ee(Me == 0) = 0;
                
        % calculate theta        
        th(ell) = 2*atan2(sqrt(1 + ee).*sin(Ee/2), sqrt(1 - ee).*cos(Ee/2));      
                
    end
    
    % Hyperbolic orbits
    if any(hyp)
        
        % extract relevant values
        eh = e(hyp); Mh = M(hyp);
        
        % apprioximate root from Mikkola's paper
        gamma = (4*eh + 1/2);
        alpha = (eh - 1)./gamma;
        beta  = Mh./(2*gamma);
        sgnb  = sign(beta);        
        z  = (beta + sgnb.*sqrt(beta.^2 + alpha.^3));
        z  = sign(z).*abs(z).^(1/3);
        s  = z - alpha./z;  
        s2 = s.^2;
        s  = s + 0.071*s.^5 ./ ((1 + 0.45.*s2).*(1 + 4*s2).*eh);   
        Eh = 3*log(s + sqrt(1 + s.^2));          
        
        % perform at least one Newton-Raphson iteration, additional ones if 
        % required by [tol]  
        Ep = inf;
        while any( abs(Eh - Ep) > tol)
            Ep = Eh;
            Eh = Eh + (Eh - eh.*sinh(Eh) + Mh) ./ (eh.*cosh(Eh) - 1);
        end
        
        % calculate theta 
        th(hyp) = 2*atan2(sqrt(eh + 1).*sinh(Eh/2), sqrt(eh - 1).*cosh(Eh/2));    
        
    end
    
    % radius vectors
    r = a.*(1 - e.^2) ./ (1 + e.*cos(th));

    % transform back into matrix
    th = reshape(th, sza);
    r  = reshape(r, sza);
    
end