function lagrange
% Plot the effective stationary potential field of a restricted 3-body system
% in the rotating frame
G = 6.67259e-11;
M1 = 10; % mass of body 1
M2 = 1; % mass of body 2
mu1 = G*M1;
mu2 = G*M2;
alpha = M2/(M1+M2);
beta = M1/(M1+M2);
r1 = [-alpha 0 0]; % position of body 1
r2 = [beta 0 0]; % position of body 2
R = alpha+beta;
omega = sqrt((mu1+mu2)/R^3); % angular velocity of rotating frame
Omega = [0 0 omega];
grid=[200 200];
range=[[-1.5*R -1.5*R];[1.5*R 1.5*R]];
X = ([1:grid(1)]-1) * (range(2,1)-range(1,1))/(grid(1)-1)+range(1,1);
Y = ([1:grid(2)]-1) * (range(2,2)-range(1,2))/(grid(2)-1)+range(1,2);
U = zeros(grid);
for i=1:grid(1)
for j=1:grid(2)
p = [X(i) Y(j) 0];
% Gravitational potentials
R1 = norm(p-r1);
R2 = norm(p-r2);
U1 = -mu1/R1;
U2 = -mu2/R2;
% Centrifugal potential
oa = cross(Omega,p);
Uc = 1/2 * dot(oa,oa);
U(i,j) = U1+U2-Uc;
end
end
contourf(X,Y,rot90(max(U,-1.5e-9)),20);
axis equal tight
hold on
% The following plot of the Lagrange points 1-3 requires the Optimisation
% toolbox. Otherwise you need to replace the fzero function with your own
% root finder.
for L=1:5
[x,y] = CalcLagrange(L);
plot(x,y,'*')
end
function [x,y] = CalcLagrange(which)
if which <= 3
switch which
case 1
s0=-1;
s1=1;
case 2
s0=1;
s1=1;
case 3
s0=-1;
s1=-1;
end
u=fzero(@lagpoly,0);
x=R*(u+beta);
y=0;
else
x=R/2*(M1-M2)/(M1+M2);
y=sqrt(3)/2*R;
if which==5
y=-y;
end
end
function of=lagpoly(u)
of = alpha*(s0+2*s0*u+(1+s0-s1)*u^2+2*u^3+u^4) - u^2*((1-s1)+3*u+3*u^2+u^3);
end
end
end