% hyper_tess.m
% Plot hyperbolic triangulation using right triangles with 
% opposite angles   pi/p and pi/q.
% colors alternate between white and a color set by hand in 'geo_patch'
% Christopher Bishop  Oct 3, 2020
function tess5(p,q)
   close all;
   % draw the limit circle
   circle = exp(2*pi*i*[0:.01:1]);
   figure; hold on; axis equal; axis off;
   plot(circle,'k','LineWidth',3);
   % define vertices of first triangle
   a=0;
   b=tanh(acosh(cos(pi/p)/sin(pi/q))/2);
   c=i*tanh(acosh(cos(pi/q)/sin(pi/p))/2);
   % set number of Eulclidean segments used to approx hyper segments
   n=30;
   % define list of triangles; three complex vertices and a parity 0,1
   triangles=[a,b,c,0];
   new_triangles=[];
   % continue while number of triangles is small, set to 100 for debugging
   % set to around 10,000 or larger to fill whole disk
   % %%%%%%%%%%%%%%%%      set number of triangles %%%%%%%%%%%%%
   while size(triangles,1) < 10000
      counter=0;
      % for each triangle, define 3 new ones by reflecting over sides
      for k=1:size(triangles,1)
         a=triangles(k,1);
         b=triangles(k,2);
         c=triangles(k,3);
         parity=triangles(k,4);
         new_triangles(counter+1,1:4)=[b,a,tri_flip(a,b,c),1-parity];
         new_triangles(counter+2,1:4)=[c,b,tri_flip(b,c,a),1-parity];
         new_triangles(counter+3,1:4)=[a,c,tri_flip(c,a,b),1-parity];
        counter=counter+3;
      end % for k 
      % add new triangles to old ones
      % uniquetol defines uniqueness for floating points up to a tolerence
      triangles = [triangles;new_triangles];
      % remove repeated triangles
      [C,IA,IC]=uniquetol([real(triangles),imag(triangles)],'ByRows',true);
      triangles=triangles(IA,:);
      % now repeat if number of triangles is not too large
   end % while 
   % plot the found triangles. To get checkerboard only plot one parity
   for k=1:size(triangles,1)
      a=triangles(k,1);
      b=triangles(k,2);
      c=triangles(k,3);
      parity=triangles(k,4);
      if parity==1
         geo_patch(a,b,c,n);
      end % if
   end % for k
   %set(gcf, 'Renderer', 'Painters');
   %print('test_save_dpi1000.eps','-r1000')
return 
end  % function tess

%    -------------    sub functions ---------------

% reflect a hyperbolic triangle over one side
function out=tri_flip(a,b,c)
f=@(z)(z-a)./(1-conj(a).*z);
finv=@(w)  (w+a)./(1+conj(a).*w);
out=finv(f(b)*conj(f(c)/f(b)));
end

% find points on geodesic segment
function out=geo(a,b,n)
f=@(z)(z-a)./(1-conj(a).*z);
finv=@(w)  (w+a)./(1+conj(a).*w);
out=finv([0:1/n:1]*f(b));
end

% plot a single geodesic segment
function geo_plot(a,b,n)
z=geo(a,b,n);
x=real(z); y=imag(z);
plot(x,y,'k');
end

% plot a geodesic triangle, fill interior black
function out=geo_patch(a,b,c,n)
edge1=geo(a,b,n);
edge2=geo(b,c,n);
edge3=geo(c,a,n);
z=[edge1,edge2,edge3];
x=real(z);
y=imag(z);
% set gray to desired value for checkerboard
gray=[ 0 0 1];  % bright blue
%gray=[.5 .5 1];  % use for light blue triangles 
gray=[1 0 0 ];  % bright red
%gray=[0 0 0];  use for black triangles
patch(x,y,gray);
end

% plot a geodesic triangle, leave interior empty
function out=geo_outline(a,b,c,n)
edge1=geo(a,b,n);
edge2=geo(b,c,n);
edge3=geo(c,a,n);
z=[edge1,edge2,edge3];
x=real(z);
y=imag(z);
plot(x,y,'k');
end

% returns poly approximation to geodesic triangle
function z=geo_bdy(a,b,c,n)
edge1=geo(a,b,n);
edge2=geo(b,c,n);
edge3=geo(c,a,n);
z=[edge1,edge2,edge3];
w=mean(z);
z=z+.001*(w-z);
end



