clc

clear

A=[ 244800            5153.65531        0.005912038265     0.984848407943

1.038062244       -1.717457876      -1.064739758       4.249105564e-9

7.422851197e-51   -8.151768125e-9   3.054738045e-7     2.237036824e-6

350.53125         2.53125            -8.381903172e-8    8.940696716e-8];

t=239050.7223;

GM=3.986005e+14;

tk=t-A(1,1);                               %规划时间

while(tk>302400 || tk

if tk>302100

tk=tk-604800;

else

tk=tk+604800;

end

end

n0=sqrt(GM/(A(1,2)^6));

n=n0+A(2,4);                               %卫星平均角速度

Mk=A(2,3)+n*tk;                            %卫星平近点角

while(Mk<0 || Mk>2*pi)

if Mk<0

Mk=Mk+2*pi;

else

Mk=Mk-2*pi;

end

end

Ek=Mk;

for i=1:3                                    %卫星偏近点角

Ek=Mk+A(1,3)*sin(Ek);

end

% cosv=(cos(Ek)-A(1,3))/(1-A(1,3)*cos(Ek));

% sinv=(sqrt(1-A(1,3)^2)*sin(Ek))/(1-A(1,3)*cos(Ek));

vk=atan(sqrt(1-A(1,3)^2)*sin(Ek)/(cos(Ek)-A(1,3)));   %卫星真近点角

if sin(Ek)>0 && (cos(Ek)-A(1,3))<0

vk=vk+pi;

elseif sin(Ek)<0 && (cos(Ek)-A(1,3))<0

vk=vk-pi;

end

w=vk+A(2,2);                                          %卫星升交点角距

uuk=A(3,4)*sin(2*w)+A(3,3)*cos(2*w);                  %摄动校正项

rrk=A(4,2)*sin(2*w)+A(4,1)*cos(2*w);

iik=A(4,4)*sin(2*w)+A(4,3)*cos(2*w);

uk=w+uuk;                       %校正后的升交点角距，卫星矢径长度，轨道倾角

rk=(A(1,2)^2)*(1-A(1,3)*cos(Ek))+rrk;

ik=A(1,4)+A(3,1)*tk+iik;

xxk=rk*cos(uk);                    %轨道平面的位置

yyk=rk*sin(uk);

wk=A(2,1)+(A(3,2)-7.2921151467e-5)*tk-A(1,1)*7.2921151467e-5;   %升交点赤经

while (wk<0 || wk>2*pi)

if wk<0

wk=wk+2*pi;

else

wk=wk-2*pi;

end

end

disp('卫星瞬时位置：')

xk=xxk*cos(wk)-yyk*cos(ik)*sin(wk)          %WGS-84

yk=xxk*sin(wk)+yyk*cos(ik)*cos(wk)

zk=yyk*sin(ik)

%计算瞬时速度

mk=n;                   %平近点角

ek=mk/(1-A(1,3)*cos(Ek));             %偏近点角

Wk=sqrt(1-A(1,3)^2)*ek/(1-A(1,3)*cos(Ek));           %升交点角距

Uuk=2*Wk*(A(3,4)*cos(2*w)-A(3,3)*sin(2*w));       %摄动校正项

Rrk=2*Wk*(A(4,2)*cos(2*w)-A(4,1)*sin(2*w));

Iik=2*Wk*(A(4,4)*cos(2*w)-A(4,3)*sin(2*w));

Uk=Wk+Uuk;                                      %校正后

Rk=(A(1,2)^2)*A(1,3)*ek*sin(Ek)+Rrk;

Ik=A(3,1)+Iik;

WK=A(3,2)-7.2921151467e-5;           %升交点赤经

% while (WK<0 || WK>2*pi)

%

% if WK<0

%     WK=WK+2*pi;

% else

%     WK=WK-2*pi;

% end

%

% end

Xxk=Rk*cos(uk)-rk*Uk*sin(uk);

Yyk=Rk*sin(uk)+rk*Uk*cos(uk);

disp('卫星瞬时速度：')

Xk=-yk*WK-(Yyk*cos(ik)-zk*Ik)*sin(wk)+Xxk*cos(wk)

Yk=xk*WK+(Yyk*cos(ik)-zk*Ik)*cos(wk)+Xxk*sin(wk)

Zk=Yyk*sin(ik)+yyk*Ik*cos(ik)

点赞

1

收藏

分享

文章举报

王特起

发布了3 篇原创文章 · 获赞 3 · 访问量 1544

私信

关注