%%%%%%%%%%%%%%%%%%%%% Schwartz Quadrilateral Constructor %%%%%%%%%%%%%%%%%%%%%
%
% QSCHWATRZ
% Usage: QSCHWATRZ(s,t,u,v,h)
% Inputs:  s,t,u,v  - real numbers with 1/s +1/t + 1/u+ 1/v < 2 
%          h  - sidelength h must satisfy:
%          cosh(h) > (cos(a)*cos(b)+1)/(sin(a)*sin(b))
%          where a = pi/s, b = pi/t
%
% Output: matrix of numbers [ -L, L, z,w] such that  L is a positive real, 
%         the hyperbolic angles at the vertices are given by
%       
%               vertex         hyperbolic angle   
%
%                 -L                pi/s             
%                  L                pi/t
%                  z                pi/u            
%                  w                pi/v       
%         and the hyperbolic length from -L to L is h

function Q = qschwartz(s,t,u,v,h)             

%% set up constants, angles and tests   
I = sqrt(-1);
eps = 1.e-15;  % mfminsearch tolerance 
a = pi/s;   
b = pi/t;                   
c = pi/u;   
d = pi/v;

angletest = 2*pi-a-b-c-d;                                   % this should be positive
LawCosIItest = cosh(h) - (cos(a)*cos(b)+1)/(sin(a)*sin(b)); % this should be positive
if angletest < 0
   error('angles fail hyperbolic test')
end;
if LawCosIItest < 0
   error('sidelength not long enough')
end;
  
%% set up euclidean data of the lines 

% euclidean sidelength
L = tanh(h/4);  

% find the centers of the vertical sides
x0 = -L; y0=0; dydx = tan(a);
if s == 2 
   zL = [(1+x0^2)/(2*x0);0];
else
   A = [2*x0, 2*y0; 2, 2*dydx];
   B = [1+x0^2+y0^2; 2*x0+2*y0*dydx];
   zL = inv(A)*B;
end;
x0 = L; y0=0; dydx = -tan(b);
if t == 2 
   zR = [(1+x0^2)/(2*x0);0];
else
   A = [2*x0, 2*y0; 2, 2*dydx];
   B = [1+x0^2+y0^2; 2*x0+2*y0*dydx];
   zR = inv(A)*B;
end;

searchoptions = optimset('TolFun',eps,'MaxFunEval',1000);
[zT,success] = fminsearch('intpencil',[0,1],searchoptions,zL,zR,d,c);

zL = zL(1)+I*zL(2); 
zR = zR(1)+I*zR(2);
zT = zT(1)+I*zT(2);

z = intlines(zR,zT);
w = intlines(zL,zT);

Q = [-L,L,z,w];

boundary('','','w');
hold on

plotline(zL,'w');
plotline(zR,'w');
plotline(zT,'w');
plot([-1,1],[0,0],'w');