function fern_simulation(N,MS,Kolor,Mrkr)

% This routine draws Barnsley's fern
% The transformations for the fern are each in the form
% T(x,y = (a*x+b*y+e, c*x+d*y+f). If the probability that a
% particular transformation T is selected is p, you can characterize
% each transformation by providing the values of
% a,b,c,d,e,f, and p, in that order.
%
% T1 -> 0,0,0,.16,0,0,.01
% T2-> .85,.04,-.04,.85,0,1.6,.85
% T3-> .2,-.26,.23,.22,0,1.6,.07
% T4-> -.15,.28,.26,.24,0,.44,.07
%
% INPUT
%	N -> The number of iterations.

% Close all existing figure windows. Delete this line
% if you do not want this action.
close all

% Initialize figure window and axes
fig1 = figure;
set(fig1,'Color','k','Position',[50   50   560   420])
ax1 = axes;


% Initialize some space to hold the sequence of points
% X1, X2, ..., XN
%X=zeros(N,2);

% The initial value X1 = [0.5,0.5].
X = [0.5,0.5];
plot(X(:,1),X(:,2),'Marker',Mrkr,'markersize',MS, ...
    'MarkerFaceColor','k','MarkerEdgeColor','k')
hold on
set(gca,'DataAspectRatio',[1 1 1],'Xlim',[-5 5], ...
    'Ylim',[0 10])
axis off

% The main loop where the iterations are performed.
k = 0;
Nc = 128;
K = colormap(cool(Nc));

while k < N
    drawnow
    %pause
    k = k + 1;
    r = rand(1,1)*100;
    if r < 1
        X = T(X,0,0,0,.16,0,0);
    elseif r < 86
        X = T(X,.85,.04,-.04,.85,0,1.6);
    elseif r < 93
        X = T(X,.2,-.26,.23,.22,0,1.6);
    else
        X = T(X,-.15,.28,.26,.24,0,.44);
    end
    if strcmp(Kolor,'rand')
        ind = ceil(rand(1,1)*Nc);
        plot(X(1),X(2),'Marker',Mrkr,'markersize',MS, ...
            'MarkerFaceColor',K(ind,:),'MarkerEdgeColor',K(ind,:))
    else
        plot(X(1),X(2),'Marker',Mrkr,'markersize',MS, ...
            'MarkerFaceColor',Kolor,'MarkerEdgeColor',Kolor)
    end
    pause(1/k)

end

%plot(X(:,1),X(:,2),'.','markersize',1)
% The transformation T

axis off

function U = T(X,a,b,c,d,e,f)
U = zeros(1,2);
U(1) = a * X(1) + b * X(2) + e;
U(2) = c * X(1) + d * X(2) + f;