Ionospheric error delay.rar_gatherq2l_suchmnm_对流层仿真_延时测量_电离层

function [pos,dtr,Orbit_parameter] = coe2pos(t, coe) % pos 卫星位置ECEK坐标和速度 [Xk; Yk; Zk;Vx;Vy;Vz] % dtr 卫星钟差相对论校正值 % Orbit_parameter 轨道参数输出，可以用于画卫星空间坐标图 % t=239050.7223; %此卫星在信号发射时刻t ( GPS时间) % %---转赋参数值------ % te = 244800; % 星历参考时刻 % a = 5153.65531.^2; % 长半轴 % e = 0.005912038265; % 偏心率 % M0 = -1.064739758; % 参考历元下平近点角 % w0 = -1.717457876; % 近地点角距 % i0 = 0.9848407943; % 轨道倾角 % Omega0 = 1.038062244; % 本周初始历元的升交点赤经 % deta_n = 4.249105564e-9; % 平运动差 % dot_i = 7.422851197e-51; % 轨道倾角变化率 % dot_Omega = -8.151768125e-9; % 升交点赤经变化率 % Cuc = 3.054738045e-7; % 升交点角距的调和改正项振幅 % Cus = 2.237036824e-6; % 升交点角距的调和改正项振幅 % Crc = 350.53125; % 卫星地心距的调和改正项振幅 % Crs = 2.53125; % 卫星地心距的调和改正项振幅 % Cic = -8.381903172e-8; % 轨道倾角的调和改正项振幅 % Cis = 8.940696716e-8; % 轨道倾角的调和改正项振幅 te = coe(1); % 星历参考时刻 a = coe(2).^2; % 长半轴 e = coe(3); % 扁率 M0 = coe(4); % 参考历元下平近点角 w0 = coe(5); % 近地点角距 i0 = coe(6); % 轨道倾角 Omega0 = coe(7); % 本周初始历元的升交点赤经 deta_n = coe(8); % 平运动差 dot_i = coe(9); % 轨道倾角变化率 dot_Omega = coe(10); % 升交点赤经变化率 Cuc = coe(11); % 升交点角距的调和改正项振幅 Cus = coe(12); % 升交点角距的调和改正项振幅 Crc = coe(13); % 卫星地心距的调和改正项振幅 Crs = coe(14); % 卫星地心距的调和改正项振幅 Cic = coe(15); % 轨道倾角的调和改正项振幅 Cis = coe(16); % 轨道倾角的调和改正项振幅 miu = 3.986005e14; % 地心引力常数 we = 7.2921151467e-5; % 地球自转角速度 F = -4.442807633e-10; % Constant, [sec/(meter)^(1/2)] %---1）计算卫星位置 % 计算归化时间 tk = t - te; if tk > 302400 tk = tk - 604800; elseif tk < -302400 tk = tk + 604800; end % 计算卫星的平均角速度 n = sqrt( miu/a.^3 ) + deta_n; % 计算平近点角 Mk = M0 + n*tk; % 计算偏近点角 Ek = 0; while 1 % 迭代 Ekold=Ek; Ek = Mk + e*sin(Ek); if abs(Ek - Ekold) < 1e-15 break; end end %Compute relativistic correction term dtr = F * e * sqrt(a) * sin(Ek); % 计算真近点角 vk = 2 * atan( sqrt( (1+e)/(1-e) ) * tan(Ek/2) ); % 计算升交点角距 fai = w0 + vk; % 计算升交点角距改正值 deta_uk = Cuc*cos(2*fai) + Cus*sin(2*fai); % 计算卫星地心向径改正值 deta_rk = Crc*cos(2*fai) + Crs*sin(2*fai); % 计算卫星轨道倾角改正值 deta_ik = Cic*cos(2*fai) + Cis*sin(2*fai); % 修正升交点角距 uk = fai + deta_uk; % 修正卫星地心相径 rk = a*(1-e*cos(Ek)) + deta_rk; % 修正卫星轨道倾角 ik = i0 + dot_i*tk + deta_ik; % 计算卫星坐标 xk = rk*cos(uk); yk = rk*sin(uk); % 计算升交点赤经 Omegak = Omega0 + (dot_Omega-we)*tk - we*te; % 将卫星坐标由轨道平面转换到ECEF坐标系 Xk = xk*cos(Omegak) - yk*sin(Omegak)*cos(ik); Yk = xk*sin(Omegak) + yk*cos(Omegak)*cos(ik); Zk = yk*sin(ik); %-----2)卫星速度----- %计算信号发射时刻的偏近点角导数 Mkd = n; Ekd = Mkd/(1- e*cos(Ek)); % 计算升交点角距导数 vkd=sqrt(1-e^2)*Ekd/(1-e*cos(Ek)) faid = vkd; % 计算各摄动校正项导数 % 计算升交点角距改正项导数 deta_ukd = 2*faid*(Cus*cos(2*fai) - Cuc*sin(2*fai)); % 计算向径改正项导数 deta_rkd =2*faid*(Crs*cos(2*fai) - Crc*sin(2*fai)); % 计算轨道倾角改正项导数 deta_ikd =2*faid*(Cis*cos(2*fai) -Cic*sin(2*fai)); % 计算升交点角距项导数 ukd = faid + deta_ukd; % 计算向径项导数 rkd = a*e*Ekd*sin(Ek) + deta_rkd; % 计算轨道倾角项导数 ikd = dot_i + deta_ikd; % 计算升交点赤经项导数 Omegakd = dot_Omega - we; % 计算卫星坐标项导数 xkd = rkd*cos(uk)-rk*ukd*sin(uk); ykd = rkd*sin(uk)+rk*ukd*cos(uk); %计算卫星在WGS-84地心地固直角坐标系中的速度 Vx = -Yk*Omegakd - (ykd*cos(ik)-Zk*ikd)*sin(Omegak)+xkd*cos(Omegak); Vy = Xk*Omegakd + (ykd*cos(ik)-Zk*ikd)*cos(Omegak)+xkd*sin(Omegak); Vz = ykd*sin(ik)+yk*ikd*cos(ik); pos = [Xk; Yk; Zk;Vx;Vy;Vz]; % Mk(Rad): Mean Anomaly % Ek(Rad): Eccentric Anomaly % vk(Rad): True Anomaly % fai(Rad): Argument of lattitude % uk(Rad):Corrcted argument of lattitude % rkd(meter): corrected radius % ikd(Rad):corrected inclination % Omegak(Rad):Corrected longitude of ascending node. % a(meter):semi-major orbit axis % e: orbit eccentricity Orbit_parameter=[Mk Ek vk fai uk rk ik Omegak a e]; % 轨道参数输出，可以用于画卫星空间坐标图 end