gps02_1.rar_alm2skyplot_antswtu_gps02_1_广播星历_星历计算

%-------------------------------------------------------------------------- % 定义常量 %-------------------------------------------------------------------------- GM_WGS84 = 3986005e+8; %WGS84引力常数 omega_e_WGS84 = 7.2921151467e-5; %WGS84地球自转角速度 %-------------------------------------------------------------------------- % 星历电文中给出的变量 %-------------------------------------------------------------------------- iode = 0.0; %电文中给出的当前参考历元的有效期 sqrt_a = 0.515365263176e+4; %电文中给出的卫星轨道椭圆长半轴的平方根 delta_n = 0.451411660250e-8; %电文中给出的平地点角改正值 toe = 0.72e+4; %电文中给出的参考时刻（周积秒） M_zero = -0.290282040486; %电文中给出的参考时刻平近点角 e1 = 0.678421219345e-2; %电文中给出的轨道椭圆偏心率 omega = -0.258419417299e+1; %电文中给出的轨道近地点角距 Cus = 0.912137329578e-5; %电文中给出的升交点赤经的改正项—正弦振幅 Cuc = 0.189989805222e-6; %电文中给出的升交点赤经的改正项—余弦振幅 Cis = 0.949949026108e-7; %电文中给出的倾角角距的改正项—正弦振幅 Cic = 0.130385160446e-7; %电文中给出的倾角角距的改正项—余弦振幅 Crs = 0.40625e+1; %电文中给出的轨道半径角距的改正项—正弦振幅 Crc = 0.201875e+3; %电文中给出的轨道半径角距的改正项—余弦振幅 i_zero = 0.958512160302; %电文中给出的参考时刻轨道倾角 i_dot = -0.253939149013e-9; %电文中给出的轨道倾角变化率 OMEGA_zero = -0.137835982556e+1; %电文中给出的参考时刻升交点赤经 OMEGA_dot = -0.819426989566e-8; %电文中给出的升交点赤经变化率 %-------------------------------------------------------------------------- % 需要计算的变量 %-------------------------------------------------------------------------- t = 0.84e+4; %观测时刻（1997年11月9号2时20分0秒的周积秒） t_k = 0.0; %观测时刻到参考时刻的归化时间（周积秒） M_k = 0.0; %观测时刻的平近点角 n = 0.0; %平均角速度 l_k = 0.0; %真近点角 phi_k = 0.0; %升角距角 C_u = 0.0; %升交点赤经摄动改正 C_i = 0.0; %升交点赤经摄动改正 C_r = 0.0; %升交点赤经摄动改正 u_k = 0.0; %经摄动改正后的升交距角 r_k = 0.0; %经摄动改正后的卫星矢径 i_k = 0.0; %经摄动改正后的轨道倾角 OMEGA_k = 0.0; %观测时刻升交点经度 x_k = 0.0; %轨道平面x坐标 y_k = 0.0; %轨道平面y坐标 X_k = 0.0; %WGS84坐标系X坐标 Y_k = 0.0; %WGS84坐标系Y坐标 Z_k = 0.0; %WGS84坐标系Z坐标 %-------------------------------------------------------------------------- % 计算观测时刻GPS卫星在WGS84坐标系的坐标 %-------------------------------------------------------------------------- n_zero = sqrt(GM_WGS84) / (sqrt_a * sqrt_a * sqrt_a); n = n_zero + delta_n; % 计算平均角速度 t_k = t - toe; while (t_k > 302400) t_k = t_k - 604800; end while (t_k < -302400) t_k = t_k + 604800; end % 计算观测时刻到参考时刻的归化时间 M_k = M_zero + n * t_k; % 观测时刻的卫星平近点角 E_k = M_k; E_k1 = M_k - 1e-8; while (abs(E_k - E_k1) > 1e-9) E_k1 = E_k; E_k = M_k + e1 * sin(E_k1); end; % 计算观测时刻的偏近点角 l_k = 2*atan(sqrt((1+e1)/(1-e1)) * tan(E_k/2)); %计算真近点角 phi_k = l_k + omega; % 计算升交距角 C_u = Cus*sin(2*phi_k) + Cuc*cos(2*phi_k); C_r = Crs*sin(2*phi_k) + Crc*cos(2*phi_k); C_i = Cis*sin(2*phi_k) + Cic*cos(2*phi_k); %计算摄动改正项 u_k = phi_k + C_u; r_k = sqrt_a*sqrt_a*(1-e1*cos(E_k)) + C_r; i_k = i_zero + i_dot*t_k + C_i; %计算摄动改正后的升交距角、卫星矢径、轨道倾角 x_k = r_k * cos(u_k); y_k = r_k * sin(u_k); %计算卫星在轨道坐标系中的位置 OMEGA_k = OMEGA_zero + (OMEGA_dot-omega_e_WGS84)*t_k - omega_e_WGS84*toe; %计算观测时刻卫星的升交点经度 X_k = x_k*cos(OMEGA_k) - y_k*cos(i_k)*sin(OMEGA_k); Y_k = x_k*sin(OMEGA_k) + y_k*cos(i_k)*cos(OMEGA_k); Z_k = y_k*sin(i_k); %计算卫星在WGS84坐标系中的位置 report = sprintf('1997年11月9号2时20分0秒卫星的WGS84坐标：\n X = %7.3f米\n Y = %7.3f米\n Z = %7.3f米',X_k,Y_k,Z_k); disp(report);