% based on form appearing in Ott (chaos in dynamical systems) page
% 216.
% theta(n+1)= mod(theta(n)+p(n),2*pi);
% p(n+1)    = p(n)+K*sin(theta(n+1));
% Eli Tziperman, for course on nonlinear dynamics and chaos, Fall 2003.

close all; clear
% nonlinearity:
K=0.75;
% number of initial conditions to cover:
N_theta=10;
N_p=10;
% number of iterations per each initial condition:
N=100;
% plotting parameter:
MarkerSize=1;

% Initialize vecctors of theta and p:
theta=zeros(N+1);
p    =zeros(N+1);

% Start from many different initial conditions to create a phase
% space picture:
for theta0=0:2*pi/N_theta:2*pi
  for p0=0:2*pi/N_p:2*pi
    theta(1)=theta0;
    p(1)=p0;
    % plot initial conditions as thicker dots:
    h=plot(theta0,p0,'.'); set(h,'MarkerSize',6); hold on;
    % limit plot to only a part of the p-axis.  extend this hy using
    % "axis([0 2*pi -3*pi 5*pi]);" to see that trajectories are not
    % trapped to this part only:
    axis([0 2*pi 0 2*pi]);

    % start from these i.c., iterate N times and plot all iterates:
    for n=1:N
      theta(n+1)= mod(theta(n)+p(n),2*pi);
      p(n+1)    = p(n)+K*sin(theta(n+1));
    end
    h=plot(theta,p,'.'); set(h,'MarkerSize',MarkerSize);
  end
end

xlabel('theta');
ylabel('p');
title('standard map, thicker dots mark initial conditions');