close all
clc
K = 80;
step = 250;
t_max = 8;
nt = 5000;
v_0x = 3;
v_0y = 0;
x0 = -10;
y0 = 3;
t = linspace(0,t_max,nt);
dt = t(2) - t(1);
vP = 0; xP0 = 2;
vP = -2; xP0 = 10;
xP = xP0 + vP * t;
yP = zeros(1,nt);
x(1) = x0; y(1) = y0;
xd(1) = x(1) - xP(1);
Rd(1) = sqrt(xd(1)^2 + y(1)^2);
vx(1) = v_0x; vy(1) = v_0y;
ax(1) = -(K / abs(Rd(1))^3) * xd(1);
ay(1) = -(K / abs(Rd(1))^3) * y(1);
x(2) = x(1) + vx(1) * dt + 0.5 * ax(1) * dt^2;
y(2) = y(1) + vy(1) * dt + 0.5 * ay(1) * dt^2;
xd(2) = x(2) - xP(2);
Rd(2) = sqrt(xd(2)^2 + y(2)^2);
vx(2) = vx(1) + ax(1) * dt;
vy(2) = vy(1) + ay(1) * dt;
for c = 3 : nt
x(c) = -K*dt^2*xd(c-1)/Rd(c-1)^3 + 2*x(c-1) - x(c-2);
y(c) = -K*dt^2*y(c-1)/Rd(c-1)^3 + 2*y(c-1) - y(c-2);
xd(c) = x(c) - xP(c);
Rd(c) = sqrt(xd(c)^2 + y(c)^2);
end
R = x + j .* y;
temp_p = zeros(1,nt+2); temp_m = zeros(1,nt+2); temp = zeros(1,nt+2);
temp_p(1:end-2) = R; temp_m(3:end)= R; temp = temp_p - temp_m;
v(1) = abs((R(2)-R(1))/dt); v(nt) = abs((R(end-1)-R(end-2))/dt);
v(2:nt-1)= abs(temp(3:end-2)./(2*dt));
ET = -100*ones(1,nt);
EK = v.^2;
EP = ET- EK;
figure(1)
x_p = real(R);
y_p = imag(R);
ns = nt /50;
title_x = 'x position (a.u.)';
title_y = 'y position (a.u.)';
lineWidth_p = 1.5;
plot(x_p,y_p,'lineWidth',lineWidth_p);
hold on
hp = plot(x_p(1:step:end),y_p(1:step:end),'o','MarkerSize',4);
set(hp,'MarkerEdgeColor',[0 0 0],'MarkerFaceColor', [0 0 0]);
xlabel(title_x);
ylabel(title_y);
axis([ XXXXXXXXXX])
axis equal
box on
hold on
hp = plot(xP(1:step:end),yP(1:step:end),'o','MarkerSize',4);
set(hp,'MarkerEdgeColor',[1 0 0],'MarkerFaceColor', [1 0 0]);
set(gca,'Xtick',[ XXXXXXXXXX]);
set(gca,'Ytick',[ XXXXXXXXXX]);
set(gca,'XLim',[-10 10]);
set(gca,'YLim',[-10 10]);
figure(2);
subplot(2,1,1);
plot(t,v,'linewidth',3);
grid on
xlabel('time t (a.u.)');
ylabel('speed of spacecraft (a.u.)');
subplot(2,1,2);
plot(t,abs(R),'linewidth',3);
grid on
xlabel('time t (a.u.)');
ylabel('separation distance (a.u.)');
figure(4);
plot(t,EK,'r','linewidth',3);
hold on
plot(t,EP,'b','linewidth',3);