计算导航卫星位置

% 读取文件并跳过文件头
clear;
fp=fopen（‘shidun0707.13N‘‘r‘）;
fseek（fp574‘bof‘）;

t = [137710300];
disp（[‘观测时间：‘num2str（t（1））‘年‘num2str（t（2））‘月‘num2str（t（3））‘日‘ ...
        num2str（t（4））‘时‘num2str（t（5））‘分‘num2str（t（6））‘秒‘]）;
for k =1:27
    toetmp = str2num（fgetl（fp））;
    PRN = toetmp（1）;
    toe = toetmp（2:7）;
    orbitInfo = zeros（74）;
	% 读取导航电文
    for row = 1:7
        if row<7
            orbitInfo（row:） = str2num（fgetl（fp））;
        else
            orbitInfo（row1:2） = str2num（fgetl（fp））;
        end
    end
    
    dt = 3600*（t（4）-toe（4））+60*（t（5）-toe（5））+（t（6）-toe（6））;
    if abs（dt）<=3600
        % 平均角速度
        GM = 3.986005*10^14;
        sqrtA = orbitInfo（24）;
        n0 = GM^0.5/sqrtA^3;
        deltan = orbitInfo（13）;
        n = n0 + deltan;
        % 平近点角
        M0 = orbitInfo（14）;
        M = M0 + n*dt;
        % 偏近点角
        E = M;
        e = orbitInfo（22）;
        while（（M+e*sin（E）-E）>10^（-6））
            E = M+e*sin（E）;
        end
        % 真近点角
        f = atan（（1-e^2）^0.5*sin（E）/（cos（E）-e））;
        % 升交点角距
        u_prime = orbitInfo（4