clear all; close all;
source = 'hw7prob4data.txt';
data = load(source);
C = length(data(:,1));

x = data(:,1);
dx = x(2)-x(1);
f = data(:,2);
subplot(311);
plot(x,f);
xlabel('x');
ylabel('projection data');
title(source);
xlim([x(1) x(end)]);
ylim([-1 1]);

% Take the fourier transform of the data
NFFT = 2^(nextpow2(length(x))+2);
Fs = 1/dx;
u = Fs/2*linspace(0,1,NFFT/2+1);
du = u(2)-u(1);
X = fft(f,NFFT)*dx;
X = real(X(1:NFFT/2+1));
subplot(312);
plot(u,X,'.-');
xlim([0 5]);
xlabel('u');
ylabel('Fourier transform of projection');

% Take inverse Hankel transform numerically
r = linspace(min(x),2*max(x),1e3);

g = zeros(size(r));
for i = 1:length(r)
    g(i) = 2*pi*sum(u.*X'.*besselj(0,2*pi*u*r(i)))*du;
end
subplot(313);
plot(r,g);
xlabel('x');
ylabel('Reconstructed slice');
ylim([-1 1]);

% Make into image
[X Y] = meshgrid(linspace(-5,5));
R = sqrt(X.^2+Y.^2);
F_orig = zeros(size(X));
for i = 1:100
    for j = 1:100
        F_orig(i,j) = interp1(r,g,R(i,j),'nearest','extrap');
    end
end
figure;
imagesc(F_orig,[-1 1]);
axis equal;
colormap gray;