www.gusucode.com > gps卫星位置预报，MATLAB编写，有动画界面源码程序 > code5/gps_final_new/DrawSatOrbit.m

    function DrawSatOrbit(PRN,T,Toe,M0,sqrt_as,es,w,delta_i,omega0,dot_omega);

i0=delta_i+0.3*pi;
i=0;

for Mk=0:0.1:2*pi
    i=i+1;
    %计算规化时间Tk
    Tk=T-Toe;
    %计算卫星的平均角速度n
    GM=3.986005e+14;                     %地球引力常数
    n0=sqrt(GM)/sqrt_as^3;
    n=n0;
    %计算信号发射时刻的平近点角Mk, 值需要在0到2pi之间
    %Mk=mod(M0+n*Tk,2*pi);
    %计算信号发射时刻的偏近点角Ek
    TempEk=Mk;
    dEk=1;
    while(abs(dEk)>0.00001)                %迭代法   Ek1=Mk+es*sin(Ek0)
        TempEkNext=Mk+es*sin(TempEk);
        dEk=TempEk-TempEkNext;
        TempEk=TempEkNext;
    end
    Ek=TempEk;
    %计算信号发射时刻的真近点角Fk，值需要在（-pi,pi]范围内。
    Fk=2*atan(sqrt((1+es)/(1-es))*tan(Ek/2));
    %计算信号发射时刻的升交点角距Pk
    Pk=w+Fk;
    %计算摄动校正后的升角点角距uk,卫星矢径长度rk和轨道倾角ik
    uk=Pk;
    rk=sqrt_as^2*(1-es*cos(Ek));
    ik=i0;
    %计算信号发射时刻卫星在轨道平面的位置（xpk,ypk）
    xpk=rk*cos(uk);
    ypk=rk*sin(uk);
    %计算信号发射时刻的升交点赤经Omegak
    dot_omega_e=7.2921151467e-05;
    omegak=omega0+(dot_omega-dot_omega_e)*Tk-dot_omega_e*Toe;
    %计算卫星在WGS-84地心地固坐标系(ECEF)中的坐标（xk，yk，zk）
    xk(1,i)=xpk*cos(omegak)-ypk*cos(ik)*sin(omegak);
    yk(1,i)=xpk*sin(omegak)+ypk*cos(ik)*cos(omegak);
    zk(1,i)=ypk*sin(ik);
end

plot3c(xk/1000,yk/1000,zk/1000,2);
hold on;