function [x,y,z] = drawOrbit(a, e, i, Omega, omega)
%% 返回空间轨道椭圆的点集（X(1)为近地点坐标） 
% a                 半长轴
% e     [0,1)       偏心率
% i     [0,pi]      轨道倾角
% Omega [0,2*pi)    升交点赤经
% omega [0,2*pi)    近地点张角（幅角）

global precision 
theta = -pi:precision:pi;
rho = a*(1-e^2)./(1-e*cos(theta));
[x,y] = pol2cart(theta,rho);
z = zeros(1,length(x));

% 变换矩阵（不是基元旋转矩阵）
M1 = makehgtform('axisrotate',[1,0,0],i);
m1 = M1(1:3,1:3);

M2 = makehgtform('axisrotate',[0,0,-1],Omega);
m2 = M2(1:3,1:3);

n = [0,0,-1]*m1*m2;
M3 = makehgtform('axisrotate',-n,omega);
m3 = M3(1:3,1:3);

for j = 1:length(x)
    tmp = [x(j),y(j),0]*m1*m2*m3;
    x(j) = tmp(1);
    y(j) = tmp(2);
    z(j) = tmp(3);
end
    
plot3([x x(1)],[y y(1)],[z z(1)],'Color',[0.4660 0.6740 0.1880]);

end