clear
close all

dz = 10e-2;
w0 = 2*pi*5e3;
c = 343.1;
x = linspace(-0.5,0.5,200);
dx = mean(diff(x));
y = linspace(-0.5,0.5,200);
dy = mean(diff(y));
[X,Y] = meshgrid(x,y);

r0 = sqrt(X.^2 + Y.^2 + dz.^2);
G0 = exp(-1i*w0/c*r0)./r0;
figure;
pcolor(real(G0));
shading interp
axis equal tight;
%%
r_ls = 0.15;
v0 = (sqrt((X+0.2).^2 + Y.^2)<r_ls) + (sqrt((X-0.2).^2 + Y.^2)<r_ls);
ixs = find(v0);

r0 = sqrt(X.^2 + Y.^2 + dz.^2);
figure;
pcolor(real(v0));
shading interp
axis equal tight;
%%
P = zeros(size(X));
X0 = X(:); Y0 = Y(:);
v0 = v0(:);
for n = 1 : numel(ixs)
    n
   r =  sqrt((X-X0(ixs(n))).^2 + (Y-Y0(ixs(n))).^2 + dz.^2);
   P = P + v0(ixs(n))*exp(-1i*w0/c*r)./r;
end
%%
figure;
pcolor(20*log10(abs(P)));
shading interp
axis equal tight;
%caxis([0,2e3])
%%
win = hann(size(P,1))*hann(size(P,2))';
Gk = fft2(win.*G0);
Pk = fft2(win.*P);
kx = ((0:size(P,2)-1)-size(P,2)/2)/size(P,2)*2*pi/dx;
ky = ((0:size(P,1)-1)-size(P,1)/2)/size(P,1)*2*pi/dy;
[Kx,Ky] = meshgrid(kx,ky);
K = w0/c;

iGk = 1./Gk;
mask = zeros(size(Gk));
mask(Kx.^2+Ky.^2<K.^2) = 1;
mask = fftshift(mask);
iGk = 1./Gk;

Vapp = Pk.*iGk.*mask;
Vapp(isnan(Vapp)) = 0;

figure;
subplot(1,2,1)
pcolor(reshape(v0,size(X,1),size(X,2)));
shading interp
axis equal tight;
subplot(1,2,2)
pcolor(ifftshift(abs(ifft2(Vapp))));
shading interp
axis equal tight;