clear all; close all;

[X,Y] = meshgrid(-512:511);
J = zeros(size(X));

num_cycles = 11;
period = 1024/num_cycles;
K = 2*pi/period;

Ntheta = 200;
theta = linspace(0,pi,Ntheta);
for i = 1:Ntheta
    kx = K*cos(theta(i));
    ky = K*sin(theta(i));
    J = J + cos(kx*X+ky*Y);
end
J = J/Ntheta;
subplot(121);
surf(X,Y,J);
view([0 90]);
shading interp;
xlim([-512 511]);
ylim([-512 511]);
axis equal;
title(sprintf('numCycles=%d, numCosines=%d',num_cycles,Ntheta));

subplot(122);
plot(-512:511,J(513,:),'.');
xlim([-512 511]);
ycont = besselj(0,K*(0:511));
hold on;
plot(0:511,ycont,'r','Linewidth',2);
legend('Numerically integrated','Analytic','Location','NorthWest');