function [x,y,r0,alpha0] = geodesic(s,theta,npts)

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

if s == 0,
    
    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)<=.9
    r0 = .5*(1./s - s);
    y0 = .5*(1./s + s); % distance of the center from zero
    
    % angle range
    alpha0 = asin( 1/y0 );
    
    alphamin = theta + pi - alpha0;
    alphamax = theta + pi + alpha0;
    
    dalpha = (alphamax-alphamin)/npts;
    alphas = alphamin+dalpha/2:dalpha:alphamax-dalpha/2;
    
    x = y0*cos(theta) + r0*cos(alphas);
    y = y0*sin(theta) + r0*sin(alphas);
else
    x=[];
    y=[];
    r0=0;
    alpha0=0;
end

end