GPS广播星历计算地固坐标系坐标

 参考连接：GPS卫星计算小程序  

prn = 1;
sqrt_a = 5153.7127704599998; %轨道半长轴平方根
delta_n = 4.1115998360700002e-009;
a0 = 0.00028600683435800001 ;
a1 = 1.7053025658199999e-012 ;
a2 = 0 ;
tgc = 4800.0 ; %观测时间 
toc = 0.00000000000000000 ; %卫星钟基准时间
M0 = 1.2263973009600000; %参考时刻的平近点角
e = 0.0053100715158500003;% 轨道离心率
omega = -1.6819446292500000; %近地点角距   rad
Cuc = -5.5264681577700003e-006;%轨道延迹方向上周期改正余弦振幅
Cus = 1.1192634701700000e-005; %轨道延迹方向上周期改正正弦振幅
Crc = 175.34375000000000; %在轨道径向方向上周期改正余弦的振幅
Crs = -105.43750000000000; %对轨道半径正弦的修正值
Cic = -9.6857547759999998e-008; %轨道倾角周期改正余弦项振幅
Cis = -7.8231096267699997e-008; %轨道倾角周期改正正弦项振幅
i0 = 0.97432927738800001; %参考时刻TOE的轨道倾角,
I = 1.8643633724999999e-010; %轨道倾角变化率
Omega_dot = -7.7124641122299999e-009; %升交点赤径在赤道平面中的长期变化
Omega0 = -2.9080127721900002; %参考时刻升交点赤径
toe=0.00000000000000000; % 星历基准时间
[X_ECEF,Y_ECEF,Z_ECEF]=SatellitePosition(prn,sqrt_a,delta_n,a0,a1,a2,tgc,toc,M0,e,omega,Cuc,Cus,Crc,Crs,Cic,Cis,i0,I,Omega_dot,Omega0,toe)

 参考连接：公式及代码，代码1，星历信息解读1、星历信息解读2

%% 广播星历解读
% GPS以1980年1月6日0时0分0秒开始计时，用周数和周内秒数来表示。
% 1-8行头文件，后面是卫星相关参数
% prn %卫星号
% sqrt_a        轨道长半轴的平方根
% delta_n       平均角速度改正项
% a0             卫星钟差
% a1             卫星钟数
% a2             卫星钟数变化率
% tgc            观测时刻对应的GPS周内秒
% toc            参考时刻对应的GPS周内秒
% M0             参考时刻TOE时的平近点角  rad
% e              轨道偏心率
% omega          近点角距   rad
% Cuc            升交点角距改正项  rad
% Cus            升交点角距改正项  rad
% Crc            轨道半径余弦改正项 m
% Crs            轨道半径正弦改正项 rad
% Cic            轨道倾斜角预先的改正项 rad
% Cis            轨道倾斜角正弦改正项 rad
% i0             参考时刻TOE的轨道倾角
% I              轨道倾角i的变化率
% Omega_dot      升交点赤经变化
% Omega0         升交点赤经
% toe             星历参考时刻
function [X_ECEF,Y_ECEF,Z_ECEF]=SatellitePosition(prn,sqrt_a,delta_n,a0,a1,a2,tgc,toc,M0,e,omega,Cuc,Cus,Crc,Crs,Cic,Cis,i0,I,Omega_dot,Omega0,toe)
%% 1.计算卫星运动的平均角速度
% 参考时刻TOE的平均角速度n0
GM = 3.986005e14 ;% m^3/s^2
n0 = sqrt(GM)/sqrt_a^3;

% 观测时刻的平均角速度
n = n0 +delta_n ; 
 
%% 2.计算信号发射时卫星的平近点角 ???什么是平近点角
delta_t = a0 + a1*(tgc - toc) + a2*(tgc - toc)^2;
t = tgc - delta_t ;%经过卫星时钟改正的时刻
tk = t - toc; %归化时间
M = M0 + n*tk; % M即为信号发射时的平近点角  rad

%% 3.计算偏近点角 和 4.计算真近点角
% E     偏近点 E = M+ e*sinE   用迭代的方法求E 
% 初始值M 迭代3 次即可
E = zeros(1,5);
for i = 1:5
    E(i) = M+ e*sin(M) ;
    M=E(i);
end
delta_temp3 = E(3)-E(2)
delta_temp2 = E(2)-E(1)

% 真近点角
num_temp = sqrt(1-e^2)*sin(E(1));
den_temp = cos(E(1)) - e;
Vk = atan2(num_temp,den_temp); %rad

%% 5.计算升交距角 和 6.计算摄动改正项
 u = omega + Vk ;
 % 计算摄动改正项
 delta_u = Cuc*cos(2*u) + Cus*sin(2*u); %升交距角的摄动改正项
 delta_r = Crc*cos(2*u) + Crs*sin(2*u); %卫星失径的摄动改正项
 delta_i = Cic*cos(2*u) + Cis*sin(2*u); %卫星轨道倾角的摄动改正项
 
 %% 7.计算摄动改正后的升交距角、卫星失径和轨道倾角
 a = sqrt_a^2 ;
 uk = u + delta_u ; %摄动改正后的升交距角
 rk = a*(1-e*cos(E(3)))+delta_r;
 ik = i0 + delta_i + I*tk;
 
 %% 8.计算卫星在轨道面坐标系中的坐标
 x = rk * cos(uk) ;
 y = rk * sin(uk) ;
 
 %% 9.计算发射时刻升交点的经度
 omega_e = 7.29211567e-5 ; % omega_e 地球自转速度 7.29211567e-5 rad/s
 L = Omega0 + Omega_dot*tk - omega_e*(tk + toe);
 
 %% 10.计算卫星在地固坐标系下坐标
 % 地心地固坐标系（Earth-Centered, Earth-Fixed，简称ECEF）
 X_ECEF = x*cos(L)-y*cos(ik)*sin(L);
 Y_ECEF = x*sin(L)+y*cos(ik)*cos(L);
 Z_ECEF = y*sin(ik);
 
end