Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Single shooter Matlab code for finding Halo orbits (*not* validated!)
function state = halofunc(X)
% Parse state vector into variables
x = X(1);
y = X(2);
z = X(3);
vx = X(4);
vy = X(5);
vz = X(6);
Phi = reshape(X(7:end), 6,6);
% Simplication
mu = Constants.mu;
r1 = sqrt((x+mu)^2 + y^2 + z^2);
r2 = sqrt((x+mu-1)^2 + y^2 + z^2);
ax = x + 2*vy - (1-mu)*(x+mu)/r1^3 - mu*(x+mu-1)/r2^3;
ay = y - 2*vx - (1-mu)*y/r1^3 - mu*y/r2^3;
az = -(1-mu)*z/r1^3 - mu*z/r2^3;
% Compute STM
A = A_matrix(X(1:6)', mu);
Phid = A*Phi;
state = [vx; vy; vz; ax; ay; az; reshape(Phid, 6^2, 1)];
clc
clear
close all
format long
X0 = [1.14 0 -0.16 0 -0.223 0;];
% Found?
X0 = [1.140000000000000 0 2.791036381471859 0 -1.282994906583491 0]
it = 0;
tol = 1e-8;
while 1
fprintf('iteration %d\n', it)
tspan = [0 50];
Phi0 = eye(6);
ic = [X0, reshape(Phi0, 1, 6^2)];
options = odeset('Events', @periodevent, 'RelTol', 1e-13, 'AbsTol', 1e-13);
[T,X] = ode45(@(t,X)halofunc(X), tspan, ic, options);
atXing = X(end, 1:6)
if atXing(4) < tol && atXing(6) < tol
disp(atXing)
disp(X0)
disp(T(end))
plotfunc(X,'Rotating', 'HW8')
break
else
% Let's fix this orbit.
Phi = reshape(X(end, 7:end), 6, 6);
dXdt = halofunc(X(end, :)');
deltaXZDot = [-atXing(4); -atXing(6)];
el1 = [Phi(4, 3) Phi(4, 5); Phi(6, 3) Phi(6, 5)];
el2 = (1/dXdt(2)) * [dXdt(4);dXdt(6)]*[Phi(2, 3) Phi(2, 5)];
rhsMat = el1 - el2;
deltaState = deltaXZDot' * inv(rhsMat);
%deltaState = inv(rhsMat) * deltaXZDot;
X0(3) = X0(3) + deltaState(1);
X0(5) = X0(5) + deltaState(2);
end
it = it + 1;
if it > 10
disp('INCOMPLETE')
plotfunc(X,'Rotating', 'HW8')
break
end
end
function [value, isterminal, direction] = periodevent(t, state)
%value = state(2);
value = t - 1.51074761956563;
isterminal = 1;
direction = -1;
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.