windfern(n, REPS)

function windfern(n, REPS)
%WINDFERN   MATLAB implementation of the Fractal Fern.
%   Michael Barnsley, "Fractals Everywhere", Academic Press, 1993.
%   WINDFERN(N) plots ten repetitions, N points per frame.
%   WINDFERN(N,REPS) plots REPS repetitions, N points per frame.
%   WINDFERN plots ten repetitions, 50000 points per frame.
%   Wind added by McKeeman

if nargin < 1, n = 50000; end
if n <= 2000, dotsize = 6; else dotsize = 1; end

% flatten curvature matrices to make available to JIT
p1  =  .85;  p2  =  .92;  p3  = .99;

A11 =  .85;  A12 =  .04;              A22 = A11;
B11 =  .20;  B12 = -.26;  B21 = .23;  B22 = .22; 
C11 = -.15;  C12 =  .28;  C21 = .26;  C22 = .24;
D22 =  .16;

b1 = 1.4;  b2 = 1.6; b3 = 0.44;

wave    = A12;
newleft = wave;
xs = zeros(1,n);
ys = zeros(1,n);

% set up plot
figure(gcf)                                       % on top
set(gcf, 'color','white','menubar','none', 'numbertitle','off',...
         'name', 'Fern in the Wind', 'doublebuff', 'on');
darkgreen = [0 2/3 0];                            % fern-ish
h = plot(xs, ys, '.','markersize',dotsize,'color',darkgreen);
axis([-4 4 0 10])                                 % set size
axis off                                          % make pretty
    
% start the demo
if nargin < 2; REPS = 10; end
for rp = 1:REPS
  oldleft = newleft;
  newleft = -wave+(rand-.5)/20;
  right   = +wave-(rand-.5)/20;
  delta   =  (right-newleft)/20;
  t1 = 0; t2 = 0;                                 % costs
  for w = [oldleft:delta:right, right:-delta:newleft]
    tstart = tic;                                 % start time compute
    x  = .5;  y = .5;
    for j = 1:n
      xs(j) = x;
      ys(j) = y;
      r = rand;
      if r < p1 
        t = A11*x + w*y;                          % rotate and scale
        y = -w*x + A22*y + b1;
        x = t;
      elseif r < p2
        b1 = 1.5;
        t = B11*x + B12*y;
        y = B21*x + B22*y + b2;
        x = t;
      elseif r < p3
        t = C11*x + C12*y;
        y = C21*x + C22*y + b3;
        x = t;
      else                                        % the stem
        y = D22*y;
        x = 0;
      end
    end
    t1 = t1 + toc(tstart);                        % end time compute
    tstart = tic();                               % start time plot
    set(h,'xdata', xs, 'ydata', ys);
    drawnow;
    pause(rand/100);                              % slow to human speed
    t2 = t2 + toc(tstart);                        % end time plot
  end
  %fprintf('rep=%d compute=%.1g secs plot=%.1g secs\n', reps, t1, t2);
end