function [f,Orbit_xyz,r,Tt]=Earthorbit(a,e,Tp,Tt,Mui,Ta)while Tt>Tp                                       %%%將當前時刻時間轉換到衛星運行周期之內    Tt=Tt-Tp;end%%%用迭代法求解開普勒方程，求解偏近點角%%%平近點角MM=2*pi/Tp*(Tt-Ta);E=MM;times=0;p=10;while p>1e-20    temp=E;    E=temp-(temp-e*sin(temp)-MM)/(1-e*cos(temp));    p=abs(E-temp);    times=times+1;    if times>10000        p=1e-20;    endend%%%%根據偏近點角E求解真近點角ff=2*atan(sqrt((1+e)/(1-e))*tan(E/2));%%%%求出衛星距地心的距離r=(a*(1-e^2))/(1+e*cos(f));Orbit_xyz=[r*cos(f) r*sin(f) 0]';