function [ X, Y, R ] = fractal_feto( N, Ps )
% [ X, Y, R ] = fractal_feto( N, Ps )
%
% exemplo de uso:
% [ X, Y, R ] = fractal_feto( 1e5 );
% figure, plot( X, Y, '.g', 'markersize', 1 )

if nargin<2
    Ps = [.85 .07 .07 .01];
end

X(N) = 0;
Y(N) = 0;
R(N) = 0;

for n = 1:N-1
    R(n+1) = aleatorio_com_probabilidade(Ps);
    switch( R(n+1) )
        case{ 1 }
            X(n+1) = 0.85*X(n) + 0.04*Y(n);
            Y(n+1) = -0.004*X(n) + 0.85*Y(n)+1.6;
        case{ 2 }
            X(n+1) = -0.15*X(n) + 0.28*Y(n);
            Y(n+1) = 0.26*X(n) + 0.24*Y(n) + 0.44;
        case{ 3 }
            X(n+1) = 0.2*X(n) - 0.26*Y(n);
            Y(n+1) = 0.23*X(n) + 0.22*Y(n) + 1.6;
        otherwise
            X(n+1) = 0;
            Y(n+1) = 0.16*Y(n);
    end    
end

if nargout==0
    plot( X, Y, '.g', 'markersize', 1 )
    clear X
end
