satellitepos.m

根据卫星星历及测者位置

function SatPos = SatellitePos(Satellite,GASTw,t)
%  计算t时刻卫星Satellite的空间位置,不考虑摄动项
%  Satellite = [ i0 Omiga0 a e w M0 ]; 　　　　其中6个分量为初始时刻t0时的初始值
%　　6个分量分别是轨道面倾角，升交点赤经，长半轴，轨道偏心率，近地点角距，平近地角
%  GASTw : 星期六午夜至星期天子夜的交换时刻格林尼治视恒星时
%  t = tk-t0;
%  i_gps = 55 * 2*pi/360; i_Galileo = 56 * 2*pi/360;
%  e = 0.02; % assumed
%  a_gps = 26562km; a_Galileo = 29584km; 
%  Ts_gps = 11h58m; Ts_Galileo =14h04m;
DtoR = 2*pi/360;
i0 = Satellite(1) * DtoR;
Omiga0 = Satellite(2) * DtoR;
a = Satellite(3);
e = Satellite(4);
w = Satellite(5) * DtoR;
M0 = Satellite(6) * DtoR;

% 计算平均角速度n
Miu = 3.986005e014;
n0 = sqrt( Miu / a^3 );
Delta_n = 0;
n = n0 + Delta_n;
% 计算归化时间tk
t0 = 0;    % 初始化时刻
tk = t - t0;
% 计算观测时刻卫星平近点角Mk
Mk = M0 + n * tk;
% 计算偏近点角Ek
% Ek=Mk+esin(Ek),故用迭代算法求出，e通常较小，所以2次迭代即可求出
Ek = Mk;
for j=1:5
    Ek = Mk + e * sin(Ek);
end
% 计算卫星矢径r
r = a * (1- e * cos(Ek));
% 计算真近点角fk
%%%%%% temp = ( sqrt(1-e^2) * sin(Ek) ) / ( cos(Ek) - e ) ;
%%%%%% fk = atan(temp);   % 该方法求得的fk只有-90度-90度
temp = sqrt((1+e)/1-e);
fk = 2 * atan( temp * tan(Ek/2) );
% 计算升交点角距Thita_k
Thita_k = fk + w;
% 计算摄动改正项Delta_u,Delta_r,Delta_i
Delta_u=0;
Delta_r=0;
Delta_i=0;
% 计算经过摄动改正的升交点角距uk，卫星矢径rk，轨道倾角ik
uk = Thita_k + Delta_u;
rk = r + Delta_r;
    Dot_i = 0;
ik = i0 + Delta_i + Dot_i*tk;
% 计算卫星在轨道平面直角坐标系中的位置
% 其中轨道平面直角坐标系X轴指向升交点
xk = rk * cos(uk);
yk = rk * sin(uk);
zk = 0;
SatPos = [xk yk zk]';  % 卫星在轨道平面直角坐标系中的位置
% 计算观测时刻升交点经度Omiga_k
Omiga = Omiga0;
we = 7.29211567e-005;        % 地球自转角速率
GAST = GASTw + we * tk;
Omiga_k = Omiga - GASTw;
% 计算卫星在地心固定坐标系统中的位置
R3 = [cos(Omiga) -sin(Omiga) 0;
      sin(Omiga) cos(Omiga) 0;
      0          0          1;];
R1 = [1          0          0;
      0          cos(ik)    -sin(ik);
      0          sin(ik)    cos(ik);];
R2 = [cos(w) -sin(w) 0;
      sin(w) cos(w)  0;
      0      0       1;];
SatPos = R3 * R1 * SatPos;  %  SatPos =[Xk Yk Zk]'