%% Apollosolutioninteractive.m  
%  14.3.2018 
%%
close all
%%
mu=1/82.45;  % Note: These values don't pass to the function Apollosys
lam=1-mu;
w0=[1.2;0;0;-1.04935751];
% w0=[1.2;0;0;-0.5];  % Another set of initial values, try your own! 
Tperiod=6.19217;
%tmax=4*pi;
Tspan=[0 Tperiod];
[T,U]=ode45(@Apollosys,Tspan,w0);
%%  Plots of orbit
hold on
x=U(:,1); y=U(:,3);     % Orbit
plot(x,y);shg  
plot(w0(1),w0(3),'*b')    % Initial point. (b=blue)
% plot(x(1),y(1),'*b')    % Same as the previous line.
plot(x(end),y(end),'ko')  % Final point (k=black)
plot(-mu,0,'*k',lam,0,'or') % plot earth and moon  (r=red)
%legend('Orbit','Initial pt.','Final pt.','Earth','Moon')
grid on
%axis equal
%% Text to plots (start with  gtext, then change to text)
%gtext('Earth')        % Get text coordinates with mouse.
 text(-.2,-.1,'Earth')
% gtext('Moon')
 text(.77,-.1,'Moon')

%% Try tolerances:
options1=odeset('AbsTol',1e-6,'RelTol',1e-8);
[T1,U1]=ode45(@Apollosys,Tspan,w0,options1);
%% New figure
%
figure
clf
hold on
x1=U1(:,1); y1=U1(:,3); % Orbit: [x,y]; 
plot(x1,y1);shg  
plot(w0(1),w0(3),'*b')    % Initial point.
% plot(x1(1),y1(1),'*b')  % Same as the previous line.
plot(x1(end),y1(end),'ko') % Final point 
% Now is periodic at the accuracy of pic.
plot(-mu,0,'*k',lam,0,'or') % plot earth and moon
%legend('Orbit','Initial pt.','Final pt.','Earth','Moon')
grid on
% Comparison:
% plot(x,y,'--m');% First orbit, default tolerances
%% "Direction field"
% Direction field doesn't quite make sense here, as we have 4 initial
% values. But interactive choise of initial points does, just fix the
% initial velocities first (and vary them).
%  Initial velocity:
dx=0;
dy=-1.04935751;
xy0 = [];  % Gather chosen initial points
n = 0;
%axis([-2 2 -2 2])
grid on
hold on
% Loop, picking up the points.
disp('Move cursor into graphics window')
disp('******************************')
disp(' ')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
tmax=10;
tmax=Tperiod;
tmax=4;
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
[T,U]=ode45(@Apollosys,[0 tmax],[xi;dx;yi;dy],options1); 
    plot(U(:,1),U(:,3))
    n = n+1;
    text(xi,yi,num2str(n))  % Place numbers at chosen initial points 
    xy0(:,n) = [xi;yi];
end
%% Initial points saved 
% Now that they are numbered, you can pick those that give interesting
% trajectories.
% plot(xy0(1,:),xy0(2,:),'r.')% Not needed because of numbers at init.pts.
%% Plot chosen trajectories:
for k=[1 3 6]
    xk=xy0(1,k);
    yk=xy0(2,k);
    [T,U]=ode45(@Apollosys,[0 tmax],[xk;dx;yk;dy],options1);
    plot(U(:,1),U(:,3),'k')
end




