% Simulation der Lagrangeschen Punkte

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% Version 3.3  2013

% Simulation der Lagrangeschen Punkte
r0 = 1;
m2 = 1;
m1 = 2; %330000; Erde - Sonne
G = 1;
M = m1 + m2;
%
x = linspace(-1.5*r0,1.25*r0,150);
y = linspace(-1.25*r0,1.25*r0,150);
[X,Y] = meshgrid(x,y);
%
r13 = sqrt(Y.^2 + (m2/M*r0 - X).^2);
r23 = sqrt(Y.^2 + (m1/M*r0 + X).^2);
rq = X.^2 + Y.^2;
omega2 = G*M/r0^3;   % omega^2
U = -G*m1./r13-G*m2./r23-omega2*rq/2;
%
Umax = max(U(:));
Umin = min(U(:));
Ucut = (49*Umax + Umin)/50;
U(U<Ucut) = NaN;
figure, surf(X,Y,U),shading interp,xlabel('x'),ylabel('y'),colormap lines, view(-9.5,84),axis tight,shg
figure, contour(X,Y,U,100),shg, axis equal
hold on
%
x1 = m2/M * r0;
x2 = -m1/M * r0;
xe = (m2-m1)/(2*M) * r0;
ye = sqrt(3)/2 * r0;
plot([x2,x1],[0,0],'k')
plot([x2,xe],[0,ye],'k')
plot([x2,xe],[0,-ye],'k')
plot([x1,xe],[0,ye],'k')
plot([x1,xe],[0,-ye],'k')
xlabel('x'),ylabel('y')
hold off
%%
%xl=(m2-m1)/(2*M)*r0;
%yl=linspace(-1,2*r0);
%r13l = sqrt(yl.^2 + (m2/M*r0 - xl).^2);
%r23l = sqrt(yl.^2 + (m1/M*r0 + xl).^2);
%rql = xl.^2 + yl.^2;
%Ul = -G*m1./r13l-G*m2./r23l-omega2*rql/2;
%figure,plot(yl,Ul),shg
