
% Iterated function system
% Ricardo Carretero Mar/03

IFS = inline('[a*x+b*y+e c*x+d*y+f]','x','y','a','b','e','c','d','f');
%IFS = inline('[a b;c d]*[x;y] + [e;f]','x','y','a','b','e','c','d','f');

a = [0.00  0.20 -0.15  0.75];
b = [0.00 -0.26  0.28  0.04];
c = [0.00  0.23  0.26 -0.04];
d = [0.16  0.22  0.24  0.85];

e = [0.00  0.00  0.00  0.00];
f = [0.00  1.60  0.44  1.60];

p = cumsum([0.10 0.08 0.08 0.74]);

Ntrans=100;
Npts=1000;
Nit=10;
NitBox=3;
x=pi/30;		% ICs
y=x*pi/4;
ms=1;
lw=2;
fs=15;
ax=[-3.3 3.3 -0.22 10.8];

Bx=[-2.5,2.95];
By=[0,10.5];
Bxx=[Bx(1) Bx(2) Bx(2) Bx(1) Bx(1)];
Byy=[By(1) By(1) By(2) By(2) By(1)];
sizeB=length(Bxx);

figure(1)
clf
axis(ax);
set(gca,'FontSize',[fs]);
xlabel('x_n'); ylabel('y_n');

Bold=zeros(4^NitBox,sizeB,2);
Bnew=zeros(4^NitBox,sizeB,2);

Bold(1,:,1) = Bxx;
Bold(1,:,2) = Byy;
colvec=jet(NitBox*10);

box=1;
if(box)
 plot(Bold(1,:,1),Bold(1,:,2),'LineWidth',lw,'Color',colvec(5,:))
 axis(ax);
 hold on
 fprintf('Press any key to continue...\n'); pause
 for(jj=1:NitBox-1)
  Nbold=4^(jj-1);
  Nbnew=4^jj;
  for(kk=1:Nbold)
   for(ii=1:4)
   col=(jj)*10+3+ii;
    BB=IFS(Bold(kk,:,1)',Bold(kk,:,2)',a(ii),b(ii),e(ii),c(ii),d(ii),f(ii));
    plot(BB(:,1),BB(:,2),'r','LineWidth',lw,'Color',colvec(col,:))
    axis(ax);
    xb=mean(BB(:,1));
    yb=mean(BB(:,2));
    text(xb,yb,[num2str(ii)],'FontSize',[18-4*jj],...
         'FontWeight','bold','Color',colvec(col,:))
    Bnew((kk-1)*4+ii,:,1)=BB(:,1);
    Bnew((kk-1)*4+ii,:,2)=BB(:,2);
    figure(1)
    fprintf('Press any key to continue...\n'); pause
   end
  end
  Bold=Bnew;
 end
 axis(ax);
end

p1 = plot(x,y,'.','EraseMode','none','MarkerSize',ms,'Color',[0 0.6 0]);
axis(ax);

X=[];
Y=[];

for i=1:Ntrans
  [M,j]=max(rand(1)<p);
  F = IFS(x,y,a(j),b(j),e(j),c(j),d(j),f(j));
  x=F(1);
  y=F(2);
end

for i=1:Nit
 X1=[];
 Y1=[];
 for i=1:Npts
  r=rand(1);
  [M,j]=max(rand(1)<p);
  F = IFS(x,y,a(j),b(j),e(j),c(j),d(j),f(j));
  x=F(1);
  y=F(2);
  X1=[X1,x];
  Y1=[Y1,y];
 end
 set(p1,'Xdata',X1,'YData',Y1)
 drawnow;
 X=[X,X1];
 Y=[Y,Y1];
end

plot(X,Y,'.','EraseMode','none','MarkerSize',ms,'Color',[0 0.6 0]);
axis(ax);
%hold on
%plot(X(1:10),Y(1:10),'b*','EraseMode','none','MarkerSize',8);
drawnow;

