function IterateSineMapAttSensi
% ======================================================================= %
% Compute and plot the iterates, x_{n+1} = S(x_n),
% where                                           S(x) = c*pi*sin(x).
% In particular, check the sensitivity of the attractor.
%
% Check the help for: IterateSineMap
% Main difference here it that this replaces the "periodicity check"
% by a sensitivity calculation.
% AGAIN:
% 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.
%     For the purposes of this scrip, it is a good idea of take n1 very
%     large, so convergence to the attractor is guaranteed even when
%     very close to a bifurcation point [when convergence can be rather
%     slow].
%     I found that, for the c's that correspond to period 2^n, with n
%     "large" [bigger than 7] I had to take n1 huge. This may be because
%     the transitions here are very closely spaced, so convergence to the
%     attractor is rarely fast.
% (2) After this is continues the iteration, with nz = max(n2, n3) more
%     ietrates. But it does TWO iterations in parallel, as follows:
%     --- Start from x_{n1},    and store the iterates in the array x.
%     --- Start from x_{n1}+xp, and store the iterates in the array y.
%     Here xp is a small perturbation.
% (3) Finally, it plots (a) x versus n for 1 leq n leq n3;
%     and               (b) x-y versus n for 1 leq n leq n3;
%     (a) provides a picture of the attractor solution;
%     (b) shows what happens with the perturbation; does it grow? How
%         fast? How much?
%     n2 and n3 are parameters to play with.
%     To actually "see" the attractor, n2 cannot be too large.
%     To see the behavior of the perturbation, n3 may need to be large.
%
%     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.
% ======================================================================= %
c_store(1)  = 0.7400000;%    Period     2 error decays to 0
c_store(2)  = 0.8400000;%    Period     4 error decays to 0
c_store(3)  = 0.8590000;%    Period     8 error decays to 0
c_store(4)  = 0.8643000;%    Period    16 error decays to 0
c_store(5)  = 0.8655700;%    Period   128 error decays to 0
% --------------------- %
%                         Theoretically the next 3 corresponds to
%                         periodic solutions. But the periods too
%                         long to trust a calculation with ~10^{-6}
%                         precision. This is in the "grey" zone fixed
%                         point calculations always yield.
%                         The errors here sometimes decay, and other
%                         times do not, but they stay below 10^{-4}.
% --------------------- %
c_store(6)  = 0.8655790;%    Period  1024
c_store(7)  = 0.86557924;%   Period  2048
c_store(8)  = 0.8655792686;% Period 16384
c_store(9)  = 0.8655900;%    Maybe chaos. Error grows to ~ 10^{-2}
c_store(10) = 0.8660000;%    Maybe chaos. Error grows to ~ 10^{-1}
%
k = 5;  c = c_store(k);
%
% --------------- Select parameter in the map.
x0 = pi*rand; %   This makes the initial data random in [0, pi].
% --------------- Select initial data for iteration.
n1 = 1000000;
% --------------- Iterations used to approach attractor
n2 = 10000;
% --------------- Iterations after n1 used for attractor plot.
n3 = 1000000;
% --------------- Iterations after n1 used for attractor sensitivity.
xp = 1e-5;
% --------------- Perturbation to attractor at n1.
%
% ======================================================================= %
% Calculation below.
nz = max(n2, n3);
x  = zeros(1, nz);
y  = zeros(1, nz);
%
for n=1:n1; x0 = c*pi*sin(x0); end
y0 = x0+xp;
for n=1:nz;
  x0 = c*pi*sin(x0); x(n) = x0;
  y0 = c*pi*sin(y0); y(n) = y0;
end
y = y-x;
%
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 c = ', num2str(c, 8), ' attractor. }'], '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
%
figure
whitebg(gcf, 'w')
plot(y(1:n3), '-b', 'LineWidth', 3); hold on
xlabel('{\bf n}', 'FontSize', 18)
ylabel('{\bf Perturbation to attractor. }', 'FontSize', 18)
title(['{\bf c = ', num2str(c, 8), ' sine map. }'], 'FontSize', 18)
set(gca, 'FontWeight', 'bold', 'FontSize', 16)
miny = min(y(1:n3));
maxy = max(y(1:n3));
axis([1, n3, 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