function   [X,h] = ifs(A, p, N, opt1)

%
% [X,h] = ifs(A, p, N, opt1)
% Chaos game for 2-D Iterated Function Systems
% Inputs:
% A = matrix of coefficients (see examples in ex_ifs.m)
% p = vector of probabilities
% N = number of iterations
% opt1 = drawing options ('yes', 'no')
% Outputs:
% X = trajectory matrix 
% h = plot handle
% Author:
% Marco Sandri - Verona - Italy - Email: sandri.marco@gmail.com
%


if nargin < 1
    error('Requires at least one input argument: matrix A');
end
if nargin < 2
    m = size(A,1);
    p = ones(1,m)/m;
end
if nargin < 3
    N = 1000;
end
if nargin < 4
    opt1='yes';
end


p = p(:)';
[m1,col]=size(A);
m2 = length(p);

%%%%%%%%%%%%%%%%%%%%%%%
%%% Error messages
%%%%%%%%%%%%%%%%%%%%%%%
if col~=6
    error('The number of columns of A must be 6');
end;
if m1~=m2
    error('The number of rows of A and the length of p must be equal');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% If sum(p) not equale 1, then normalize p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if sum(p)~=1 
    p = p/sum(p);
end;
p = cumsum(p);


X = zeros(N,2);
X(1,:)=rand(1,2);

trans=100;
for k=1:N+trans-1
    n = rand(1);
    n_mappa = find(diff([0 (p > n)]));
    T = reshape(A(n_mappa,:),3,2)';
    X(k+1,:)=(T*[X(k,:) 1]')';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Graphics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch(opt1)
case('yes')
    h=plot(X(trans+1:end,1),X(trans+1:end,2),'.',...
        'Markersize',3,'Color',[0 1 0]);
    set(gcf,'Color',[0 0 0]);
    set(gca,'Color',[0 0 0]);
end