%
% 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
%
%  common factor for all points 
double(rho);
format long, coefL3;
format long, coefL2;
format long, coefL1;


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)

% 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)

% 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)

L4 =[-0.5+rho  sqrt(3)/2.]

L5 =[-0.5+rho  -sqrt(3)/2.]