function fern ()
clf; hold on;
drawFern (10, 5000);
axis equal;

function drawFern (n, m)
for i=1:n
    v = random ('Uniform', 0, 1, 2, 1);
    vs = computePoints (v, m);
    h = plot (vs(:,1), vs(:,2), 'r.');
    set (h, 'MarkerSize', 1);
end

function vs = computePoints (v, n)
t1 = [0.0 1.6]';
r1 = rotate (-2.5, -2.5, 0.85, 0.85);
t2 = [0.0 1.6]';   
r2 = rotate (49,   49,   0.3,  0.34);
t3 = [0.0 0.44]';
r3 = rotate (120,  -50,   0.3,  0.37);
t4 = [0.0 0.0]';
r4 = rotate (0,    0,   0.0,  0.16);

p = [81 9 9 1];
cp = cumsum (p);

x = random ('Uniform', 0, sum(p), 1, n);
vs = ones (n, 2);

for i=1:length(x)
    if     (x(i) < cp(1))
        v = r1 * v + t1;
    elseif (x(i) < cp(2))
        v = r2 * v + t2;
    elseif (x(i) < cp(3))
        v = r3 * v + t3;
    else
        v = r4 * v + t4;
    end
    vs(i,1) = v(1);
    vs(i,2) = v(2);
end


function res = rotate (phi, psi, r, s)
phi = pi * phi / 180;
psi = pi * psi / 180;
res = [r*cos(phi) -s*sin(psi); r*sin(phi) s*cos(psi)];