如何击落美国间谍卫星？

% 2d kepler problem

clear all; close all; clc;

T = 20;

xstart = [2;0];

xend = [0 ;2];

iter = 10;

v00 = [1;0.7];

vincrement = 0.001;

relax = 1;

% euler-lagrange equation

N2 = 1e5;

dt = T/N2;

for ss = 1 : iter

x0 = xstart;             % initial positioN2

v0 = v00;

y = [x0; v0];

ylist = zeros(4, N2+1);

Elist = zeros(1, N2+1);

ylist(:,1) = y;

Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

for s = 2: (N2+1)

ytemp = y;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K1;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K2;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + K3;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

ylist(:,s) = y;

Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

end

diff = ylist(1:2, end) - xend;

%         h3 = figure;

plot(ylist(1,:), ylist(2,:), xstart(1), xstart(2),'*', xend(1), xend(2),'o')

axis equal

pause(0.1)

Mov(ss) = getframe;

x0 = xstart;             % initial position

v0 = v00 + [vincrement;0];

y = [x0; v0];

ylist = zeros(4, N2+1);

Elist = zeros(1, N2+1);

ylist(:,1) = y;

Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

for s = 2: (N2+1)

ytemp = y;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K1;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K2;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + K3;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

ylist(:,s) = y;

Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

end

diff1 = ylist(1:2, end) - xend;

x0 = xstart;             % initial position

v0 =v00+ [0; vincrement];

y = [x0; v0];

ylist = zeros(4, N2+1);

Elist = zeros(1, N2+1);

ylist(:,1) = y;

Elist(1,1) = 0.5*norm(v0)^2 - 1/norm(x0);

for s = 2: (N2+1)

ytemp = y;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K1 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K1;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K2 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + 0.5*K2;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K3 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

ytemp = y + K3;

r = sqrt(ytemp(1)^2 + ytemp(2)^2);

K4 = dt*[ytemp(3); ytemp(4); -ytemp(1)/r^3; -ytemp(2)/r^3];

y = y + (K1+ 2*K2 + 2*K3 + K4)/6;

ylist(:,s) = y;

Elist(1,s) = 0.5*norm(y(3:4))^2 - 1/norm(y(1:2));

end

diff2 = ylist(1:2, end) - xend;

%  v01 =  v00 -[vincrement , 0 ; 0, vincrement] * ([diff1 - diff, diff2- diff]\diff)

v00 = v00 - relax * ([diff1 - diff, diff2 - diff]\([vincrement , 0 ; 0, vincrement]* diff))

end

% movie(Mov)

movie2gif(Mov, 'lambert.gif','DelayTime',0.3)

http://wap.sciencenet.cn/blog-100379-1181810.html

全部精选博文导读

GMT+8, 2021-7-25 02:02