planeWave.m

%% Calculate wavefield generated by plane waves in the frequency domain
% --------------------------------------------------
% Author: Achilles Kappis
% e-mail: axilleaz@protonmail.com
%
% Date: 30/07/2024 (DD/MM/YYYY)
%
% Copyright: MIT
% --------------------------------------------------
% Functionality: Calculate the wavefield generated by plane waves at
%                specified directions in the frequency domain.
% --------------------------------------------------
% Input
% 
% sDir [numeric]: The direction of incidence of the plane waves. This can
%                 be either an Nx2 matrix with N being the number of angles
%                 of incidence in [azimuth, elevation] format, expressed in
%                 radians, or an Nx3 matrix with the entried corresponding
%                 to the wavenumbers along each of the Cartesian axes. See
%                 "notes" for more information.
% 
% rPos [numeric]: The position of the receivers in Cartesian coordinates.
%                 The variable must be a matrix with dimensions Mx3 where
%                 M is the number of receivers and the columns represent
%                 the x, y and z coordinates respectively.
% 
% f [numeric]: The frequencies for which to calculate the wave field. It
%              must be a scalar or vector.
% 
% Q [numeric] (Optional): The source strength of the plane waves. This
%                         variable must be either a scalar corresponding to
%                         all sources having the same strength, or a matrix
%                         with dimensions NxK, where K is the number of
%                         frequencies of interest. The strengths can be
%                         complex numbers. [Default: 1].
% 
% c [numeric] (Optional): The speed of sound in m/s. This variable is not
%                         used if the wavenumbers are provided in the sDir
%                         argument. [Default: 343].
% 
% --------------------------------------------------
% Output
% 
% waveField [numeric]: The complex pressures at the position of the
%                      receivers. The variable is a matrix of dimensions
%                      MxNxF where M is the number of receivers, N the
%                      number of sources and F the number of frequencies.
% 
% --------------------------------------------------
% Notes
% 
% If the variable "sDir" provides the wavenumbers along the Cartesian axes,
% it is the responsibility of the user to make sure that the condition
% k^2 = kx^2 + ky^2 + kz^2 holds true.
% 
% --------------------------------------------------
function waveField = planeWave(sDir, rPos, f, Q, c)
    % ====================================================
    % Check for number of arguments
    % ====================================================
    narginchk(3, 5);
    nargoutchk(0, 1);

    % ====================================================
    % Validate input arguments
    % ====================================================
    % Validate mandatory arguments
    validateattributes(sDir, "numeric", {'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Source positions", 1);

    % Check form of sDir (directions or wavenumbers are provided)
    if sum(size(sDir, 2) ~= 2) > 0
        error("The size of the directions matrix must be either Nx2 (incidence angles), or Nx3 (wavenumbers along the Cartesian axes)");
    end

    validateattributes(rPos, "numeric", {'2d', 'ncols', 3, 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Receiver positions", 2);
    validateattributes(f, "numeric", {'vector', 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Frequencies", 3);

    % Validate optional arguments
    if nargin > 3 && ~isempty(Q)
        validateattributes(Q, "numeric", {'2d', 'finite', 'nonempty', 'nonnan'}, mfilename, "Source strength(s)", 5);

        % Make sure Q has the correct dimensions
        if ~isscalar(Q) && sum((size(Q) ~= [size(sDir, 1), length(f)])) > 0
            error("Source strength(s) must be either a scalar or its length must match the number of sources.");
        end
    else
        Q = 1;
    end

    if nargin > 4 && ~isempty(c)
        validateattributes(c, "numeric", {'scalar', 'finite', 'positive', 'nonnan', 'nonempty', 'real'}, mfilename, "Speed of propagation", 6);
    else
        c = 343;
    end

    % ====================================================
    % Pre-process
    % ====================================================
    % Make sure frequencies is a column vector
    f = f(:);

    % Calculate wavenumber from angles of indidence to wavenumber (if needed)
    if size(sDir, 2) == 2
        [kx, ky, kz] = sph2cart(sDir(:, 1), sDir(:, 2), ones(size(sDir, 1), 1));
        k = [kx, ky, kz] .* (2 * pi * permute(f, [3, 2, 1])/c);
    end

    % Create a source strength vector (if needed)
    if isscalar(Q)
        Q = ones(size(k, 1), length(f)) * Q;
    end

    % =============================================
    % Calculate wavefield
    % ============================================
    % Go through the frequencies
    for fIdx = length(f):-1:1
        waveField(:, :, fIdx) = Q(:, fIdx).' .* exp(-1j * rPos * k(:, :, fIdx));
    end
end