（三）基于MATLAB计算卫星位置

（二）基于MATLAB计算卫星位置

参数由读取导航文件（.n）得到。

GM=3.986005e+14;   %万有引力常数
OMEGAedot=7.2921151467e-5;   %地球自转角速度 （rad）
  %(1)计算规划时间 
    if  tk >302400
        tk=tk-604800;
    end
    if  tk <-302400
        tk=tk+604800;
    end;
    %(2)计算卫星平均角速度并校正
    n0=sqrt(GM)/Nav(j).sqrtA^3;
    n=n0+ Nav(j).dn;
    %(3)计算平近点角Mk
    Mk=Nav(j).M0+n*tk;  
    %(4)计算偏近点角Ek
    %EK在电离层校正时会用到
    Ek=Mk;
    while 1
      Ek_=(Mk+Nav(j).ec*sin(Ek));
      if (abs(Ek_-Ek)<10^(-12))
          Ek=Ek_;
          break;
      else
         Ek=Ek_;
      end;
    end; 
    if Ek < 0
       Ek = Ek + 2*pi;
    end
    %(5)计算真近点角
    Vk = atan(sqrt(1-Nav(j).ec^2)*sin(Ek)/(cos(Ek)-Nav(j).ec));
    if (cos(Ek)-Nav(j).ec)/(1-Nav(j).ec*cos(Ek)) < 0
        Vk = Vk - pi;
    end
    %(6)计算升交点角距
    u=Vk+Nav(j).omega;
    %(7)摄动校正项
    delta_uk = Nav(j).cus*sin(2*u) + Nav(j).cuc*cos(2*u);
    delta_rk = Nav(j).crs*sin(2*u) + Nav(j).crc*cos(2*u);
    delta_ik = Nav(j).cis*sin(2*u) + Nav(j).cic*cos(2*u);
    %(8)代入校正项
    uk = u + delta_uk;
    rk = Nav(j).sqrtA.^2*(1-Nav(j).ec*cos(Ek)) + delta_rk;
    ik = Nav(j).i0 + Nav(j).idt*tk + delta_ik;
    %(9)卫星在轨道平面位置
    xk = rk*cos(uk);
    yk = rk*sin(uk);
    %(10)升交点赤经
    OMEGAk = Nav(j).OMEGA0 + (Nav(j).OMEGAdot - OMEGAedot)*tk - OMEGAedot*Nav(j).toe;
    %(11)坐标系变换至WGS-84坐标系
    Xk = xk*cos(OMEGAk) - yk*cos(ik)*sin(OMEGAk);
    Yk= xk*sin(OMEGAk) + yk*cos(ik)*cos(OMEGAk);
    Zk = yk*sin(ik);
    %(12)地球自转修正
    X=Xk*cos(OMEGAedot*Dt)+Yk*sin(OMEGAedot*Dt);
    Y=-Xk*sin(OMEGAedot*Dt)+Yk*cos(OMEGAedot*Dt);
    Z=Zk;
    SatPos=[X,Y,Z];