function [x,y,r0,alpha0] = geodesic(s,theta,npts,R)
%GEODESIC: computes a discretized version of the unique curve 
% \gamma (s,\theta,t) in the Poincare disk of radius R >= 1 such that 
% \gamma (s,\theta, 0) = s (\cos\theta, \sin\theta)
% \dot\gamma (s,\theta,0) = (-\sin\theta, \cos\theta)
%
% OUTPUTS: 
%     - x,y: the coordinates of the points of that curve included in the
%     unit disk
%     - r0: when s is not 0, radius of the circle
%     - alpha0: when s is not 0, 1/2 of the aperture of the arc-circle 
%     describing the curve

if nargin == 3,
    R = 1;
end

if R<1,
    error('R should be greater than 1.')
end


if abs(s) > 1,
    error(['s = ',num2str(s),', too large !']);
end

if s == 0, % construct a straight line of direction (-\sin\theta, \cos\theta)
    
    x = -sin(theta)*(-1+1/npts:2/npts:1-1/npts);
    y = cos(theta)*(-1+1/npts:2/npts:1-1/npts);
    
elseif abs(s) < 1 % construct an arc-circle

    r0 = .5*(R*R/s - s); % radius of the circle
    rho0 = .5*(R*R/s + s); % distance of the center from zero
    
    % range of angles of the arc-circle
    alpha0 = abs(acos( (rho0^2 + r0^2 - 1)/(2*rho0*r0) ));
    
    alphamin = theta + pi - alpha0;
    alphamax = theta + pi + alpha0;
    
    dalpha = (alphamax-alphamin)/npts;
    alphas = alphamin+dalpha/2:dalpha:alphamax-dalpha/2;
    
    x = rho0*cos(theta) + r0*cos(alphas);
    y = rho0*sin(theta) + r0*sin(alphas);
else
    x=[];
    y=[];
    r0=0;
    alpha0=0;
end

end