% fractal "lite"
% a basic IFS fractal drawing program
% using the random method


ifs = input('Enter IFS data matrix -> ');

npts = input('Enter number of points to compute -> ');

[r,c] = size(ifs);
if c~=6
	error('Matrix must have 6 columns');
end;

% find a point on the fractal ...
pt = [0;0];
A = [ ifs(1,1:2); ifs(1,3:4) ];
b = [ ifs(1,5:6)]';

% ... by iterating the first function 20 times
for k=1:20
	pt = A*pt+b;
end;

% next we compute the weight of each function

detlist = zeros(1,r);
for k=1:r
	detlist(k) = abs(ifs(k,1)*ifs(k,4) - ifs(k,2)*ifs(k,3));
end;

% if all determinants are zero, set detlist to all 1's
if all(detlist < 3*eps)
	detlist = ones(1,r);
end;

% find all the zeros in detlist...
zlist = find(detlist < 3*eps);
least = min(detlist(find(detlist>3*eps)));

% ... and set them to a small value
detlist(zlist) = (least/10)* ones(size(zlist));
 
% convert detlist into a probability vector
detlist = detlist/sum(detlist);

% and compute a selection vector 
% by computing cumulative sums
selector = cumsum(detlist);


% we precompute our selections

rlist = rand(1,npts);
choice = zeros(1,npts);

for k=1:npts
	choice(k) = min(find(rlist(k) < selector));
end;

pts = zeros(2,npts); %  a place to hold the points

for k=1:npts
	p = choice(k);
	A = [ifs(p,1:2); ifs(p,3:4)];
	b = [ifs(p,5:6)]';
	pt = A*pt+b;
	pts(:,k) = pt;
end;
	
	
% plot the points
plot(pts(1,:),pts(2,:),'.')
axis([0 1 0 1]);
axis('square');