clc; clear; close all;

% consider the unit disk as the domain. R will be the radius of the
% Poincare model
R = 1.1;

%----- plot a bunch of geodesics for a fixed direction
s = -0.9:0.05:0;
nth1 = 6; dth1 = 2*pi/nth1;
thetas1 = dth1/2:dth1:2*pi-dth1/2;

figure; subplot(1,2,1);
plot(exp(1i*[0:2*pi/100:2*pi, 0]), 'b'); hold on;


for j = 1:length(s),
    for k = 1:nth1,
        [x,y] = UWgeodesic(s(j),thetas1(k),30,R);
        plot(x,y,'r'); hold on
    end
end
axis equal; axis([-1 1 -1 1]);
title(['geodesics, R =',num2str(R)]);

%----- backprojection geometry

x0 = [0.4, -.1]; y0 = [0.8, -.1];

nth1 = 20; dth1 = 2*pi/nth1;
thetas1 = dth1/2:dth1:2*pi-dth1/2;

subplot(1,2,2); plot(exp(1i*[0:2*pi/100:2*pi, 0]), 'b'); hold on;

for j = 1:length(x0),
    for k = 1:nth1,
        s0 = UWhyperS(x0(j),y0(j),thetas1(k),R);
        [x,y] = UWgeodesic(s0,thetas1(k),30,R);
        plot(x,y,'r'); hold on
    end
end
plot(x0,y0,'ko');
axis equal; axis([-1 1 -1 1]);
title('Backprojection geometry')
