clear
a2=-0.733333;
a1bar=0.004594;
omega=45;
V_max=5.006;
omega_c=22;
phi_0=0;
psi_0=9.1743;
phi_old=phi_0;
psi_old=psi_0;
phi_e=.7173;
psi_e=2*pi;
phi=1:150;
psi=1:150;
P=[8.6394 4.6968;4.6968 5.9586];
   for t=1:150
x=[phi_old-phi_e psi_old-psi_e]';
V=x'*P*x;
if V>V_max
omega=45;
else
omega=omega_c;
end
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)*(-rand);
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
plot(phi,psi/(2*pi),'.')
plot(phi_c,eval(psi_upper))
plot(phi_c,eval(psi_lower))
plot(phi_e,psi_e,'+')
plo_lev
%grid;
title('An orbit of the controlled ball map','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)
hold off
figure
plot([0 1:1:150],[phi_0 phi(1:150)],'.')
grid;
title('A plot of the phase angle versus time of the controlled ball map','Fontsize',12);
xlabel('Iteration','Fontsize',12)
ylabel('Phase angle of plate at time of bounce','Fontsize',12)
set(gca,'Fontsize',12)
figure
plot([0 1:1:150],[psi_0/(2*pi) psi(1:150)/(2*pi)],'.')
grid;
title('A plot of the phase change versus time of the controlled ball map','Fontsize',12);
xlabel('Iteration','Fontsize',12);
ylabel('Phase change between bounces/(2*pi)','Fontsize',12)
set(gca,'Fontsize',12)