shg;R0=1;												%以地球半徑為一個單位
a=12*R0;b=9*R0;T0=2*pi;								%T0是軌道周期
T=5*T0;dt=pi/100;t=[0:dt:T]';
f=sqrt(a^2-b^2);										%地球與另一焦點的距離
th=12.5*pi/180;										%衛星軌道與x-y平面的傾角
E=exp(-t/20);											%軌道收縮率
x=E.*(a*cos(t)-f);y=E.*(b*cos(th)*sin(t));z=E.*(b*sin(th)*sin(t));
plot3(x,y,z,'g')										%畫全程軌線
[X,Y,Z]=sphere(30);X=R0*X;Y=R0*Y;Z=R0*Z;		%獲得單位球坐標
grid on,hold on,surf(X,Y,Z),shading  interp	%畫地球
x1=-18*R0;x2=6*R0;y1=-12*R0;y2=12*R0;z1=-6*R0;z2=6*R0;
axis([x1 x2 y1 y2 z1 z2])							%確定坐標范圍
view([117 37]),comet3(x,y,z,0.02),hold off		%設視角、畫運動軌線