GPS_gps星历_计算卫星位置_星历计算_卫星星历_星历

%程序功能：根据所提供的星历参数,计算此卫星在信号发射时刻t( GPS时间)239050.7223s %时的时空位置 %GPS广播星历参数共有16个，其中包括1个参考时刻，6个对应参考时刻的开普勒轨道参数 %和9个反映摄动力影响的参数。这些参数通过GPS卫星发射的含有轨道信息的导航电文传递给用户。 %时间参数 pra1 = 244800; %te---星历参考时刻，即星历表参考历元(s) %开普勒六参数 pra2 = -1.064739758; %M0---按参考历元te计算的平近点角（弧度） pra3 = 0.005912038265; %e---轨道偏心率 pra4 = 5153.65531; %sqrt(a)长半轴平方根 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);