function [f_2D] = zoneplate2D(half_width, freq, sampling)

% using linear chirp in 1D case

if (nargin == 0)
   half_width = 51
   freq = 40000
   sampling = 0.0005 
elseif (nargin == 1)
   freq = 40000
   sampling = 0.0005 
else
   sampling = 0.0005 
end

t = -(sampling*half_width):sampling:(sampling*half_width);
t_len = length(t);

X = repmat(t,t_len,1);
Y = X';

f_2D = cos( (X.^2 + Y.^2)*freq );

imshow(f_2D,[-1 1])