% spcg.m
% This draws the limit set of special parabolic commutator groups
% as described on Box 21 on p. 229.
%
% SYNTAX:    spcg(tra,trb,sgn,fine)
%
% tra=tr(a), trb=tr(b)
%  sgn = { 0 : take the positive square root of discriminant for tr(ab)
%        { 1 : take the negative square root of discriminant for tr(ab)
%  fine = { 0 : do a coarse plot
%         { 1 : do a fine plot
% =====================================================================
% Given tra, trb, the program computes trab to be one of the solutions
% to the quadratic (depending on choice of sgn):
%   x^2 - tra trb x + tra^2 + trb^2 = 0
% and computes
%   z0 = (trab - 2) trb / (trb trab - 2tra + 2i trab)
% and generator matrices a,b as follows:
%   a = [       tra/2      (tra trab - 2trb + 4i) / (2trab + 4) / z0 ]
%       [ (tra trab - 2trb - 4i)z0 / (2trab - 4)        tra/2        ]
%   b = 1/2 [ trb - 2i    trb  ]
%           [   trb    trb + 2i]
% The program proceeds to plot the limit set of <a,b>.
% For example, you may wish to try (2,2,1,0) or (1.91+.05i,1.91+.05i,0,0)

function spcg(tra,trb,sgn,fine)

global epsilon epsilonmin fp gens level lvlmax oldpoint oldppoint endpt numplots numb4plot plots;

switch fine
   case{0}
      disp('Coarse setting.');
      eps=.01;epsm=.001;lvlm=100;numb4=100; 
   case{1}
      disp('Fine setting.');
      eps=.002;epsm=.0001;lvlm=100;numb4=1000; 
end

set(0,'RecursionLimit',lvlm+100)
epsilon=eps; % distance must be less that epsilon in order to plot
epsilonmin=epsm; % distance must be greater than epsmin to plot
lvlmax=lvlm; % maximum level to consider
numb4plot=numb4; % number of points to compute before plotting
level=0;
plots=1;

% Initialize axes
axis equal;
hold on;

% 1. Initialize generators
trab=(tra*trb+(-1)^sgn*sqrt(tra^2*trb^2-4*(tra^2+trb^2)))/2;
z0=(trab-2)*trb/(trb*trab-2*tra+2i*trab);
a=[tra/2 (tra*trab-2*trb+4i)/(2*trab+4)/z0;
   (tra*trab-2*trb-4i)*z0/(2*trab-4) tra/2]
b=[trb-2i trb; trb trb+2i]/2 % generators a & b
A=inv(a); B=inv(b);
gens=cell(4,1);
gens{1}=a;gens{2}=b;gens{3}=A;gens{4}=B;
invs=[3 4 1 2]; % inverse of 1 is 3, etc.
repet=cell(4,3); % See p. 218-219
for k=1:4
   begpt(k)=fixpt(gens{mod4(k+1)}*gens{mod4(k+2)}*gens{mod4(k+3)}*gens{k});
   endpt(k)=fixpt(gens{mod4(k-1)}*gens{mod4(k-2)}*gens{mod4(k-3)}*gens{k});
   repet{k,2}=gens{k};
end

% Repetends to consider
repet{1,1}=b*A*B*a;repet{1,3}=B*A*b*a;
repet{2,1}=A*B*a*b;repet{2,3}=a*B*A*b;
repet{3,1}=B*a*b*A;repet{3,3}=b*a*B*A;
repet{4,1}=a*b*A*B;repet{4,3}=A*b*a*B;
for i=1:4
   for j=1:3
      fp(i,j)=fixpt(repet{i,j}); %Compute fixed points to be considered
   end
end
   
numplots=0;   
% 2. Explore the tree!
for k=1:4
   oldpoint=begpt(k);
   exploreTree(gens{k},k);
end

% End of main program segment. Supplementary functions follow.

% Box 17: Recursive depth-first search (p. 152) with adaptation (p. 183)
function exploreTree(X,l);
global epsilon epsilonmin fp gens level lvlmax oldpoint oldppoint endpt numplots numb4plot plots;
if numplots>=numb4plot*plots
   numplots
   plots=plots+1;
   drawnow; % updates the Matlab plot
            % (otherwise, Matlab only draws at the end.)
end

for k=l-1:l+1
   tag=mod4(k);
   Y=X*gens{tag};
   newpoint=mob(Y,endpt(tag));
   dist=abs(newpoint-oldpoint);
   for j=1:3
      z(j)=mob(Y,fp(tag,j));
   end
   termNow=0;
   if dist<epsilon & abs(z(1)-z(2))<=epsilon & abs(z(2)-z(3))<=epsilon 
      termNow=1; %let's terminate at the following 'if'-statement!
   end
   if termNow==1 | level>=lvlmax %then draw!
      if numplots==0
         oldppoint=newpoint;
      end
      if dist>=epsilonmin & level>1
         drawline(oldppoint,newpoint);
         numplots=numplots+1;
         oldppoint=newpoint;
      end
      oldpoint=newpoint;
   else %not yet: explore further!
      level=level+1;
      exploreTree(Y,tag); %recursive call
      level=level-1;
   end
end % End of function exploreTree(X,l);

function drawline(point1,point2)
x=linspace(real(point1),real(point2),2);
y=linspace(imag(point1),imag(point2),2);
line(x,y);
line(-x,-y); %applies 180 degree rotation automatically

function w=fixpt(T); % Note 3.3 on p.78, to determine fixed point
a=T(1,1); b=T(1,2); c=T(2,1); d=T(2,2);
if c==0
   w=inf; %modified in Ch8
else
   w=(a-d+sqrt( (a+d)^2-4 ) )/(2*c);
end



function w=mob(T,z); %Mobius on point (Box 8 on p.75)
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

function m=mod4(n); % = {1,2,3,4} mod 4
m=mod(n-1,4)+1;