% bfs.m
% This draws the Schottky circles, based on the algorithm described
% in Box 12: Breadth-first search (p. 115).
%
% SYNTAX:   bfs.m(levstart,levend)
%
% levstart = level to start of plotting
% levend = level to end
% Here, level refers to the Schottky circle level.

function bfs(levstart,levend);

% 0. Initialize axes
axis equal;
axis ([-2.5 2.5 -2.5 2.5]);
hold on;
colors=['r' 'y' 'g' 'b' 'm' 'c']; % colors array

% 1. Initialize generators and circles ready.
theta=pi/4;
b=[1 cos(theta); cos(theta) 1]/sin(theta); % generators
a=[1 i*cos(theta); -i*cos(theta) 1]/sin(theta);
gens=cell(4,1);
gens{1}=a;gens{2}=b;gens{3}=inv(a);gens{4}=inv(b);
invs=[3 4 1 2]; % inverse of 1 is 3, etc.
group=cell(1000000,1); % to store words, eg. aaBA
circ(1).cen=sec(theta)*i;
circ(1).rad=tan(theta);
for j=2:4
   circ(j).cen=circ(j-1).cen*-i;
   circ(j).rad=tan(theta);
end

% 2. Enumeration loop
num=[1 5]; % indices where levels 1 and 2 begin
% lev=1
for i=1:4
   group{i}=gens{i}; % add word to list
   tag(i)=i; % tag with the last generator used
end
   
inew=5; % initialize counter
for lev=2:levend-1
   for iold=num(lev-1):num(lev)-1 % words from the previous level
      for j=1:4
         if j~=invs(tag(iold)) % avoid cancellations
            group{inew}=group{iold}*gens{j}; % add word to list
            tag(inew)=j; % tag with the last generator used
            inew=inew+1;
         end
      end
   end
   num(lev+1)=inew; % record end of level
end



% 3. Output loop
if levstart==1
   for j=1:4
      disc(circ(j),'r');
   end
end

for lev=max([1 levstart-1]):levend-1
   for i=num(lev):num(lev+1)-1
      clr=colors(mod(lev,6)+1); % get color from colors array
      for j=1:4
         if j~=invs(tag(i)) % avoid cancellations
            newcirc=mobius_on_circle(group{i}, circ(j));
            disc(newcirc,clr);
         end
      end
   end
   drawnow;
end

% End of main program segment. Supplementary functions follow.

% Box 10 on p.91
function D=mobius_on_circle(T,C);
if T(2,2)/T(2,1)+C.cen==0
   z=inf;
else
   z= C.cen - (C.rad*C.rad)/(T(2,2)/T(2,1)+C.cen)';
end
D.cen=mobius_on_point(T,z);
D.rad=abs(D.cen-mobius_on_point(T,C.cen+C.rad));

% Box 8 on p.75
function w=mobius_on_point(T,z);
if z==inf
   if T(2,1)==0
      w=inf;
   else
      w=T(1,1)/T(2,1);
   end
else
   if T(2,1)*z+T(2,2)==0
      w=inf;
   else
      w=(T(1,1)*z+T(1,2))/(T(2,1)*z+T(2,2));
   end
end

% disc - Draws a disc.
function disc(circ, colour, nsides)
if nargin == 2 % number of arguments is only 2
   nsides = 64;
end
nsides = round(nsides);  % make sure it is an integer
a = [0:pi/nsides:(2+1/nsides)*pi];
fill(circ.rad*cos(a)+real(circ.cen),circ.rad*sin(a)+imag(circ.cen),colour);