Skip to content

Instantly share code, notes, and snippets.

@Hugo-Diaz-N
Last active Oct 19, 2018
Embed
What would you like to do?
Code to model a spacecraft is given by two 2nd-order equations for the position of the body (x(t), y(t)) --- Numerical Method: Störmer-Verlet (Dual) -- Matlab
%---------------- Setting Problem --------------------%
x0 = 0.994; % x(0) Initial position
y0 = 0.0; % y(0) Initial position
dx0 = 0.0; % x'(0) Initial velocity
dy0 = -2.031732629557; % y'(0) Initial velocity
%-----------------------------------------------------%
%---------------- Setting Model ----------------------%
mu1 = 0.012277471; % Moon mass, M_m/(M_e+M_m)
mu2 = 1-mu1; % Earth mass, M_e/(M_e+M_m)
%-----------------------------------------------------%
D1 = @(x,y) ( (x+mu1).^2+y.^2 ).^(1.5);
D2 = @(x,y) ( (x-mu2).^2+y.^2 ).^(1.5);
r = @(x,y) 1-( mu1./D2(x,y)+mu2./D1(x,y) );
g = @(x,y) mu1*mu2*(1./D2(x,y)-1./D1(x,y));
f = @(q,v) [ r(q(1),q(2))*q(1)+g(q(1),q(2))+2*v(2);...
r(q(1),q(2))*q(2)-2*v(1)];
l = @(q) [ r(q(1),q(2))*q(1)+g(q(1),q(2));...
r(q(1),q(2))*q(2)];
%-----------------------------------------------------%
dt = 0.0001;
t = 0;
T = 11.12455; % period
N = ceil( T/dt)+2;
A = (1/(dt^2+1))*[1 dt;-dt 1];
q = zeros(2,N); v = zeros(2,N);
q(:,1)=[x0;y0]; v(:,1)=[dx0;dy0];
for n=2:N
vn12 = v(:,n-1);
qn12 = q(:,n-1);
qn = qn12+(dt/2)*vn12;
v(:,n) = A*(vn12+(dt/2)*(f(qn,vn12)+l(qn)));
q(:,n) = qn+(dt/2)*v(:,n);
end
plot(q(2,:),q(1,:),'-.b','LineWidth',2) % plot(y,x)
hold on
plot(0,-mu1,'ro',0,mu2,'bo','MarkerSize', 12)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment