%lagrangian_points

function lagrangian_points

%% ASTRONOMICAL DATA %%

%Astronomical masses /kg
m_earth = 5.972e24;
m_moon = 7.34767309e22;
m_sun=1.99e30;
m_jupiter = 1.898e27;
m_pluto = 1.30900e22;
m_charon = 1.55e21;

%Gravitational force constant /m^3 kg^-1 s^-2
G = 6.67e-11;

%Earth-Sun distance /m
AU = 149597870700;

%Earth to moon distance /m
LD = 384400e3;

%Pluto to Charon distance /m
PC = 19640e3;

%Sun to Jupiter distance /m
SJ = 5.202*AU;

%

%% INPUTS %%

%Masses /kg;
M1 = m_pluto;
M2 = m_charon;

%Mass separation /m
a = PC;

%Range of plot
xmin = -2*a;
xmax = 2.2*a;
ymin = -2*a;
ymax = 2*a;

%Gravitational potential scale
gscale = G*M1/a;

%Minimum phi/gscale to plot
phi_scaled_min = -2.3;

%Number of points for surface
N = 1000;

%Fontsize for surface
fsize = 10;

%

%Compute matrix of x,y values
x = linspace(xmin,xmax,N);
y = linspace(ymin,ymax,N);
[x,y] = meshgrid(x,y);

%Determine angular speed of M2 rotationa about M1 /rads^-1
omega = sqrt( (M1 + M2)*G )/(a^1.5);

%Compute gravitational potential + term to take into account rotating
%reference frame
phi = -G*M1./sqrt( x.^2 + y.^2 ) - G*M2./sqrt( (x-a).^2 + y.^2 ) -...
    0.5*(omega^2)*( ( x - a/(1+M1/M2) ).^2 + y.^2 );
phi( phi/gscale < phi_scaled_min ) = NaN;

%Determine staionary points
[xs,ys,phis,maxormin] = stationary_points(x,y,phi);

%

%Plot surface of phi(x,y)
figure('name','Two body gravitational potential','color',[1 1 1],'renderer','opengl')
s = surf(x/a,y/a,phi/gscale);
hold on;
if isempty(xs)==0
    for n=1:numel(xs)
        if maxormin(n)==1
            %Local maxima
            plot3(xs(n)/a,ys(n)/a,phis(n)/gscale,'w+');
        elseif maxormin(n)==-1
            %Local minima
            plot3(xs(n)/a,ys(n)/a,phis(n)/gscale,'m+');
        else
            %Saddle point
            plot3(xs(n)/a,ys(n)/a,phis(n)/gscale,'k+');
        end          
    end
end
caxis([phi_scaled_min,0])
shading interp
colormap prism
axis equal
axis tight
view([-43,54]);
grid on
box on
set(gca,'fontsize',fsize);
xlabel('x','fontsize',fsize)
ylabel('y','fontsize',fsize)
title(['x scale is ',num2str(a,4),' m, M1=',num2str(M1,4),' kg, M2=',num2str(M2,4),' kg'],'fontsize',fsize)
colorbar
interp_colormap(N);
print(gcf,'lagrangian 3D.png','-dpng','-r300')

%

%Plot 2D image of phi(x,y)
figure('name','Two body gravitational potential','color',[1 1 1],'renderer','opengl')
pcolor(x/a,y/a,phi/gscale);
hold on;
if isempty(xs)==0
    for n=1:numel(xs)
        if maxormin(n)==1
            %Local maxima
            plot(xs(n)/a,ys(n)/a,'w+');
        elseif maxormin(n)==-1
            %Local minima
            plot(xs(n)/a,ys(n)/a,'m+');
        else
            %Saddle point
            plot(xs(n)/a,ys(n)/a,'k+');
        end          
    end
end
caxis([phi_scaled_min,0])
shading interp
colormap prism
axis equal
axis tight
grid on
box on
set(gca,'fontsize',fsize);
xlabel('x','fontsize',fsize)
ylabel('y','fontsize',fsize)
title(['x scale is ',num2str(a,4),' m, M1=',num2str(M1,4),' kg, M2=',num2str(M2,4),' kg'],'fontsize',fsize)
colorbar
interp_colormap(N);
print(gcf,'lagrangian 2D.png','-dpng','-r300')

%Export .stl file of surface
stlwrite( 'lagrange.stl',surf2patch(s,'triangles'));

%%

%stationary_points
% Function which determines stationary points on a 2D surface defined by
% matrices x,y,z, where z = z(x,y). Outputs are xs,ys,zs coordinates, with
% maxormin yield: -1 minima, 1 maxima, 0 saddle
function [xs,ys,zs,maxormin] = stationary_points(x,y,z)

%Step through matrix z and determine whether a point is a local maxima,
%local minima or saddle
xs = [];
ys = [];
zs = [];
maxormin = [];
dim = size(z);
for n=1:numel(z)
    %Determine (row,column) coordinates
    [row,col] = ind2sub( dim,n );
    
    %Determine whether current location is a stationary point
    if  (row==1) || (row==dim(1)) || (col==1) || (col==dim(2))
        %Edge of z matrix, ignore
    else
        if (z(row,col) >= z(row-1,col)) && (z(row,col) >= z(row+1,col)) && (z(row,col) >= z(row,col-1)) && (z(row,col) >= z(row,col+1))
            %Local maxima
            maxormin = [maxormin,1];
            xs = [xs,x(row,col)];
            ys = [ys,y(row,col)];
            zs = [zs,z(row,col)];
        elseif (z(row,col) <= z(row-1,col)) && (z(row,col) <= z(row+1,col)) && (z(row,col) <= z(row,col-1)) && (z(row,col) <= z(row,col+1))
            %Local minima
            maxormin = [maxormin,-1];
            xs = [xs,x(row,col)];
            ys = [ys,y(row,col)];
            zs = [zs,z(row,col)];
        elseif (z(row,col) <= z(row-1,col)) && (z(row,col) <= z(row+1,col)) && (z(row,col) >= z(row,col-1)) && (z(row,col) >= z(row,col+1))
            %Saddle
            maxormin = [maxormin,0];
            xs = [xs,x(row,col)];
            ys = [ys,y(row,col)];
            zs = [zs,z(row,col)];
        elseif (z(row,col) >= z(row-1,col)) && (z(row,col) >= z(row+1,col)) && (z(row,col) <= z(row,col-1)) && (z(row,col) <= z(row,col+1))
            %Saddle
            maxormin = [maxormin,0];
            xs = [xs,x(row,col)];
            ys = [ys,y(row,col)];
            zs = [zs,z(row,col)];
        end
    end   
end

%%

%interp_colormap
% Function which interpolates current colourmap to yield better graduated
% shading. N is number of possible colours.
function new_map = interp_colormap(N)

%Get current colourmap
map = colormap;

%Initialise new colormap
new_map = ones(N,3);

%Get size of current colormap and initalise red,green,blue vectors
dim = size(map);
R = ones(1,dim(1));
G = ones(1,dim(1));
B = ones(1,dim(1));
RR = ones(1,N);
GG = ones(1,N);
BB = ones(1,N);

%Populate these with current colormap
R(:) = map(:,1);
G(:) = map(:,2);
B(:) = map(:,3);

%Interpolate to yield new colour map
x = linspace( 1, dim(1), N );
RR = interp1( 1:dim(1), R, x );
GG = interp1( 1:dim(1), G, x );
BB = interp1( 1:dim(1), B, x );
new_map(:,1) = RR(:);
new_map(:,2) = GG(:);
new_map(:,3) = BB(:);

%Set colormap to be new map
colormap( new_map );

%End of code


