clear;
%
%   Fractal Fern
%
%   the transformations:
%     
%          0  0     0            0.2   -0.26  0
%     t1=  0  0.16  0       t2=  0.23   0.22  1.6      t3=  ...
%          0  0     1            0      0     1
%
t1=[ 0     0    0;  0     0.16  0;    0 0 1]; 
t2=[ 0.2  -0.26 0;  0.23  0.22  1.6;  0 0 1]; 
t3=[-0.15  0.28 0;  0.26  0.24  0.44; 0 0 1]; 
t4=[ 0.85  0.04 0; -0.04  0.85  1.6;  0 0 1]; 

points=100000;              % number of points in the fern
ran_num=rand(points,1);     % find points random numbers
xyf=zeros(2,points);        % initialize memory for the fern

xy=[.5; .4; 1];             % a row vector as initial condition
%
%   xy is most likely not a point on the attractor (the fern)
%   perform 1000 iteration as transient to reach the attractor   
%
ran_trans=rand(1000,1);     % find 1000 random numbers
for k=1:1000
    if ran_trans(k)<=0.01 xy=t1*xy; end;
    if ran_trans(k)>0.01 & ran_trans(k)<=0.08 xy=t2*xy; end;
    if ran_trans(k)>0.08 & ran_trans(k)<=0.15 xy=t3*xy; end;
    if ran_trans(k)>0.15 xy=t4*xy; end;
end
%
%   Now calculate the points of the fern
%
for k=1:points
    if ran_num(k)<=0.01 xy=t1*xy; end;
    if ran_num(k)>0.01 & ran_num(k)<=0.08 xy=t2*xy; end;
    if ran_num(k)>0.08 & ran_num(k)<=0.15 xy=t3*xy; end;
    if ran_num(k)>0.15 xy=t4*xy; end;
    xyf(1:2,k)=xy(1:2);
end
figure('Position',[200 200 600 600]);       % create a figure
scatter(xyf(1,:),xyf(2,:),2,'g','filled');  % plot points size 2 green
set(gca,'Color','k');                       % axis background black    
    
    
    
    