Add planeWave function to calculate plane wave fields in the frequency domain...

Add planeWave function to calculate plane wave fields in the frequency domain at specified locations and frequencies

Sound Fields/MATLAB/Functions/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
\ No newline at end of file