source: proiecte/SolarSim/Matlab/keplerToCarthezian.m @ 152

function [R,V,v] = keplerToCarthezian(I,Omega,w,a,ecc,M,tol)

    if nargin < 6 || nargin > 7
        error('Incorrect number of inputs')
    elseif nargin == 6
        tol = 10^-8;
    end
    
    mu_sun = 132712440018; %km^3/s^2

    E = CalcEA(M,ecc);
    nu = nuFromM(M,ecc,tol);
    r = a * (1 - ecc*ecc)/(1+ecc*cos(nu));
    p = a * (1 - ecc*ecc);
    
    rx = r * (cos(Omega)*cos(w+nu) - sin(Omega)*cos(I)*sin(w+nu));
    ry = r * (sin(Omega)*cos(w+nu) + cos(Omega)*cos(I)*sin(w+nu));
    rz = r * sin(I) * sin(w+nu);

    v = sqrt( mu_sun * (2/r - 1/a) );
    
%    vx = -sqrt(mu_sun/p) * (cos(Omega)*(sin(w+nu)+ecc*sin(w))-sin(Omega)*cos(I)*(cos(w+nu)+ecc*cos(w)));
%    vy = -sqrt(mu_sun/p) * (sin(Omega)*(sin(w+nu)+ecc*sin(w))-cos(Omega)*cos(I)*(cos(w+nu)+ecc*cos(w)));
%    vz =  sqrt(mu_sun/p) * sin(I)*(cos(w+nu)+ecc*cos(w));

    if (I > 0)
        A = [ sin(I)*sin(Omega) -sin(I)*cos(Omega) ; rx ry ];
        B = [ -cos(I) ; -rz ];
        rez = inv(A) * B;
        a = rez(1);
        b = rez(2);

        vz = v / sqrt(1 + a*a + b*b);
        vx = a * vz;
        vy = b * vz;
    else
        vx = v / sqrt(1+(rx*rx)/(ry*ry));
        vy = - rx/ry * vx;
        vz = 0;
    end
    
    R = [rx ry rz];
    V = [vx vy vz];
    
end