function s = hyperS(x,y,theta,R)
%HYPERS s(x,y,theta)
% for a given point (x,y) in the unit disk and \theta\in\S^1, returns the
% unique geodesic of second coordinate \theta that passes through (x,y)
% can handle several thetas

if nargin == 3, 
    R = 1; 
end

if x*x+y*y > 1,
    error(['(x,y) = (',num2str(x),',',num2str(y),') is outside the unit disk.']);
end

p = x.*cos(theta) + y.*sin(theta);

idx = find(p ~=0 );

s = zeros(size(theta)); % s=0 wherever p=0

beta = (R*R + x.*x+y.*y)./2./p(idx);
s(idx) = beta - sign(beta).*sqrt(beta.*beta - R*R);

%-- version handling only one theta
 
% if p == 0
%     s = 0;
% else
%     beta = (R^2+x.*x+y.*y)./2./p;
%     s = beta - sign(beta).*sqrt(beta.*beta - R^2);
% end

end
        