function p = bem_field_point(xf, cen, nrm, ar, k, p_surf, q_surf)
%BEM_FIELD_POINT  Evaluate the acoustic pressure at an exterior field point.
%
%   P = BEM_FIELD_POINT(XF, CEN, NRM, AR, K, P_SURF, Q_SURF)
%
%   Once the surface pressure P_SURF and normal derivative Q_SURF are known,
%   the pressure at any exterior point XF follows from the representation
%   integral (Helmholtz integral equation with the c-coefficient = 1 in the
%   exterior):
%
%       p(xf) = - integral_S [ G q  -  (dG/dn) p ] dS
%
%   evaluated here with the same constant-element quadrature used in the
%   solve (one Gauss point at each element centroid).
%
%   Inputs:
%       XF      - (1 x 3) exterior field-point coordinates
%       CEN     - (N x 3) element centroids
%       NRM     - (N x 3) outward unit normals
%       AR      - (N x 1) element areas
%       K       - wavenumber
%       P_SURF  - (N x 1) surface pressure (from bem_acoustic_solve)
%       Q_SURF  - (N x 1) surface dp/dn
%
%   Output:
%       P       - complex acoustic pressure at XF
%
%   For a SCATTERING problem this returns the SCATTERED pressure; add the
%   incident wave p_inc(xf) separately to obtain the total field.
%
%   See also BEM_ACOUSTIC_SOLVE, BEM_RADIATE_SPHERE, BEM_SCATTER_SPHERE.

    r  = cen - xf(:).';
    rn = sqrt(sum(r.^2, 2));

    g     = exp(-1j * k * rn) ./ (4 * pi * rn);
    rdotn = sum(r .* nrm, 2);
    dgdn  = g .* (rdotn ./ rn) .* (-1j * k - 1.0 ./ rn);

    p = -sum( (g .* q_surf(:) - dgdn .* p_surf(:)) .* ar );
end
