From 5c0313ed6d439d9c3d33f90f4abed5cb241cbffc Mon Sep 17 00:00:00 2001
From: ZaellixA <axilleaz@protonmail.com>
Date: Tue, 30 Jul 2024 15:42:03 +0100
Subject: [PATCH] Add planeWave function to calculate plane wave fields in the
frequency domain at specified locations and frequencies
---
Sound Fields/MATLAB/Functions/planeWave.m | 120 ++++++++++++++++++++++
1 file changed, 120 insertions(+)
create mode 100644 Sound Fields/MATLAB/Functions/planeWave.m
diff --git a/Sound Fields/MATLAB/Functions/planeWave.m b/Sound Fields/MATLAB/Functions/planeWave.m
new file mode 100644
index 0000000..2761008
--- /dev/null
+++ b/Sound Fields/MATLAB/Functions/planeWave.m
@@ -0,0 +1,120 @@
+%% 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
--
GitLab