function s = hyperS(x,y,theta)
%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)

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);

if p == 0
    s = 0;
else
    beta = (1+x*x+y*y)/2/p;
    s = beta - sign(beta)*sqrt(beta*beta-1);
end

end
        