
global k m L0 L g
format compact

% paramaters:

tf=200;
k=1;
m=0.5;
g=1;
L0=1;
L=m*g/k+L0

% ICs:

x0 = 1, y0 = L/2, xp0 = 0; yp0 = 0;   % regular
x0 = 0.02, y0 = L/1.15, xp0 = 0; yp0 = 0; % chaotic
x0 = 1, y0 = 4*L, xp0 = 0; yp0 = 0;   % regular high energy

K0 = 0.5*(xp0.^2+yp0.^2);
V0 = ext_pend_potential(x0,y0);
E0 = K0 + V0;

tspan=tf;
Y0 = [x0;y0;xp0;yp0];

options= odeset('Reltol',1.e-2,'Abstol',1.e-3,'MaxStep',0.1);
%options= odeset('Reltol',1.e-8,'Abstol',1.e-9,'MaxStep',0.01);
%options=[odeset('OutputFcn',@odephas2)];
%options=[];

[T,Y]=ode45('ext_pend_ODEs',tspan,Y0,options);

Torbit = length(T);
x=Y(:,1);
y=Y(:,2);
xp=Y(:,3);
yp=Y(:,4);

K = 0.5*(xp.^2+yp.^2);
V = ext_pend_potential(x,y);
E = K + V;

fs=14;

figure(1);clf
subplot(3,1,1);set(gca,'FontSize',fs);
plot(T,x);
axis tight
xlabel('t');ylabel('x(t)');
subplot(3,1,2);set(gca,'FontSize',fs);
plot(T,y)
axis tight
xlabel('t');ylabel('y(t)');
subplot(3,2,5);set(gca,'FontSize',fs);
plot(x,y);
axis tight
xlabel('x(t)');ylabel('y(t)');
subplot(3,2,6);set(gca,'FontSize',fs);
plot(xp,yp)
axis tight
xlabel('x''(t)');ylabel('y''(t)');

monitor=0;
if(monitor)
 figure(2);clf
 subplot(3,1,1);set(gca,'FontSize',fs);
 plot(T,E-E0)
 axis tight
 ylabel('E(t)-E0');
 subplot(3,1,2);set(gca,'FontSize',fs);
 plot(T,K,T,K0*ones(Torbit,1),'r')
 axis tight
 ylabel('K(t)');
 subplot(3,1,3);set(gca,'FontSize',fs);
 plot(T,V,T,V0*ones(Torbit,1),'r')
 axis tight
 xlabel('t');ylabel('V(t)');
end

orbit=1;
mymovie=1;
skip=10;
Npts=20;

if(mymovie)
 figure(3);clf
 mx=min(x);
 Mx=max(x);
 my=min(y);
 My=max(y);
 pendulum_shot(x(1),y(1),Npts,1.2*[mx Mx my My]);
 fprintf('Press any key to see movie');
 pause;
 for ii=1:skip:length(T);
  xx=x(ii);
  yy=y(ii);
  pendulum_shot(xx,yy,Npts,1.2*[mx Mx my My]);
  if(orbit)
   hold on
   plot(x(1:skip:ii),y(1:skip:ii),'-k')
   drawnow
   hold off
  end
 end
end
