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

clear;
clf;
cla;
hold off;


n=100;
t_max = 10;
time_steps=1000;


% GENERAL SETUP

% Generating the boundary of the region to track
circle=@(s)([1.6+0.4.*cos(s);0.4.*sin(s)]); % parametrised circle
spoints=linspace(-pi,pi,n); % points in s-parameter space
pppoints=circle(spoints); % points in phase-plane

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

% SYMPLECTIC SOLVER

% Setup
pstr='';
for i =1:n
    pstr=[ pstr sprintf('-sin(q(%u));',i ) ];
end
pstr = [ '[' pstr ']' ];

dpdt=inline(pstr,'q');

qstr='';
for i =1:n
    qstr=[ qstr sprintf('p(%u);',i ) ];
end
qstr = [ '[' qstr ']' ];

dqdt=inline(qstr,'p');

% Solve
q0=pppoints(1,:); % q=y is horizontal/x-axis
p0=pppoints(2,:); % p=y' is vertical/y-axis

[symp.q,symp.p]=seep(dqdt,dpdt,tspan,q0,p0); % seep is in the package [1]

vo = VideoWriter('blob.avi','Uncompressed AVI'); % video object
%Settings
vo.FrameRate=30;
%vo.Quality=100;
open(vo); % open up the video file

for i=1:time_steps
   % NON-SYMPLECTIC
   % q=y is horizonal/x-axis, p=y' is vertical/y-axis
   
   % SYMPLECTIC
   % q=y is horizonal/x-axis, p=y' is vertical/y-axis
   symp.xvert=symp.q(i,:);
   symp.yvert=symp.p(i,:);
  
   % PLOTTING
   plot(symp.xvert,symp.yvert,'-b');
   set(gca,'Color',[1 1 1]);
   title('Phase space evolution; y''''+sin(y) = 0')
   xlabel(gca,'q (= y)');
   ylabel('p (= y'')');
   axis square;
   axis([-2 2 -2 2]);
   text('Position',[-1,-1.8], 'String','Symplectic Euler','Color',[0 0 1]);
   text(1.5,-1.9,sprintf('t:%0.3g',tspan(i)));
   cF=getframe(gcf); % current frame
   writeVideo(vo,cF);
end
close(vo);