clear all;
g=9.81;    
l=1.0;   %pendulum length
A=0.01;     %shake amplitude m
oc=sqrt(g/A);
omega=oc*0.1;   %shaking 
                %omega=2*pi/100;   %shaking 
makeslowmovie=1;
makefastmovie=0;

dt=0.01;tmax=12;
n=tmax/(dt)+1;
tspan=linspace(0,tmax,n);
y0 = [pi/10 0.001];  %initial angle and thetadot
[t,y] = ode45(@(t,y) odefcny(t,y,g,l,A,omega), tspan, y0);

s=size(y)
tlen=s(1)
figure(1);clf;
subplot(2,1,1);
plot(t,y(:,1));grid on;
ylabel('\theta');
xlabel('t');
title('\omega<<\omega_c');

if (makeslowmovie==1)
figure(2);clf;
%vidfile = VideoWriter('./figs/testmovie.mp4','MPEG-4');
%vidfile = VideoWriter('./figs/Slowshake.avi','Uncompressed AVI');
vidfile = VideoWriter('./figs/Slowshake.avi');
vidlen=20; %in seconds
vidfile.FrameRate = floor(tlen/vidlen)
open(vidfile);
theta=linspace(0,2*pi,100);
xr=l*cos(theta);yr=l*sin(theta);
for in=1:tlen
yp=l*cos(y(in,1));xp=l*sin(y(in,1));
plot(xr,yr,'k--');
hold on;
plot(xp,yp,'.r','MarkerSize',55);grid on;axis('square');
hold off;

Fr = getframe(gcf); 
writeVideo(vidfile,Fr);
end
close(vidfile);
end %if (makeslowmovie==1)


figure(1);
omega=oc*100;   %shaking 
[t,y] = ode45(@(t,y) odefcny(t,y,g,l,A,omega), tspan, y0);
subplot(2,1,2);
plot(t,y(:,1));grid on;
ylabel('\theta');xlabel('t');
title('\omega>>\omega_c');

print('-dpng','./Vertical_shake.png')
s=size(y)
tlen=s(1)

if (makefastmovie==1)
figure(2);clf;
%vidfile = VideoWriter('./figs/testmovie.mp4','MPEG-4');
%vidfile = VideoWriter('./figs/Fastshake.avi','Uncompressed AVI');
vidfile = VideoWriter('./figs/Fastshake.avi');
vidlen=20; %in seconds
vidfile.FrameRate = floor(tlen/vidlen)
open(vidfile);
theta=linspace(0,2*pi,100);
xr=l*cos(theta);yr=l*sin(theta);
for in=1:tlen
yp=l*cos(y(in,1));xp=l*sin(y(in,1));
plot(xr,yr,'k--');
hold on;
plot(xp,yp,'.r','MarkerSize',55);grid on;axis('square');
hold off;

Fr = getframe(gcf); 
writeVideo(vidfile,Fr);
end
close(vidfile);
end %if (makeslowmovie==1)