% bakers_map.m
% Jeremie Korta
% November 19, 2003
% The following program simulates the baker's map in Ott (page 77 2nd ed)
% parametrized by splitting parameters, alpha and beta as well as
% stretching parameters lambda_a and lambda_b.  Several features of the map
% may be explored, such as the strange attractor, which any irrational
% initial condition between 0 and 1 will eventually explore entirely, as
% well as divergence of ICs and map of a circle or the unit square itself. 

figure(1)
clf        % Uses Figure 1... this clears it first.

% splitting parameter:
alpha = 0.5000000001; % small pertubation needed for plotting attractor
beta = 1-alpha;

% stretching parameters:
lambda_a = 1/3;
lambda_b = 1/3;

% map random points (c = 0), a filled circle (c = 1)
% generate a picture of the strange attractor (c = 2)
% sensitivity of initial conditions (c = 3)
% or just follow the mapping of a single point (c = 4):
c = 4;

if c == 2
    z = [1 1]/sqrt(2); % initial (irrational!) position
    N = 20; % remove transient behavior
    for i = 1:N
        z(1) = (z(2) < alpha).*(lambda_a*z(1)) + (z(2) >= alpha).*(lambda_b*z(1)+(1-lambda_b));
        z(2) = (1/alpha)*(z(2) < alpha).*z(2) + (1/beta)*(z(2) >= alpha).*(z(2)-alpha);
    end
    iter = 10000;  % keep track of next few iterates
    x = zeros(iter,2);
    x(1,:) = z;
    for i = 2:iter
        x(i,1) = (x(i-1,2) < alpha).*(lambda_a*x(i-1,1)) + (x(i-1,2) >= alpha).*(lambda_b*x(i-1,1)+(1-lambda_b));
        x(i,2) = (1/alpha)*(x(i-1,2) < alpha).*x(i-1,2) + (1/beta)*(x(i-1,2) >= alpha).*(x(i-1,2)-alpha);    
    end
    plot(x(:,1),x(:,2),'.b','MarkerSize',3)
    s=['Baker''s Map Invariant Set with \alpha = ',num2str(alpha),', \lambda_a = ', ...
            num2str(lambda_a),' and \lambda_b = ', num2str(lambda_b)];
    title(s) 
    
elseif c == 0
    % number of mappings:
    n = 4;
    x = zeros(n,100000);
    y = zeros(n,100000);
    x(1,:) = rand(1,100000);
    y(1,:) = rand(1,100000);
    for i = 2:n
        x(i,:) = (y(i-1,:) < alpha).*(lambda_a*x(i-1,:)) + (y(i-1,:) >= alpha).*(lambda_b*x(i-1,:)+(1-lambda_b));
        y(i,:) = (1/alpha)*(y(i-1,:) < alpha).*y(i-1,:) + (1/beta)*(y(i-1,:) >= alpha).*(y(i-1,:)-alpha);
    end
    for i=1:n
        subplot(2,2,i); plot(x(i,:),y(i,:),'.b','MarkerSize',3)
        xlim([0, 1])
        ylim([0, 1])
    end
    
elseif c == 1
    n = 16;
    x = zeros(n,10000);
    y = zeros(n,10000);
    % circle radius:
    r = 0.01;
    % circle position:
    z = [0.6325, 0.619];
    % uniformly distributed vector of radii:
    v = r*rand(1,10000);
    % uniformly distributed vector of angles:
    angle = 2*pi*rand(1,10000);
    x(1,:) = z(1) + v.*cos(angle);
    y(1,:) = z(2) + v.*sin(angle);   
    for i = 2:n
        x(i,:) = (y(i-1,:) < alpha).*(lambda_a*x(i-1,:)) + (y(i-1,:) >= alpha).*(lambda_b*x(i-1,:)+(1-lambda_b));
        y(i,:) = (1/alpha)*(y(i-1,:) < alpha).*y(i-1,:) + (1/beta)*(y(i-1,:) >= alpha).*(y(i-1,:)-alpha);
    end
    for i=1:n
        subplot(4,4,i); plot(x(i,:),y(i,:),'.b','MarkerSize',3)
        xlim([0, 1])
        ylim([0, 1])
    end
    
elseif c == 3
    n = 16;
    x = zeros(n,2);
    y = zeros(n,2);
    x(1,:) = [1 1.01]/sqrt(2);
    y(1,:) = [1 1.01]/sqrt(2);
    for i = 2:n
        x(i,:) = (y(i-1,:) < alpha).*(lambda_a*x(i-1,:)) + (y(i-1,:) >= alpha).*(lambda_b*x(i-1,:)+(1-lambda_b));
        y(i,:) = (1/alpha)*(y(i-1,:) < alpha).*y(i-1,:) + (1/beta)*(y(i-1,:) >= alpha).*(y(i-1,:)-alpha);
    end
    for i=1:n
        subplot(4,4,i)
        plot(x(i,:),y(i,:),'.r','MarkerSize',7)
        xlim([0, 1])
        ylim([0, 1])
    end
    figure(2)
    z = sqrt((x(:,2)-x(:,1)).^2 + (y(:,2)-y(:,1)).^2);
    plot(z);
    
elseif c == 4
    n = 16;  % shouldn't carry out calculation much more, otherwise numerical
             % truncation becomes significant
    x = zeros(n,1);
    y = zeros(n,1);
    %x(1,1) = 3/4; % a period-2 orbit (on the attractor)
    %y(1,1) = 1/3;
    x(1,1) = 1/26; % a period-3 orbit
    y(1,1) = 4/7;
    for i = 2:n
        x(i,1) = (y(i-1,1) < alpha).*(lambda_a*x(i-1,1)) + (y(i-1,1) >= alpha).*(lambda_b*x(i-1,1)+(1-lambda_b));
        y(i,1) = (1/alpha)*(y(i-1,1) < alpha).*y(i-1,1) + (1/beta)*(y(i-1,1) >= alpha).*(y(i-1,1)-alpha);
    end
    for i=1:n
        subplot(4,4,i)
        plot(x(i,1),y(i,1),'.r','MarkerSize',10)
        xlim([0, 1])
        ylim([0, 1])
    end
end