function [r,v,rx,ry,rz,vx,vy,vz] = orbelm_orbelm2xv(ao,ecc,inc,Omega,omega,MA,M);

%ao = semi-major axis
%inc = inclination
%MA = mean anomaly 
%ecc = eccentricity
%Omega = right ascension of ascending node
%omega = argument of perigee
%E = eccentric anamoly
%M = Body mass

%use the sun if M not input
if ~exist('M','var');
    M = [];
end;
if isempty(M);
    M = 1.98892e30; %kg
    %M = 5.6851e26; %kg (Saturn)
end;

%solve for eccentric anamoly
thresh = 1e-15;
maxcnt = 1000;
for j=1:numel(ao.*ecc.*inc.*Omega.*MA.*M);
    if numel(MA) > 1;
        Eo = MA(j);
        MA_temp = MA(j);
    else;
        Eo = MA(1);
        MA_temp = MA(1);
    end;
    if numel(ecc) > 1;
        ecc_temp = ecc(j);
    else;
        ecc_temp = ecc;
    end;
    E_temp = Eo - (Eo-ecc_temp.*sin(Eo)-MA_temp)./(1-ecc_temp.*cos(Eo));
    cnt = 1;
    while ((abs(E_temp - Eo) > thresh) & (cnt < maxcnt));
        Eo = E_temp;
        E_temp = Eo - (Eo-ecc_temp.*sin(Eo)-MA_temp)./(1-ecc_temp.*cos(Eo));
        cnt = cnt+1;
    end;
    if cnt==maxcnt; 
        disp(['Reached Max Count For Index ' num2str(j)]);
    elseif cnt==1;        
        disp(['Performed No Iterations for Index ' num2str(j)]);
    end;
    E(j) = E_temp;   
end;

%compute vectors P and Q in the PQW frame
Px = cos(omega).*cos(Omega) - sin(omega).*cos(inc).*sin(Omega);
Py = cos(omega).*sin(Omega) + sin(omega).*cos(inc).*cos(Omega);
Pz = sin(omega).*sin(inc);
Qx = -sin(omega).*cos(Omega) - cos(omega).*cos(inc).*sin(Omega);
Qy = -sin(omega).*sin(Omega) + cos(omega).*cos(inc).*cos(Omega);
Qz =  sin(inc).*cos(omega);

%definet the constant mu in Keplers third law
G = 6.67300e-11; %m^3./hs./s^2
mu = G.*M;

%calculate the cartesian position r and velocity v
rx = ao.*(cos(E)-ecc).*Px + ao.*sqrt(1-ecc.^2).*sin(E).*Qx;
ry = ao.*(cos(E)-ecc).*Py + ao.*sqrt(1-ecc.^2).*sin(E).*Qy;
rz = ao.*(cos(E)-ecc).*Pz + ao.*sqrt(1-ecc.^2).*sin(E).*Qz;
vx = (1./(1-ecc.*cos(E))).*sqrt(mu./ao).*(-ecc.*sin(E).*Px + sqrt(1-ecc.^2).*cos(E).*Qx);
vy = (1./(1-ecc.*cos(E))).*sqrt(mu./ao).*(-ecc.*sin(E).*Py + sqrt(1-ecc.^2).*cos(E).*Qy);
vz = (1./(1-ecc.*cos(E))).*sqrt(mu./ao).*(-ecc.*sin(E).*Pz + sqrt(1-ecc.^2).*cos(E).*Qz);

r = zeros(numel(rx),3);
r(:,1) = rx;
r(:,2) = ry;
r(:,3) = rz;
v = zeros(size(r));
v(:,1) = vx;
v(:,2) = vy;
v(:,3) = vz;

