%% runKepler.m  script for 2-body problem
% 
%%
%
% $$x''(t)=-x(t)/(x(t)^2+y(t)^2)^{3/2}$$
% $$y''(t)=-y(t)/(x(t)^2+y(t)^2)^{3/2}$$
%
close all
tinit=0;
tFin=2*pi;
tSpan=[tinit tFin];
% u=[x, x', y, y']'
% Initial values:
uN=[1 0 0 .4]';  %  x(0)=1, x'(0)=0, y(0)=0, y'(0)=.4  % Init velo North
uNE=[1 .4 0 .4]';  %  x(0)=1, x'(0)=.4, y(0)=0, y'(0)=.4  % Init velo NE
%uW=[1 -.4 0 0]';  %  x(0)=1, x'(0)=-.4, y(0)=0, y'(0)=0  % Init velo W
%uE=[1 .4 0 0]';  %  x(0)=1, x'(0)=.4, y(0)=0, y'(0)=0  % Init velo E
% Odefun: Kepler.m
%% Solve using ode23
u0=uN;
u0=uNE;
%u0=uW;   %Failure, stop early
%u0=uE;  %Failure
[T,U]=ode23(@Kepler,tSpan,u0);
% 
% Orbit (x(t),y(t)) is (u1(t),u3(t)):
plot(U(:,1),U(:,3),'--.')
axis equal
grid on
hold on
plot(0,0,'*g','MarkerSize',15)  % Earth
plot(0,0,'og','MarkerSize',15)  % Earth
plot(u0(1),u0(3),'xr','MarkerSize',10)  % Initial point.
%%  Animations
% comet(U(:,1),U(:,3))  % Easy, no options, here too fast
%% More flexible, a bit fast
h=animatedline('Marker','.','Color','b','Markersize',10);
for k=1:length(U(:,1))
   addpoints(h,U(k,1),U(k,3))
   drawnow
end

%% Draw line segments, most flexible, needs no new function, just plot.
%clf
% [T,U]=ode23(@Kepler,[0 1],uE);
% axis([-1.5 1.5 -1.5 1.5])

for j=1:length(U(:,1))
    plot([0 U(j,1)],[0 U(j,3)])
    pause(0.05)   % adjust
end  
%
%%
length(T)
N=floor(length(T)/2)
for j=1:N
    plot([U(j,1)],[U(j,3)],'*r')
    pause(0.2)   % adjust
end 

%% Step lengths, adaptivity ode23
steps=diff(T)
format long
[min(steps) max(steps)]
% ans = 0.000012800000000   0.403309852408788
figure
plot(T(2:end), diff(T))

%% Same with ode45
clear; close all
%%
tinit=0;
tFin=2*pi;
tSpan=[tinit tFin];
% u=[x, x', y, y']'
% Initial values:
u0=[.4 0 0 2]';  %  x(0)=.4, x'(0)=0, y(0)=0, y'(0)=2
% Odefun: Kepler.m
%% Solve using ode45
figure
[T,U]=ode45(@Kepler,tSpan,u0);
% 
% Orbit (x(t),y(t)) is (u1(t),u3(t)):
plot(U(:,1),U(:,3),'--.')
axis equal
grid on
hold on
plot(0,0,'*','MarkerSize',10)  % Earth
plot(0,0,'o','MarkerSize',10)  % Earth
%% More flexible, a bit fast
figure
axis([-.2 1.2 -.4 .5])
h=animatedline('Marker','.','Color','b','Markersize',10);
for k=1:length(U(:,1))
    addpoints(h,U(k,1),U(k,3))
    drawnow
end
%% Stepsizes
figure
steps=diff(T);
format long
[min(steps) max(steps)]
% ans = 0.000012800000000   0.403309852408788
figure
plot(T(2:end), diff(T))
%%
sol=ode45(@Kepler,tSpan,u0)
%% Evaluating the sol at arbitrary t-points
%
% 1) Tspan=linspace(0,2*pi);
% 2) sol structure
% 3) spline interpolation
Tspan=[0 2*pi]
u0=[1 .2 0 .4]
sol=ode23(@Kepler,Tspan,u0)
%{
sol = 

  struct with fields:

     solver: 'ode23'
    extdata: [1×1 struct]
          x: [1×54 double]
          y: [4×54 double]
      stats: [1×1 struct]
      idata: [1×1 struct]
%}
x0=u0(1);y0=u0(3);
T=linspace(Tspan(1),Tspan(2),500);
% u=solfun(sol,T)    % WHY?
u=deval(sol,T)
x=u(1,:);
y=u(3,:);
d=(x-x0).^2+(y-y0).^2
dmax=max(d)
T(d==dmax)
max(T)
ind=1:length(T);
maxind=ind(d==dmax)
dmin=min(d(maxind:end))
minind=ind(d==dmin)
%%
Tmin=T(d==dmin)
max(T)
%%
plot(x,y,'*')
%%
figure(2)
plot(x(1:minind),y(1:minind),'r.')
hold on
p%% Eccentricity:
max(U(:,[1 3]))-min(U(:,[1 3]))
e=ans(2)/ans(1)
%{
ans =
    2.0010    1.6002
%}lot(x(minind),y(minind),'*k')

%%
grid on
figure(1)
plot(x(end),y(end),'k*')
%%
plot(0,0,'o')
axis([-.2 1.2 -.2 .4])
for j=1:round(minind/2)
    plot([x(j)],[y(j)],'*r')
    pause(0.02)   % adjust
end

minind/2


%%
Y=deval(sol,0:2);
%{
Y =
   0.400000000000000  -0.628533592651708  -1.339519375606119
                   0  -0.982363224767453  -0.465860957390635
                   0   0.800083909830288   0.538714881331783
   2.000000000000000  -0.022591401529652  -0.410042830267516
%}
T=linspace(0,2*pi);Y=deval(sol,T);
plot(Y(1,:),Y(3,:))
shg
% Solution is now rowwise.
%% Example of time dpendent params: (doc ode45)
%
% $$ y'' =A/B ty$$
% 
% Write as 1st order syst:
% y1=y, y2=y'=y1'
%  Then:
% y1' = y2
% y2' = A/B*ty1
%
% odefcn.m represents this system of equations as a function that accepts 
% four input arguments: t, y, A, and B.

% function dydt = odefcn(t,y,A,B)
% dydt = zeros(2,1);
% dydt(1) = y(2);
% dydt(2) = (A/B)*t.*y(1);

% Solve the ODE using ode45. Specify the function handle such that it passes
% in the predefined values for A and B to odefcn.

% A = 1;
% B = 2;
% tspan = [0 5];
% y0 = [0 0.01];
% [t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);
%%
close all
options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)
[T,U]=ode23(@Kepler,tSpan,u0,options);

%{
218 successful steps
0 failed attempts
655 function evaluations
%}

%%
options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)
[T,U]=ode45(@Kepler,tSpan,u0,options);
%{
33 successful steps
4 failed attempts
223 function evaluations
%}
%% P. 349: Same Apollo-problem with many more questions.

