clear all;
close all;

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

x = -n:(n-1);
[X Y] = meshgrid(x);
R = sqrt(X.^2+Y.^2);
f = 1/(pi*(a/2)^2)*((R/a)<0.5);

F = fft2(f);
F = exp(1j*2*pi*(n+1)*((X-1)+(Y-1))/(2*n)).*F;
F = fftshift(F);

% Analytical transform
u = x/(2*n);
Fa_0deg = 1/pi*besselj(1,2*pi*(a/2)*u)./(a/2*u);

subplot(121);
f_0deg = F(n+1,:);
plot(u,f_0deg,'.');
xlim([min(u) max(u)]);
ylabel('0-deg slice through 2D FFT');
hold on;
plot(u,Fa_0deg,'r');
title(sprintf('a = %d, n = %d',a,n));
legend('FFT','Analytic','Location','NorthEast');


subplot(122);
f_45deg = diag(F);
plot(sqrt(2)*u,f_45deg,'.');
xlim(sqrt(2)*[min(u) max(u)]);
hold on;
ylabel('45-deg slice through 2D FFT');
Fa_45deg = 1/pi*besselj(1,2*pi*(a/2)*sqrt(2)*u)./(a/2*sqrt(2)*u);
plot(sqrt(2)*u,Fa_45deg,'r');