function fern(Nstep)
	%% generates fern like structure

	if nargin < 1, Nstep=15000; end

	%% transformation 	
	% x = a*x + b*y + e;
	% y = c*x + d*y + f
	% transformation matrix
	%	1	2	3	4	5	6	7
	%	a	b	c	d	e	f	p
	% column p is a probability of choosing given row for transformation	
	%transformation_matrix = [ ...
		%[0 	0 	0.00 	0.27 	0.5 	0 	0.02];...
		%[-0.139	0.263	0.246 	0.224	0.57 	-0.036	0.15]; ...
		%[0.17 	-.215 	0.222	0.176 	0.408 	0.0893	.13]; ...
		%[0.781	0.034 	-0.032 	0.739	0.1075	0.27 	0.70]; ...
		%];
	transformation_matrix = [ ...
		[0       0       0       0.16   0   0      0.01];...
		[0.85    0.04    -0.04   0.85   0   1.6    0.85];      ...
		[0.2     -0.26   0.23    0.22   0   1.6    0.07];      ...
		[-0.15   0.28    0.26    0.24   0   0.44   0.07];      ...
		];
	p=[0,0]; % initial point



	points=zeros(Nstep,2); % (x,y)  coordinates array

	%% make a choice 
	function choice = probabilistic_choice(probabilities)
		%% given array of probabilities choose one scenario
		Np = length(probabilities);
		acum_probabilities=zeros(1,Np);

		acum_probabilities(1) = probabilities(1);
		for k=2:Np
			acum_probabilities(k) = acum_probabilities(k-1) + probabilities(k);
		end

		r=rand(1,1);

		%% which choice
		for k=1:Np
			if ( r <= acum_probabilities(k) )
				choice = k;
				return;
			end
		end

		%% we should never be here if all probabilities add up to 1
		% but let's call it "omitted" probability
		choice=k+1;
	end



	function pnew = new_point(pold, transform)
		a=transform(1);
		b=transform(2);
		c=transform(3);
		d=transform(4);
		e=transform(5);
		f=transform(6);

		x=pold(1);
		y=pold(2);

		x = a*x + b*y + e;
		y = c*x + d*y + f;

		pnew=pold; % to create array with similar geometry to pold
		pnew(1)=x;
		pnew(2)=y;
	end


	for i = 1:Nstep
		choice=probabilistic_choice( transformation_matrix(:,7) );
		transform=transformation_matrix(choice,:);

		p=new_point(p, transform);
		points(i,:)= p;
	end

	figure();
	plot( points(:,1), points(:,2), 'g.');
end


		
