%
%  w02planetX.m -- 15 jan 04 -- djm
%
%  - in search of lagrange points

%  coordinates

L = 15;  xr =  -L:0.05:L;  xi = xr;

[zr zi] = meshgrid(xr,xi);  zz = zr + i*zi;

%  first integral

GMs = 10.0;  Gme =  1.0;  R = 10;
As  = - R*Gme/(GMs+Gme);
ae  =   R*GMs/(GMs+Gme);
om2 = (GMs+Gme)/(R^3);

pot = om2*abs(zz).^2 + 2*GMs./abs(zz-As) + 2*Gme./abs(zz-ae);

force = om2*zz ... 
      - GMs*(zz-As)./(abs(zz-As).^3) - Gme*(zz-ae)./(abs(zz-ae).^3);

%  contourplot

figure(1);  clf;  hold on
contour(xr,xi,pot,[3.2:0.05:4]);
axis([-1 1 -1 1]*L);  
axis square;  colorbar
hold on
plot(As,0,'k*');  plot(ae,0,'k*')

ptitle = ['\bf potential energy & zero force contours'];
title(ptitle)
xlabel('\bf Re(x)')
ylabel('\bf Im(x)')

contour(xr,xi,real(force),[0 0],'k');
contour(xr,xi,imag(force),[0 0],'k--');