Nx=3;
Ny=6;

xmax=10;
dx=0.1;
sigma=0.1;
a1=2.*xmax./Nx;
a2=2.*xmax./Ny;
bx=1
by=0.5;
cx=-1
cy=1;

[x,y]=meshgrid(-xmax:dx:xmax);
A=0.*x;
B=0.*x
C=0.*x
nx=ceil(xmax/a1);
ny=ceil(xmax/a2);

for a=-nx:nx;
    for b=-ny:ny
        xpos=a1.*a;
        ypos=a2.*b;
        A=A+exp(-((x-xpos).^2+(y-ypos).^2)./sigma.^2);
        B=B+exp(-((x-xpos-bx).^2+(y-ypos-by).^2)./sigma.^2);
         C=C+exp(-((x-xpos-cx).^2+(y-ypos-cy).^2)./sigma.^2);
    end
end
D=A+0.0.*B+0.0.*C;
kappa=0.0;
D=D.*exp(-kappa.*y);
subplot(1,2,1)
imagesc(D)
axis equal
xlim([1 size(A,1)])
ylim([1 size(A,2)])
subplot(1,2,2)
imagesc(abs(fftshift(fft2(D))).^2)
axis equal
xlim([80 120])
ylim([80 120])
