% REFERENCES:
% [1] "Symplectic Integrators" from the MATLAB File Exchange
%      http://www.mathworks.com/matlabcentral/fileexchange/7686

clear;
clf;
cla;
hold off;


n=1;
t_max = 35;
time_steps=350;


% GENERAL SETUP
q0=1.6;
p0=0.5;

% Time interval
tspan = linspace(0,t_max,time_steps);

% SYMPLECTIC SOLVER

% Setup

dpdt=@(q)(-sin(q));

dqdt=@(p)(p);

tic,
[se.q,se.p]=seep(dqdt,dpdt,tspan,q0,p0); % seep is in the package [1]
toc
fprintf(' (seep)\n');

% EXPLICIT (NON-SYMPLECTIC) Euler
[ee.q,ee.p]=ee(dqdt,dpdt,tspan,q0,p0);

% INBUILT SOLVER

% Setup
odefun=@(t,y)([y(2);-sin(y(1))]);
y0=[q0;p0];

% Solver options
opt=odeset('Stats','on');

% Solve
tic,
[~,y]=ode45(odefun,tspan,y0,opt); 
toc 
fprintf(' (ode45)\n');

% Convert to p,q format
% evens=2*(1:n);
% odds=2*(1:n)-1;
% ns.q=y(:,odds);
% ns.p=y(:,evens);
ode45ns.q=y(:,1);
ode45ns.p=y(:,2);

%plot(nonsymp.q,nonsymp.p,'-r',symp.q,symp.p,'-b');


vo = VideoWriter('short.avi','Uncompressed AVI'); % video object
%Settings
vo.FrameRate=30;
%vo.Quality=100;
open(vo); % open up the video file
% 
% 
% Axis settings
set(gca,'Color',[1 1 1]);
set(gca,'nextplot','replacechildren');
axis square;
axis([-2 2 -2 2]);
title('Phase space evolution; y''''+sin(y) = 0')
xlabel(gca,'q (= y)');
ylabel('p (= y'')');

text('Position',[0.5,0], 'String','Symplectic Euler','Color',[0 0 1]);
text('Position',[0,0], 'String','ode45','Color',[1 0 0]);
text('Position',[-1,0], 'String','Explicit Euler','Color',[0 0 0]);

% 
%lsh=plot(ns.q,ns.p,'.r',s.q,s.p,'.b');
%set(lsh(1),'MarkerSize',1.0);
%set(lsh(2),'MarkerSize',1.0);

for i=1:(time_steps-1)
    line(ode45ns.q(i:(i+1)),ode45ns.p(i:(i+1)),'Color',[1 0 0]);
    line(se.q(i:(i+1)),se.p(i:(i+1)),'Color',[0 0 1]);
    line(ee.q(i:(i+1)),ee.p(i:(i+1)),'Color',[0 0 0]);
%     if i == 1
%         % this is a hack. 
%         % to plot separate line series, we need matrices with >1 row
%         % below, we delete the extra row
%         h=plot(ns.q(1:2,:),ns.p(1:2,:),'-r',s.q(1:2,:),s.p(1:2,:),'-b');
%         ns.h=h(1:n);
%         s.h=h((n+1):2*n);
%         axis square;
%         axis([-2 2 -2 2]);
%         
%         % delete the values for the second time step
%         for j=1:length(h) 
%             hand=h(j);
%             xd=get(hand,'XData');
%             yd=get(hand,'YData');
%             xd(end)=[];
%             yd(end)=[];
%             set(hand,'XData',xd);
%             set(hand,'YData',yd);            
%         end
%     else
%         for j=1:length(ns.h) 
%             hand=ns.h(j);
%             xd=get(hand,'XData');
%             yd=get(hand,'YData');
%             xd(i)=ns.q(i,j);
%             yd(i)=ns.p(i,j);
%             set(hand,'XData',xd);
%             set(hand,'YData',yd);            
%         end
%         for j=1:length(s.h) 
%             hand=s.h(j);
%             xd=get(hand,'XData');
%             yd=get(hand,'YData');
%             xd(i)=s.q(i,j);
%             yd(i)=s.p(i,j);
%             set(hand,'XData',xd);
%             set(hand,'YData',yd);            
%         end 
%     end
    h=text(1.5,-1.9,sprintf('t:%0.1f',tspan(i+1)));
    cF=getframe(gcf); % current frame
    writeVideo(vo,cF);
    delete(h); % clear text, otherwise overtyping produces illegibility
end
% % 
close(vo); % close video object/file