%
% Program to compute the three collinear Lagrangian points
%
% Note that the two equilateral L4 and L5 are given by 
%             X = rho - 0.5, Y = sqrt(3)/2
%

% input rho for earth-moon
%
%     rho = 0.01215;
rho =input(' input the mass ratio (M2/(M1+M2),  rho =  ');
%
%  common factor for all points 
double(rho);


c1 = 1.0;
c2 = 2 - 4*rho;
c3 = 1 -6*rho + 6*rho*rho;


% L3 points

c4 = -1 - 2*rho *(1-rho)^2 + 2*(1-rho)*rho^2;
c5 = rho^2*(1-rho)^2 +2*rho^2 -2*(1-rho)^2;
c6 = -rho^3 - (1-rho)^3;
coefL3=[c1 c2 c3 c4 c5 c6];
r3 = roots(coefL3);
%
nsize = max(size(r3));
for i=1:nsize,
    if (isreal(r3(i))==1)
        Lpoint = r3(i);
    end;
end;  
disp('  ');
disp(['Lagrangian L3 point  =   ', num2str(Lpoint)]);


% L2 points

c4 = 1 -2*rho -2*rho *(1-rho)^2 + 2*(1-rho)*rho^2;
c5 = rho^2*(1-rho)^2 +2*rho^2 +2*(1-rho)^2;
c6 = -rho^3 + (1-rho)^3;
coefL2=[c1 c2 c3 c4 c5 c6];
%
r2 = roots(coefL2);
nsize = max(size(r2));
for i=1:nsize,
    if (isreal(r2(i))==1)
        Lpoint = r2(i);
    end;
end;  
disp('  ');
disp(['Lagrangian L2 point  =   ', num2str(Lpoint)]);


% L1 points

c4 = 1  - 2*rho *(1-rho)^2 + 2*(1-rho)*rho^2;
c5 = rho^2*(1-rho)^2 -2*rho^2 +2*(1-rho)^2;
c6 = rho^3 + (1-rho)^3;
coefL1=[c1 c2 c3 c4 c5 c6];
%
r1 = roots(coefL1);
nsize = max(size(r1));
for i=1:nsize,
    if (isreal(r1(i))==1)
        Lpoint = r1(i);
    end;
end;  
disp('  ');
disp(['Lagrangian L1 point  =   ', num2str(Lpoint)]);

L4 =[-0.5+rho  sqrt(3)/2.];
L5 =[-0.5+rho  -sqrt(3)/2.];
disp('  ');
disp(['Lagrangian L4 point  =   (', num2str(L4), ')']);
disp('  ');
disp(['Lagrangian L5 point  =   (', num2str(L5), ')']);