function fig7_13

%  using velocity Verlet on the system y' = f(t,y)  with   y(0) = y0

t0=0;
y0=[pi/4 0];
tmax=4000;
k=1/2;
n=round(tmax/k)
t=linspace(t0,tmax,n);

[y,v]=verlet(t,y0,k,n);

% determine the number of zero crossings
% k=1/2 => 1009 crossing for 0<t<3266
% nzeros=0;
% for it=2:n
%     if y(it)*y(it-1)<0
%         nzeros=nzeros+1;
%     end
% end
%
% nzeros_verlet=nzeros

figure(1)
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.035,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

plot(t,y,'-b','LineWidth',2)
hold on
T1=40;
axis([0 T1 -0.85 0.85])
plot([0 T1],[pi/4  pi/4],'--r','LineWidth',2)
set(gca,'YTick',[-pi/4 0 pi/4]);
set(gca,'YTickLabel',{' ','0',' '})
say=['-\pi/4'];
text(-3,-pi/4,say,'FontSize',20)
say=['\pi/4'];
text(-2.6,pi/4,say,'FontSize',20)
set(gca,'XTick',[0 10 20 30 40])
box on
grid on
xlabel('t-axis')
ylabel('Theta')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

exportgraphics(gca,'/Users/mark/Desktop/ver1.eps')

figure(2)
clf
% get(gcf)
set(gcf,'Position', [24 796 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',-0.01,'MB',0.18,'MR',-0.001,'ML',0.03,'P',0.04)

tr=linspace(0,tmax,100);
aa=1e-4*log(7*pi/3);
for ir=1:100
    R(ir)=0.25*pi*exp(-aa*tr(ir));
end

plot(t,y,'-b','LineWidth',0.01,'color','#00BFFF')
hold on
plot([0 tmax],[pi/4  pi/4],'--r','LineWidth',2)
set(gca,'YTick',[-pi/4 0 pi/4]);
% set(gca,'YTickLabel',{'-pi/4','0','pi/4'})
set(gca,'YTickLabel',{' ','0',' '})
say=['-\pi/4'];
text(-300,-pi/4,say,'FontSize',20)
say=['\pi/4'];
text(-250,pi/4,say,'FontSize',20)
box on
grid on
xlabel('t-axis')
ylabel('Theta')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

%exportgraphics(gca,'/Users/mark/Desktop/ver2.eps')

% right-hand side of DE
function z=f(t,y)
z=[y(2) -sin(y(1))];

% Verlet method
function [ypoints, vpoints] =verlet(t,y0,k,n)
ypoints=y0(1);
vpoints=y0(2);
y=y0(1);  v=y0(2);
for i=2:n
    Y=y+k*v-0.5*k^2*sin(y);
    V=v-0.5*k*sin(y)-0.5*k*sin(Y);
    ypoints=[ypoints, Y];
    vpoints=[vpoints, V];
    y=Y;  v=V;
end











