function [ps] = planewave_cyl2d(x, a, k, m)
% PLANEWAVE_CYL2D Reflection of a plane wave from a cylinder in 2D
%   PS = PLANEWAVE_CYL2D(X, A, K, M) Computes the pressure field in points
%       X (N x 2 matrix) scattered by a cylinder of radius A located at the 
%       origin of the coordinate system, the incident field is a plane wave
%       traveling in the x direction with a wave number K. The parameter M
%       defines the order of the approximation (default: M = 100).
%
% References
%   [1] P. M. Morse, K. U. Ingard. Theoretical Acoustics. McGraw-Hill
%   (1968) pages: 401-403

% Copyright 2011-2021 Bence Oltean, Peter Fiala, Peter Rucz

% Assign default order
if nargin < 4
    m = 100;
end

% Convert coordinates to polar
[phi, r] = cart2pol(x(:,1), x(:,2));

g = [
    atan(-besselj(1,k*a)/bessely(1,k*a)), ...
    atan((besselj((1:m)-1,k*a)-besselj((1:m)+1,k*a)) ./ ...
         (bessely((1:m)+1,k*a)-bessely((1:m)-1,k*a)))
    ];

A = -1i*[1, 2*1i.^(1:m)].*exp(-1i*g).*sin(g);

ps = bsxfun(@times, cos(phi*(0:m)), A) .* ...
        besselh(repmat(0:m, size(r,1), 1), repmat(k*r, 1, m+1));

ps = conj(sum(ps,2));

end % of function planewave_cyl2d
