clear all;
close all;

a = 32;
n = 2^nextpow2(4*a);

r = 0:(n-1);
dr = r(2)-r(1);
f = 1/(pi*(a/2)^2).*((r/a)<0.5);
plot(r,f);

q = linspace(0,(n-1)/(2*n),n);
for i = 1:n
    F(i) = 2*pi*sum(r.*f.*besselj(0,2*pi*r*q(i)))*dr;
end

plot(q,F,'.-');
xlabel('q');
xlim([min(q) max(q)]);
grid on;
title('Problem 3(a): Numerical Hankel transform, a = 32');