clear all;
clc;
global u
u=3.9860044e5;
%定义轨道根数OE
i=0.3673;
OMEGA=5.2888;
omega=0.1198;
p=1.1071e4;
e=0.6270;
tao=-951.9128;
%已知轨道根数反解r和v
a=p/(1-e^2); %反解a
M=(u/a^3)^0.5*(-tao); %反解M
syms E
E=solve(M-E+e*sin(E)); %反解E
cosf=(cos(E)-e)/(1-e*cos(E));
r=a*(1-e^2)/(1+e*cosf);
T=2*pi*(a^3/u)^0.5; %计算周期用来判断sinf符号
if M>T/2 %计算sinf
sinf=-(1-cosf^2)^0.5
else
sinf=(1-cosf^2)^0.5
end
syms vr vs xp yp vxp vyp
vr=(u/(a*(1-e^2)))^0.5*e*sinf;
vs=(u/(a*(1-e^2)))^0.5*(1+e*cosf);
xp=r*cosf;
yp=r*sinf;
VJ=[vr;vs;0];
VP=[vxp;vyp;0];
XP=[xp;yp;0];
A1=[cosf sinf 0;-sinf cosf 0;0 0 1];
VP=A1\VJ;
A21=[cos(omega) sin(omega) 0;-sin(omega) cos(omega) 0;0 0 1]; %求转换矩阵
A22=[1 0 0;0 cos(i) sin(i);0 -sin(i) cos(i)];
A23=[cos(OMEGA) sin(OMEGA) 0;-sin(OMEGA) cos(OMEGA) 0;0 0 1];
A2=A21*A22*A23;
V=A2\VP; %反解速度
vx=double(V(1,1));
vy=double(V(2,1));
vz=double(V(3,1));
X=A2\XP; %反解位置
x=double(X(1,1));
y=double(X(2,1));
z=double(X(3,1));
R0=[x y z vx vy vz];
%用代数方程画出Orbit1
syms MM
for ii=1:101
t(ii)=((ii-1)*T)/100;
MM(ii)=(u/(a^3))^0.5*(t(ii)-tao);
syms EE
fun=MM(ii)-EE+e*sin(EE);
EE=solve(fun); %直接调用solve求E
cosff(ii)=(cos(EE)-e)/(1-e*cos(EE)); %求cosf
%求直角坐标中卫星坐标
rr(ii)=a*(1-e^2)/(1+e*cosff(ii));
xx(ii)=rr(ii)*cosff(ii);
yy(ii)=a*(1-e^2)^0.5*sin(EE); %rsinf=bsinE
RR=[xx(ii);yy(ii);0]; %求Orbit3坐标
RRR=A2\RR;
XXX(ii)=RRR(1,1);
YYY(ii)=RRR(2,1);
ZZZ(ii)=RRR(3,1);
end
%画Orbit1
figure(1)
hold on;
grid on %画网格
title('Orbit1') %主题
xlabel('x/km') %x轴含义和单位
ylabel('y/km') %y轴含义和单位
plot(xx(:),yy(:),'r-',xx(1),yy(1),'or');hold on;
h1=line('Color','r','LineStyle','-','Marker','o','Markersize',5,'EraseMode','xor');
for kk=1:101
set(h1,'XData',xx(kk),'YData',yy(kk)); %轨迹
drawnow; %画出当前一帧
pause(0.05) %暂停
end
%用微分方程画出Orbit2
hh=T/100; %积分步长(秒)
options=odeset('RelTol',1e-12,'AbsTol',1e-15); %积分的误差选项
[t,iR]=ode45('Differential',[0:hh:T],R0,options); %积分结果放在iR中
%画Orbit2
figure(2)
hold on;
plot3(iR(:,1),iR(:,2),iR(:,3),'b',iR(1,1),iR(1,2),iR(1,3),'ob');
title('Orbit2');
xlabel('X/km') %x轴含义和单位
ylabel('Y/km') %y轴含义和单位
zlabel('Z/km') %z轴含义和单位
grid on;
row=length(iR); %求数组的长度
for nn=1:row %求Orbit4坐标
RRRR=[iR(nn,1);iR(nn,2);iR(nn,3)];
RRRRR=A2*RRRR;
xxxx(nn)=RRRRR(1,1);
yyyy(nn)=RRRRR(2,1);
end
h2=line('Color','b','LineStyle','-','Marker','o','Markersize',5,'EraseMode','xor');
for jj=1:row
set(h2,'XData',iR(jj,1),'YData',iR(jj,2),'ZData',iR(jj,3)); %卫星
drawnow; %画出当前一帧
pause(0.05) %暂停
end
%画Orbit3
figure(3)
hold on;
grid on %画网格
title('Orbit3') %主题
xlabel('X/km') %X轴含义和单位
ylabel('Y/km') %Y轴含义和单位
zlabel('Z/km') %Z轴含义和单位
plot3(XXX(:),YYY(:),ZZZ(:),'r-',XXX(1),YYY(1),ZZZ(1),'or');hold on;
h3=line('Color','r','LineStyle','-','Marker','o','Markersize',5,'EraseMode','xor');
for pp=1:101
set(h3,'XData',XXX(pp),'YData',YYY(pp),'ZData',ZZZ(pp)); %轨迹
drawnow; %画出当前一帧
pause(0.05) %暂停
end
%画Orbit4
figure(4)
hold on;
plot(xxxx(:),yyyy(:),'b',xxxx(1),yyyy(1),'ob');
title('Orbit4');
xlabel('x/km') %x轴含义和单位
ylabel('y/km') %y轴含义和单位
grid on;
h4=line('Color','b','LineStyle','-','Marker','o','Markersize',5,'EraseMode','xor');
for qq=1:row
set(h4,'XData',xxxx(qq),'YData',yyyy(qq)); %卫星
drawnow; %画出当前一帧
pause(0.05) %暂停
end
%Orbit2&3比较
figure(5)
hold on;
grid on; %画网格
title('Orbit2&3') %主题
xlabel('X/km') %x轴含义和单位
ylabel('Y/km') %y轴含义和单位
zlabel('Z/km') %y轴含义和单位
plot3(iR(:,1),iR(:,2),iR(:,3),'b',iR(1,1),iR(1,2),iR(1,3),'ob');hold on;
plot3(XXX(:),YYY(:),ZZZ(:),'r-',XXX(1),YYY(1),ZZZ(1),'or');hold on;
h5=line('Color','b','LineStyle','-','Marker','o','Markersize',5,'EraseMode','xor');
h6=line('Color','r','LineStyle','-','Marker','o','Markersize',10,'EraseMode','xor');
for ss=1:row
set(h5,'XData',iR(ss,1),'YData',iR(ss,2),'ZData',iR(ss,3)); %轨迹
set(h6,'XData',XXX(ss),'YData',YYY(ss),'ZData',ZZZ(ss));
drawnow; %画出当前一帧
pause(0.05) %暂停
end
%Orbit1&4比较
figure(6)
hold on;
grid on; %画网格
title('Orbit1&4') %主题
xlabel('x/km') %x轴含义和单位
ylabel('y/km') %y轴含义和单位
plot(xx(:),yy(:),'r-',xx(1),yy(1),'or');hold on;
plot(xxxx(:),yyyy(:),'b',xxxx