clear
a2=-0.733333;
a1bar=0.004594;
omega=45;
phi_old=0;
psi_old=9.1743;
phi=[1:10000];
psi=[1:10000];
   for t=1:10000
phi(t)=phi_old+psi_old;
      while phi(t)>2*pi
phi(t)=phi(t)-2*pi;
       end
phi_old=phi(t);
psi(t)=-a2*psi_old+a1bar*(omega^2)*cos(phi_old);
    if psi(t)<0
psi(t)=-psi(t);
    end
psi_old=psi(t);
end
phi_c= 0:.001:2*pi;
   psi_upper='(2*a1bar*(45^2)+a1bar*(45^2)*cos(phi_c))/(2*pi)';
   dummy=cos(phi_c);
   if dummy <0
      dummy=0;
   end
   psi_lower='a1bar*45^2*dummy/(2*pi)';
   axis([0 2*pi 0 5])
   hold on
%   pause
   plot(phi,psi/(2*pi),'.')
   plot(phi_c,eval(psi_upper),phi_c,eval(psi_lower))
grid;
title('An orbit of the ball map (11.117)','Fontsize',12)
xlabel('Phase angle of plate at time of bounce','Fontsize',12)
ylabel('Phase change between bounces/(2*pi)','Fontsize',12)
set(gca,'Fontsize',12)
