function G = fbt(r,g)
% Function to compute Hankel transform using backprojection
% method.  This version for nth-order transforms.
% Correct up to a sign.  For n odd, multiply result by
% (-1)^(n-1)/2 .
% To use:  G = fbt(r,g);
%
% by D. Holmgren, Mar 3 95 holmgren@phobos.astro.uwo.ca
rg = fftshift(r.*g);
F = fft(rg);
F = real(F);
nl = length(r);
G = zeros(nl,1);
xnu = linspace(0.,1.,nl)';
n = input('Order n: ');
nth = input('Number of theta points? ');
theta = linspace(0,pi/2,nth)';
w = ones(nth,1); w(1) = 0.5; w(nth) = 0.5;
for k = 1:nl
 for l = 1:nth
  rct = xnu(k)*cos(theta(l));
  Fi = table1([xnu F],rct);
  G(k) = G(k) + w(l)*Fi*cos(n*theta(l));
 end
end
G = (pi/nth) .* G;