function g = G_grad(x,y,k)
%
% function g = G_grad(x,y,k)
%
% This function evaluates the gradient of the 2D free-field Green's function, G(x,y), for the
% Helmholtz equation with respect to the components of x.
%
% On entry:
%
% x    is the length 2 position vector of the observation point
% y    is the position vector of the source (same size as x)
% k    is the wave number
%
% On exit:
%
% g    is the value of the gradient of the free-field Green's function
%
R = norm(x-y);
g = (x-y)*(i*k/4) * besselh(1,k*R)/R;