Matlab根据广播星历表计算卫星坐标

    clear,clc;
    
    %读取卫星星历文件
    filename='20040130.04N.txt';
    fid=fopen(filename);
    
    %按行读取星历文件的每一个元素，从PRN7号卫星开始
    data=textscan(fid,'%f %f %f %f %f %f %f %f %f
    %f\n','HeaderLines',11);
    
    %将读取的数组转换成数值型矩阵
    B1=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B2=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B3=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B4=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B5=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B6=cell2mat(data);
    data=textscan(fid,'%f %f %f %f\n','HeaderLines',1);
    B7=cell2mat(data);
    data=textscan(fid,'%f\n','HeaderLines',1);
    B8=cell2mat(data);
    
    %定义常量GM,Ve
    GM=3.986004e14;
    Ve=7.2921151467e-5;
    
    %求出轨道长半轴a
    sqrta=B3(4);
    a=sqrta^2;
    e=B3(2);
    
    %求出平均角速度n0
    n0=sqrt(GM/(a^3));
    XYZ=[];
    
    %计算从2004年1月30日16：00开始，每隔1分钟的卫星坐标
    for i=0:1:20
    %计算当前时刻的儒略日
    JulianDay=fix(365.25*(2004-1))+fix(30.6001*(1+13))+30+(16+i/60+0/3600)/24+1720981.5;
    %将儒略日转换成GPS时
    	t2=(JulianDay-2444244.5)/7-fix((JulianDay-2444244.5)/7);
    t=t2*604800;
    
    %计算卫星钟差改正后的当前时刻的GPS时
    	t=t-(B1(8)+B1(9)*(t-B4(1))+B1(10)*(t-B4(1))^2);
    	
    %计算参考时刻的儒略日
    	JulianDayoe=fix(365.25*(2004-1))+fix(30.6001*(1+13))+30+(8+0/60+0/3600)/24+1720981.5;
    
     %将儒略日转换成GPS时
       	t2oe=(JulianDayoe-2444244.5)/7-fix((JulianDayoe-2444244.5)/7);
    toe=t2oe*604800;
    
    %计算时间差tk
    tk=t-toe;
    
    %计算改正平角速度
    n=n0+B2(3);
    
    %计算平近点角Ek
    Mk=B2(4)+n*tk;
    
    %迭代计算，求出偏近点角Ek
    Ek=Mk;
    Ek0=inf;
    
    while
    		abs(Ek0-Ek)>1e-100
        Ek0=Ek;
        Ek=Mk+e*sin(Ek0);
    end
    
    %计算真近点角
    	Sinfk=(sqrt(1-e^2)*sin(Ek))/(1-e*cos(Ek));
    Cosfk=(cos(Ek)-e)/(1-e*cos(Ek));
    	fk=atan2(Sinfk,Cosfk);
    
    %计算升交角距
    phik=fk+B5(3);
    
    %从卫星星历文件中读取各个方向的振幅Crs,Cuc,Cus,Cic,Cis,Crc
    Crs=B2(2);
    Cuc=B3(1);
    Cus=B3(3);
    Cic=B4(2);
    Cis=B4(4);
    Crc=B5(2);
    
    %计算卫星轨道摄动项改正数
    deltauk=Cus*sin(2*phik)+Cuc*cos(2*phik);
    	deltark=Crs*sin(2*phik)+Crc*cos(2*phik);
    	deltaik=Cis*sin(2*phik)+Cic*cos(2*phik);
    
    %计算改正后的升交角距uk
    uk=phik+deltauk;
    
    %计算改正后的向径rk
    	rk=a*(1-e*cos(Ek))+deltark;
    
    %计算改正后的倾角ik
    ik=B5(1)+deltaik+B6(1)*tk;
    
    %计算观测瞬间升交点的经度Lk
    	Lk=B4(3)+(B5(4)-Ve)*tk-Ve*B4(1);
    
    %计算卫星在轨道平面内的坐标x,y,z
    	xyz=[rk*cos(uk);rk*sin(uk);0];
    
    %计算旋转矩阵的值
    	RzRx=[cos(Lk),-sin(Lk)*cos(ik),sin(Lk)*sin(ik);sin(Lk),cos(Lk)*cos(ik),-cos(Lk)*sin(ik);0,sin(ik),cos(ik)];
    
    %计算卫星在协议地球坐标系中的位置
    	XYZ=vpa([XYZ;1e-3*(RzRx*xyz)'],10);
    end
    
    XYZ