if ~exist('OCTAVE_VERSION')
	error('Function definitions in scripts available only in Octave. Edit the source code, extract functions to separate files and rerun');
end

clear;
disp('Problem 1 from midterm 2014');

global r1;
global r2;
global h;

r1 = @(x,y,mu) sqrt((x+mu).^2 + y.^2); % these are squares of dist
r2 = @(x,y,mu) sqrt((x-(1-mu)).^2 +y.^2);

function w = omega(x,y,mu)
% compute omega at rectangular gridpoints (x,y) (e.g. returned by meshgrid)
global r1;
global r2;
R1 = r1(x,y,mu);
R2 = r2(x,y,mu);
w = -0.5*((1-mu)*R1.^2 + mu*R2.^2) -(1-mu)*(1./R1) - (mu)*(1./R2); 
end

function dwx = dxomega(x,y,mu)
global r1;
global r2;
global h; % for FD approx only
R1 = r1(x,y,mu);
R2 = r2(x,y,mu);
dxR1 = (x+mu)./R1;
dxR2 = (x-(1-mu))./R2;
if h == 0
	% exact formula: takes effort to derive
	dwx = -(1-mu)*R1.*dxR1-mu*R2.*dxR2 + (1-mu)*(1./R1.^2).*dxR1 + (mu)*(1./R2.^2).*dxR2;
else
	% quick and dirty approximation
	dwx = (omega(x+h,y,mu) - omega(x,y,mu))/h;
end
end

function dwy = dyomega(x,y,mu)
global r1;
global r2;
global h; % for FD approx only
R1 = r1(x,y,mu);
R2 = r2(x,y,mu);
dyR1 = y./R1;
dyR2 = y./R2;
if h == 0
	dwy = -(1-mu)*R1.*dyR1-mu*R2.*dyR2 + (1-mu)*(1./R1.^2).*dyR1 + (mu)*(1./R2.^2).*dyR2;
else 
	dwy = (omega(x,y+h,mu) - omega(x,y,mu))/h;
end
end

function dX = gravity(t,U,mu)
x = U(1); y = U(2); vx = U(3); vy = U(4);
dx = vx;
dvx = 2*vy - dxomega(x,y,mu);
dy = vy;
dvy = -2*vx - dyomega(x,y,mu);
dX = [dx;dy;dvx;dvy];
end

h = 1.0*sqrt(eps);

mu = 0.012277471;
% plot 
N = 500;
A = 1.5; B = 1.5; 
x = linspace(-A,A,N);
y = linspace(-B,B,N);
[X,Y] = meshgrid(x,y);
Z = omega(X,Y,mu);
contour(X,Y,Z);
title('\Omega(x,y)');
ZX = dxomega(X,Y,mu);
ZY = dyomega(X,Y,mu);
contour(X,Y,ZX.^2+ZY.^2);

mu = 0.012277471;
x0 = 0.994; vx0 = 0;
y0 = 0; 

vy0 = -2.0015;
T =  17.065;

if h ~= 0 
	warning('Using FD approx for omega derivatives');
end
t = linspace(0,T,1000);
U = lsode(@(U,t)gravity(t,U,mu), [x0,y0,vx0,vy0]', t);
plot(U(:,1),U(:,2)); hold on;
plot([mu,1-mu],[0,0], 'r*'); hold off;

% restiricting to 1-D problem:
f1 = @(x) norm(gravity(0,[x;0;0;0],mu),Inf);
M = 1000; x = linspace(-1.5,1.5, M); f1v = [];
for i=1:M
	f1v(i) = f1(x(i));
end
semilogy(x,f1v);

function x = findlagrange(X,mu)
opts.TolFun = 1e-12;
[x, fx, info] = fsolve(@(U)gravity(0,U,mu), X(:), opts);
norm(fx)
if info ~= 1 
	warning(sprintf('x(1) = %g can be inaccurate: fsolve exit code: %d\n', x(1), info));
end
end

u = []; 
u = [u,findlagrange([0.836, 0, 0, 0]',mu)];
u = [u,findlagrange([1.156, 0, 0, 0]',mu)];
u = [u,findlagrange([-1, 0, 0, 0]',mu)];
u = [u,findlagrange([0, -1, 0, 0]',mu)];
u = [u,findlagrange([0, +1, 0, 0]',mu)];

u
plot(U(:,1),U(:,2)); hold on;
plot([mu,1-mu],[0,0], 'r*'); 
plot(u(1,:),u(2,:), 'g+'); 
hold off;
legend('satellite','planets','lagrange points');

warning('Accuracy falsification and aliasing check - omitted here');
