function tiling_poincare_disc() 
 
 

 clear all ,close all, 
hold on
 
center=[0,0];
radius=1;  
THETA=linspace(0,2*pi,2000);
RHO=ones(1,2000)*radius;
[X,Y] = pol2cart(THETA,RHO);
X=X+center(1);
Y=Y+center(2); 
fill(X,Y,'r'),   axis square , grid on
 
    axis([-1 1 -1 1]) 
    
    a=0.25;
    b=0.398;
     u1=0;u2=b;
     v1=-b*a*(sqrt(10+2*sqrt(5))); v2=b*a*(sqrt(5)-1);
     w1=-b*a*(sqrt(10-2*sqrt(5))); w2=-b*a*(sqrt(5)+1);
     y1=b*a*(sqrt(10-2*sqrt(5))); y2=-b*a*(sqrt(5)+1);
     x1=b*a*(sqrt(10+2*sqrt(5))); x2=b*a*(sqrt(5)-1);
   
     
    z1=u1+i*u2;  Z1=z1;
    z2=v1+i*v2;  Z2=z2; 
    z3=w1+i*w2;  Z3=z3;
    z4=y1+i*y2;  Z4=z4;
    z5=x1+i*x2;  Z5=z5;
    
    zz1=LINE(z1,z2);  zz2=LINE(z2,z3);  zz3=LINE(z3,z4);  zz4=LINE(z4,z5);  zz5=LINE(z5,z1); plot([zz1,zz2,zz3,zz4,zz5])
    
   
level=4; %level of recursion

tree(z1,z2,z3,z4,z5,level);
tree(z5,z1,z2,z3,z4,level);
tree(z4,z5,z1,z2,z3,level);
tree(z3,z4,z5,z1,z2,level);
tree(z2,z3,z4,z5,z1,level);
axis equal tight 
return

function tree(z1,z2,z3,z4,z5,n)
hold on
if n>0
          
    [z1,z2,z3,z4,z5]=pent(z1,z2,z3,z4,z5);        
    zz1=LINE(z1,z2);  zz2=LINE(z2,z3);  zz3=LINE(z3,z4);  zz4=LINE(z4,z5);  zz5=LINE(z5,z1); plot([zz1,zz2,zz3,zz4,zz5])
    
   
               
      tree(z1,z2,z3,z4,z5,n-1);
tree(z5,z1,z2,z3,z4,n-1);
tree(z4,z5,z1,z2,z3,n-1);
tree(z3,z4,z5,z1,z2,n-1);
tree(z2,z3,z4,z5,z1,n-1);
end

return
 

function [z] = LINE(a,b)

lambda=((1-a*conj(b))/(abs(1-a*conj(b))))*sqrt((1-a*conj(a))/(1-b*conj(b)));
    
z=-1:0.01:1;
    
z=(z.*(a-lambda*b) - (a+lambda*b))./(z.*(1-lambda)-(1+lambda));
    
function [z1,z2,z3,z4,z5] = pent(c,d,e,a,b)

z1=(b*conj(a)-a*conj(b)+(a-b+b*conj(b)*a-a*b*conj(a))*conj(c))/...
       (conj(a)-conj(b)+b*conj(b)*conj(a)-a*conj(a)*conj(b)+(a*conj(b)-conj(a)*b)*conj(c));

z2=(b*conj(a)-a*conj(b)+(a-b+b*conj(b)*a-a*b*conj(a))*conj(d))/...
       (conj(a)-conj(b)+b*conj(b)*conj(a)-a*conj(a)*conj(b)+(a*conj(b)-conj(a)*b)*conj(d));
  
z3=(b*conj(a)-a*conj(b)+(a-b+b*conj(b)*a-a*b*conj(a))*conj(e))/...
       (conj(a)-conj(b)+b*conj(b)*conj(a)-a*conj(a)*conj(b)+(a*conj(b)-conj(a)*b)*conj(e));
  
       
       z4 = a;
       z5 = b;
       
           
    
 