%%%%%%%%%%%%%%%%%%%%% Parametrization of Geodesic Segment %%%%%%%%%%%%%%%%%%%%
%
% GEOSEG
% Usage: GEOSEG(seg)
% Inputs:  seg = [x y] -  x,y are points inside or on the unit circle  
%
% Outputs: the hyperbolic geodesic arc connecting x and y is constructed
%          and is returned as a vector of complex points in the plane  

function w = geoseg(seg)             

I = sqrt(-1);
T  = 0:(1/(narcpnts-1)):1; 
oT = ones(size(T))-T;
x = seg(1); y = seg(2);     
   z = x*T+y*oT;                  % construct segment from x to y 
   m = (x+y)/2;                   % midpoint of 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       
   L = lft([x m y], [x mc y]);    % LFT mapping chord to arc
   w = lftact(L,z);               % map chord to arc 

   clear T;
   clear oT;
   clear z;