function IterateSineMap
% ======================================================================= %
% Compute and plot the iterates, x_{n+1} = S(x_n),
% where                                           S(x) = c*pi*sin(x).
% Note that S(x) maps the interval [0, pi] onto
% itself, provided 0 <= c <= 1.
% THE RELEVANT PARAMETERS ARE AT THE START OF
% THE SCRIPT, WHERE THEY ARE EASY TO CHANGE.
%
% Specifically:
% (1) The script starts with some specified x_0, and computes x_n up to
%     x_{n1}. It does this to get over the "transient" period, and reach
%     its attractor behavior.
% (2) Exactly what n1 is is something you need to play with. It should
%     not be small, but in many cases it does not need to be very large,
%     say: something like n1 = 300. HOWEVER, near the critical values of
%     c where the behavior changes (bifurcations) the convergence to the
%     attractor can be quite slow (we will see why this happens later).
%     Because of this,
%
%         I RECOMMEND that you take n1 HUGE. With modern computers you
%         will not be able to tell the difference between n1 = 300, and
%         n1 = 3000000. And the second choice will avoid puzzles [i.e.:
%         a change in initial data gives you a different answer, and
%         neither looks recognizable as a special behavior].
%         
% (3) Then it continues computing the x_n, from x_{n1+1} to x_{n1+n2}.
%     However it now saves the x_n, and at the end it makes a plot of
%     x_n versus n [i.e.: a line connecting the points (n, x_n) with
%     straight segments]. This line is SOLID RED.
% (4) Again, n2 is a parameter that you need to play with. It has to be
%     large enough that you can see what the sequence of x_n is doing
%     [e.g.: it is periodic with some period, or looks random, or ...].
%     But it cannot be too large, because then you will not be able to
%     see anything in the plot!
%     Note: within MatLab you can use the zoom feature to see details,
%           but this will not work once you save the plot as a pdf for
%           the written answer to the problem set.
%
% FINALLY:
% Because in many cases the behavior is periodic (or close to), the code
% has a feature that allows you to detect this.  This is done via the
% parameter n3. If you specify a value for n3 that is bigger than zero
% (say 1, 2, ...), then the script will also plot (n, x(n+n3}-x{n}), as
% a SOLID BLUE line (in a different window). If the behavior is periodic
% of period n3 [i.e.: x_{n+n3} = x_n], then this second line will be 0
% everywhere.
% FURTHER: you may see this blue line being small somewhere only,
% meaning that you have near periodic behavior there but not elsewhere.
% This is what characterizes "intermittent chaos". <=====================
%
% REMARKS:
% 1) n3 = 1 [period one] corresponds to x_n constant.
% 2) If x_n is periodic of, say period 3, it will also be periodic of
%    period 6. Thus, to be sure the period is 6, you have to rule out
%    periods 2, and 3. So you will need to check n3 = 1, 2, 3, 6 ...
%    though many times it is obvious it is not constant [so no need for
%    n3 =1] and it may also be obvious that, say, it is not period 3.
%    Then you only check n3 = 2, 6.
% ======================================================================= %
c  = 0.94;
x0 = pi*rand; % this makes the initial data random in [0, pi]
n1 = 3000000;
n2 = 30;
n3 = 3;
%
% ======================================================================= %
% Calculation below.
x  = zeros(1, n2+n3);
%
for n=1:n1; x0 = c*pi*sin(x0); end
for n=1:n2+n3; x0 = c*pi*sin(x0); x(n) = x0; end
%
figure
whitebg(gcf, 'w')
plot(x(1:n2), '-r', 'LineWidth', 3); hold on
plot(x(1:n2), 'ob', 'MarkerSize', 8, 'LineWidth', 2);
xlabel('{\bf n}', 'FontSize', 18)
ylabel('{\bf x}', 'FontSize', 18)
title(['{\bf Iterates for c = ', num2str(c, 8), ' }'], 'FontSize', 18)
set(gca, 'FontWeight', 'bold', 'FontSize', 12)
minx = floor(1000*min(x))/1000;
maxx = ceil(1000*max(x))/1000;
axis([1, n2, minx-0.0005, maxx+0.0005])
grid on; zoom on
%
if n3 > 0
   figure
   whitebg(gcf, 'w')
   y = x(1+n3:n2+n3)-x(1:n2);
   plot(y, '-b', 'LineWidth', 3); hold on
   xlabel('{\bf n}', 'FontSize', 18)
   ylabel(['{\bf shift by ', num2str(n3), ' difference}'], 'FontSize', 18)
   title(['{\bf Iterates for c = ', num2str(c, 8), ' }'], 'FontSize', 18)
   set(gca, 'FontWeight', 'bold', 'FontSize', 12)
   miny = min(y);
   maxy = max(y);
   % axis([1, n2, miny-0.0005, maxy+0.0005])
   axis([1, n2, miny-1e-100, maxy+1e-100])
   grid on; zoom on
   P = get(gcf, 'Position');
   set(gcf, 'Position', P + [0+0.25*P(3), 0-0.25*P(4), 0, 0])
end
% END