衛星星曆是描述衛星運動軌道的資訊。也可以說衛星星曆就是一組對應某一時刻的軌道引數及其變率。有了衛星星曆就可以計算出任意時刻的衛星位置及其速度。GPS衛星星曆分為預報星曆和後處理星曆。預報星曆又稱廣播星曆。

GPS廣播星曆引數共有16個，其中包括1個參考時刻，6個對應參考時刻的開普勒軌道引數和9個反映攝動力影響的引數。這些引數通過GPS衛星發射的含有軌道資訊的導航電文傳遞給使用者。

1.星曆參考時刻 ：

2.長半軸平方根：

3.偏心率：

4.參考曆元下平近點角：

5.近地點角距：

6.軌道傾角：

7.本週初始曆元的升交點赤經：

8.平運動差（由精密星曆計算得到的衛星平均角速度與按給定引數計算所得的平均角速度之差）：

9.軌道傾角變化率（弧度/秒）：

10.升交點赤經變化率（弧度/秒）：

11.緯度幅角的餘弦調和項改正的振幅（弧度）：

12.緯度幅角的正弦調和項改正的振幅（弧度）：

13.軌道半徑的餘弦調和項改正的振幅（m）：

14.軌道半徑的正弦調和項改正的振幅（m）：

15.軌道傾角的餘弦調和項改正的振幅（弧度）：

16.軌道傾角的正弦調和項改正的振幅（弧度）：

以下為matla程式：

%the homework 1 of chapater 1
%student ：Taelon
%time ：2018/9/25
%程式功能：根據所提供的星曆引數,計算此衛星在訊號發射時刻t( GPS時間)239050.7223s
%時的時空位置

%GPS廣播星曆引數共有16個，其中包括1個參考時刻，6個對應參考時刻的開普勒軌道引數
%和9個反映攝動力影響的引數。這些引數通過GPS衛星發射的含有軌道資訊的導航電文傳遞給使用者。
%時間引數
pra1 = 244800;           %te---星曆參考時刻，即星曆錶參考曆元(s)

%開普勒六引數
pra2 = 5153.65531;       %sqrt(a)長半軸平方根
pra3 = 0.005912038265;   %e---軌道偏心率
pra4 = -1.064739758;     %M0---按參考曆元te計算的平近點角（弧度）
pra5 = -1.717457876;     %w0---近地點角距
pra6 = 0.9848407943;     %i0---按參考曆元計算的軌道傾角(弧度)
pra7 = 1.038062244;      %Ω0---本週初始曆元的升交點赤經(弧度）

%軌道攝動九引數
pra8  = 4.249105564e-9;      % Δ?---平運動差（弧度）
pra9  = 7.422851197e-51;     % dot_i---軌道傾角變化率
pra10 = -8.151768125e-9;     % dot_Omega---升交點赤經變化率
pra11 = 3.054738045e-7;      % Cuc---升交點角距的調和改正項振幅
pra12 = 2.237036824e-6;      % Cus---升交點角距的調和改正項振幅
pra13 = 350.53125;           % Crc---衛星地心距的調和改正項振幅
pra14 = 2.53125;             % Crs---衛星地心距的調和改正項振幅
pra15 = -8.381903172e-8;     % Cic---軌道傾角的調和改正項振幅
pra16 = 8.940696716e-8;      % Cis---軌道傾角的調和改正項振幅

t=239050.7223;  %此衛星在訊號發射時刻t ( GPS時間)

miu = 3.986005e14;      % 地心引力常數
we = 7.2921151467e-5;   % 地球自轉角速度
F = -4.442807633e-10;   % Constant, [sec/(meter)^(1/2)]

% 計算歸化時間
te = pra1;
tk = t - te;
% 計算衛星的平均角速度
a = pra4.^2;  
n0 = sqrt(miu/a.^3);
n = n0 + pra8;
% 計算平近點角
Mk = pra2 + n * tk;
% 計算偏近點角Ek
%迭代計算：相鄰兩次計算差之絕對值值<1e-15時結束迭代計算,Ek的迭代初始值為0
Ek0=0;Ek=Mk;
while abs(Ek0-Ek)>1e-15
    Ek0 = Ek;
    Ek=Mk + pra3 * sin(Ek0);
end
% 計算衛星鐘差相對論校正值 Compute relativistic correction term

% 計算真近點角:取值在（-pi,pi]
v1=sqrt(1-pra3^2)*sin(Ek);
v2=cos(Ek)-pra3;
vk=atan(v1/v2);
if v1>0 & v2<0
Vk = pi + vk;%第二象限
elseif v1<0 & v2<0
Vk = vk - pi;%第三象限 
end
% 計算升交點角距
Faik = Vk + pra7;
% 計算升交點角距改正值
Sigmauk=pra11*cos(2*Faik)+pra12*sin(2*Faik);
% 計算衛星地心向徑改正值
Sigmark=pra13*cos(2*Faik)+pra14*sin(2*Faik);
% 計算衛星軌道傾角改正值
Sigmaik=pra15*cos(2*Faik)+pra16*sin(2*Faik);
% 修正升交點角距
uk=Faik+Sigmauk;
% 修正衛星地心相徑
rk=a*(1-pra3*cos(Ek))+Sigmark;
% 修正衛星軌道傾角
ik=pra6+Sigmaik+tk*pra9;
% 計算衛星軌道平面直角座標
xk=rk*cos(uk);
yk=rk*sin(uk);
% 計算升交點赤經
Omegak=pra5+(pra10-we)*tk-we*te;
% 將衛星座標由軌道平面轉換到ECEF座標系
Xk=xk*cos(Omegak)-yk*cos(ik)*sin(Omegak);
Yk=xk*sin(Omegak)+yk*cos(ik)*cos(Omegak);
Zk=yk*sin(ik);