clc; clear; close all;

% consider the unit disk as the domain, and assume the image is supported
% in the disk of radius 0.9

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

figure(1);
plot(exp(i*[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);
        plot(x,y,'r'); hold on
    end
end
axis equal; axis([-1 1 -1 1]);
title('geodesics')

%----- backprojection geometry

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

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

figure(2);
plot(exp(i*[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));
        [x,y] = geodesic(s0,thetas1(k),30);
        plot(x,y,'r'); hold on
    end
end
plot(x0,y0,'ko');
axis equal; axis([-1 1 -1 1]);
title('Backprojection geometry')
