%%%%%%%%%%%%%%%%%%%%%%%%%%%% Half Plane  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% HPMAT
% Usage: HPMAT(tri)
% Inputs:  tri - a triangle [P Q R]  
%   
% Output:  a matrix L used to determine if an arbitrary point z is on the 
%          same side of the geodesic [P Q] as R 
%          the point z will be on the same side if lftact(L,z) has positive 
%          imaginary part

function L = hpmat(tri)             

   I = sqrt(-1);
   x = tri(1);
   y = tri(2);
   z = tri(3);

   m = (x+y)/2;                   % midpoint of Euclidean segment
   d = I*(x-y)/2;                 % vector orthogonal to vector to segment
   u = abs(m);                    % calculate some scalars 
   v = abs(d);                     
   rt = 2*real(m*conj(d))/(u^2-v^2-1);    
   s = rt/(1+sqrt(1+rt^2));              
   mc = m+s*d;                    % (euclidean) midpoint of arc       
   L1 = lft([x mc y], [-1 0 1]);   % LFT mapping H-segment to E-segment
   L3 = [I 0; 0 -I];
   w  = imag(lftact(L,z));
   if (w < 0) L = L3*L; end;