function pos=SP(tk,brdc)
GM=3.986005e+14;
we=7.2921151467e-5;%地球自转平均角速度
a1=brdc(2,4);
a=a1*a1; % 长半轴
e=brdc(2,2);
w=brdc(4,3);
Crs=brdc(1,2);Cuc=brdc(2,1);Cus=brdc(2,3);
Cic=brdc(3,2);Cis=brdc(3,4);Crc=brdc(4,2);
% (1)已知长半轴a,计算平均角速度n
n0=sqrt(GM)/brdc(2,4)^3; % 平均角速度
n=n0+brdc(1,3); %改正平角速度
% (2)求偏近点角Ek
Mk=brdc(1,4)+n*tk; %平近点角
%迭代计算
Ek=Mk;Ek1=inf;
while abs(Ek1-Ek)>1e-13
Ek1=Ek;
Ek=Mk+e*sin(Ek1);
end
% (3)计算真近点角f
v1=sqrt(1-e^2)*sin(Ek);
v2=cos(Ek)-e;
fk=atan2(v1,v2); %真近点角
r0=brdc(2,4)^2*(1-e*cos(Ek));
% (4)根据已知的近地点角距w求升交角距theta
u0=fk+w; %升交距
fai=2*u0;
%(5)求改正数
duk=Cus*sin(fai)+Cuc*cos(fai); %轨道摄动项改正数
drk=Crs*sin(fai)+Crc*cos(fai);
dik=Cis*sin(fai)+Cic*cos(fai);
uk=u0+duk;
rk=r0+drk;
ik=brdc(4,1)+dik+brdc(5,1)*tk;
X=[rk*cos(uk);rk*sin(uk);0];
%计算升交点经度
OMk=brdc(3,3)+brdc(4,4)*tk-we*(brdc(3,1)+tk);
Rx=[1 0 0; 0 cos(ik) -sin(ik);0 sin(ik) cos(ik)];
Rz=[cos(OMk) -sin(OMk) 0;sin(OMk) cos(OMk) 0;0 0 1];
pos=Rz*Rx*X;
end
卫星 | 利用matlab根据星历读取卫星位置
最新推荐文章于 2024-07-15 11:45:50 发布