diff --git a/Sound Fields/MATLAB/Functions/rectRoomPtSrc.m b/Sound Fields/MATLAB/Functions/rectRoomPtSrc.m
new file mode 100644
index 0000000000000000000000000000000000000000..0247ce92fc58aee3cc0c63b201bb60cf200bb2e3
--- /dev/null
+++ b/Sound Fields/MATLAB/Functions/rectRoomPtSrc.m
@@ -0,0 +1,338 @@
+%% Calculate wavefield generated by point sourcesinside a rectangular room
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 29/05/2025 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate the wavefield generated by point sources at
+% specified locations inside a rectangular room of specified
+% dimensions and absorption.
+% --------------------------------------------------
+% Input
+%
+% sPos [numeric]: The position of the sources in Cartesian coordinates. The
+% variable should be a matrix with dimensions 3xN where N
+% is the number of sources and the rows represent the x, y
+% and z coordinates respectively.
+% IMPORTANT: The coordinates of the sources must NOT lie
+% outside the room (set in "roomSize" argument).
+%
+% rPos [numeric]: The position of the receivers in Cartesian coordinates.
+% The variable should be a matrix with dimensions 3xM where
+% M is the number of receivers and the rows represent the
+% x, y and z coordinates respectively.
+% IMPORTANT: The coordinates of the receivers must NOT lie
+% outside the room (set in "roomSize" argument).
+%
+% roomDims [numeric]: The dimensions of the room. This can either be a
+% scalar, in which case the room will correspond to a
+% cubic enclosure, or a vector with three elements
+% corresponding to the x, y and z dimensions of the
+% room. The values must be positive scalars.
+%
+% f [numeric]: The frequencies for which to calculate the wave field. It
+% must be a positive scalar or vector of positive values.
+%
+% nModes [numeric] (Optional): The (approximate) number of modes to be used
+% for the calculation of the sound field. This
+% value can either be a positive scalar
+% denoting the total number of modes that will
+% be used or a 3xN matrix where the first
+% dimension will contain the mode index along
+% each dimension and N corresponds to the
+% number of total modes. See notes for the
+% mode indices used when a scalar is provided.
+% [Default: 16^3 = 4096].
+%
+% dampRatio [numeric] (Optional): The damping ratio of the boundaries. This
+% can be a scalar denoting that the damping
+% ratio is identical for all frequencies or
+% a vector with the number of elements
+% matching the length of the "f" argument.
+% The damping ratio at the natural
+% frequencies of the room will be
+% calculated through linear interpolation
+% from the damping ratio values provided.
+% [Default: 0.1].
+%
+% Q [numeric] (Optional): The source strengths. This argument can be a
+% scalar indicating that all sources have the same
+% strength per frequency, a vector with N elements,
+% each corresponding to the source of each source
+% per frequency, or a matrix with dimensions NxF,
+% where each element will correspond to the
+% strength of each source for each frequency.
+% [Default: 1].
+%
+% c [numeric] (Optional): The speed of sound in m/s. [Default: 343].
+%
+% rho [numeric (Optional): The density of the medium in kg/(m^3).
+% [Default: 1.204].
+%
+% --------------------------------------------------
+% Output
+%
+% modeIdx [numeric]: The mode of the indices in the order they are placed
+% in the output arguments below (this helps to identify
+% contribution of each mode). This is a 3xK matrix with
+% the first dimension denoting the indices of each mode
+% and K denoting the number of modes.
+%
+% waveField [numeric]: The complex pressures at the position of the
+% receivers. The variable is an array of dimensions
+% MxNxFxK where M is the number of receivers, N the
+% number of sources and F the number of frequencies.
+%
+% srcWaveField [numeric]: The complex pressures at the positions of the
+% receivers with the contribution of the individual
+% modes being summed up. This should correspond to
+% the contribution of each source to each receiving
+% position per frequency. This is an MxNxF array.
+%
+% broadWaveField [numeric]: The complex pressures at the positions of the
+% receivers per mode with the contribution per
+% frequency being summed up. This corresponds to
+% the broadband (in the frequencies provided)
+% contribution of each mode at each position per
+% source. This is an MxNxK array.
+%
+% srcBroadWaveField [numeric]: The complex pressures at the positions of
+% the receivers with the contribution of each
+% mode and frequency summed up. This is an MxN
+% array and corresponds to the "total" (from a
+% frequency perspective) transfer function
+% from each source to each receiver
+%
+% sumWaveField [numeric]: The complex pressures at the positions of the
+% receivers with the contribution of the sources
+% being summed up. This corresponds to the sum of
+% the contributions of each source to each mode per
+% frequency and is an MxFxK array.
+%
+% sumSrcWaveField [numeric]: The complex pressures at the positions of the
+% receivers with the contributions of the
+% sources and the modes being summed up. This
+% corresponds to the total complex pressure at
+% each receiver position per frequency due to
+% all the sources and modes and is an MxF array.
+%
+% sumBroadWaveField [numeric]: The complex pressures at the positions of
+% the receivers with the contributions of the
+% sources summed up for all frequencies. This
+% correspond to the broadband (in the
+% frequencies provided) contribution of each
+% mode due to all sources and is an MxK array.
+%
+% totalWaveField [numeric]: This is the wavefield at each receiver position
+% with the contribution of each mode, source and
+% frequency summed up. In a sense this
+% corresponds to the total complex pressure at
+% each receiver position and is a vector with M
+% elements.
+%
+% --------------------------------------------------
+% Notes
+%
+% The actual number of modes calculated when a scalar is provided for
+% "nModes" is equal to ceil(nthroot(nModes, 3))^3 - 1. This value will
+% result in a set of mode indices that will provide all available
+% combinations of mode indices up to ceil(nthroot(nModes, 3)). As an
+% example, if nModes = 5 we would get ceil(nthroot(5, 3))^3 - 1 == 7 and
+% the modes used would be
+%
+% (1, 0, 0) | (0, 1, 0) | (0, 0, 1)
+% (1, 1, 0) | (1, 0, 1) | (0, 1, 1)
+% (1, 1, 1)
+%
+%
+%
+% Dependencies: - twoPtDist(): To calculate the distances between sources
+% and receivers.
+% --------------------------------------------------
+function [modeIdx, waveField, srcWaveField, broadWaveField, srcBroadWaveField, sumWaveField, sumSrcWaveField, sumBroadWaveField, totalWaveField] = rectRoomPtSrc(sPos, rPos, roomDims, f, nModes, dampRatio, Q, c, rho)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(4, 9);
+ nargoutchk(0, 9);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(sPos, "numeric", {'2d', 'nrows', 3, 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Source positions", 1);
+ validateattributes(rPos, "numeric", {'2d', 'nrows', 3, 'finite', 'nonnan', 'nonempty', 'real'}, mfilename, "Receiver positions", 2);
+ validateattributes(roomDims, "numeric", {'vector', 'finite', 'nonnan', 'nonempty', 'real', 'positive'}, mfilename, "Room dimensions", 3);
+
+ if ~isscalar(roomDims) && numel(roomDims) ~= 3
+ error("The dimensions argument must either be a scalar or a vector with three positive elements.");
+ end
+
+ validateattributes(f, "numeric", {'vector', 'finite', 'nonnan', 'nonempty', 'real', 'positive'}, mfilename, "Frequencies", 4);
+
+ % Validate optional arguments (and set their values)
+ if nargin > 8 && ~isempty(rho)
+ validateattributes(rho, "numeric", {'scalar', 'finite', 'positive', 'nonnan', 'nonempty', 'real'}, mfilename, "Medium density", 9);
+ else
+ rho = 1.204;
+ end
+
+ if nargin > 7 && ~isempty(c)
+ validateattributes(c, "numeric", {'scalar', 'finite', 'positive', 'nonnan', 'nonempty', 'real'}, mfilename, "Speed of propagation", 8);
+ else
+ c = 343;
+ end
+
+ if nargin > 6 && ~isempty(Q)
+ if isscalar(Q)
+ validateattributes(Q, "numeric", {'scalar', 'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
+ elseif isvector(Q)
+ validateattributes(Q, "numeric", {'vector', 'numel', size(sPos, 2), 'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
+ else
+ validateattributes(Q, "numeric", {'2d', 'size', [size(sPos, 2), length(f)] ,'nonempty', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Source strengths", 7);
+ end
+ else
+ Q = 1;
+ end
+
+ if nargin > 5 && ~isempty(dampRatio)
+ if isscalar(dampRatio)
+ validateattributes(dampRatio, "numeric", {'scalar', 'finite', 'nonempty', '<=', 1, 'real', 'nonnan', 'positive'}, mfilename, "Damping ratio", 6);
+ else
+ validateattributes(dampRatio, "numeric", {'vector', 'numel', length(f), 'finite', 'nonempty', '<=', 1, 'real', 'nonnan', 'positive'}, mfilename, "Damping ratio", 6);
+ end
+ else
+ dampRatio = 0.1;
+ end
+
+ if nargin > 4 && ~isempty(nModes)
+ if ~isscalar(nModes)
+ validateattributes(nModes, "numeric", {'2d', 'nrows', 3, 'nonempty', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Modes", 5);
+ else
+ validateattributes(nModes, "numeric", {'scalar', 'nonempty', 'integer', 'real', 'finite', 'nonnan', 'nonnegative'}, mfilename, "Number of modes", 5);
+ end
+ else
+ nModes = 4096;
+ end
+
+ % Check validity of source and receiver positions
+ if sum(~isbetween(sPos(1, :), -roomDims(1)/2, roomDims(1)/2)) > 0 || ...
+ sum(~isbetween(sPos(2, :), -roomDims(2)/2, roomDims(2)/2)) > 0 || ...
+ sum(~isbetween(sPos(3, :), -roomDims(3)/2, roomDims(3)/2)) > 0
+ error("Sources outside the room where found.");
+ end
+
+ if sum(~isbetween(rPos(1, :), -roomDims(1)/2, roomDims(1)/2)) > 0 || ...
+ sum(~isbetween(rPos(2, :), -roomDims(2)/2, roomDims(2)/2)) > 0 || ...
+ sum(~isbetween(rPos(3, :), -roomDims(3)/2, roomDims(3)/2)) > 0
+ error("Receivers outside the room where found.");
+ end
+
+
+ % ====================================================
+ % Pre-process
+ % ====================================================
+ % Convert frequencies to circular frequencies
+ f = 2 * pi * f;
+
+ % Create room dimensions vector
+ if isscalar(roomDims)
+ roomDims = roomDims * ones(3, 1);
+ end
+
+ % Create matrix with mode indices
+ if isscalar(nModes)
+ modeIdx = ceil(nthroot(nModes, 3)) - 1;
+ modeIdx = combinations(0:modeIdx, 0:modeIdx, 0:modeIdx);
+ modeIdx = table2array(modeIdx);
+ modeIdx = modeIdx(2:end, :);
+ else
+ modeIdx = nModes;
+ end
+ modeIdx = modeIdx./roomDims; % Calculate the ratio of the index to the dimensions of the room (used in the calculations)
+
+ % Calculate normalisation factors
+ normFact = sqrt(prod((modeIdx > 0) + 1, 2));
+
+ % Create vector with damping ratios
+ if isscalar(dampRatio)
+ dampRatio = dampRatio * ones(length(f), 1);
+ end
+
+ % Create matrix with source strengths
+ if isscalar(Q)
+ Q = Q * ones(size(sPos, 2), length(f));
+ elseif isvector(Q)
+ Q = repmat(Q(:), 1, length(f));
+ end
+
+
+ % ====================================================
+ % Calculate sound field
+ % ====================================================
+ % Calculate eigenfrequencies
+ wEig = pi * c * sqrt(sum(modeIdx.^2, 2));
+
+ % Inter/extrapolate damping ratio at the eigenfrequencies
+ if length(f) > 1
+ dampRatio = interp1(f, dampRatio, wEig, "linear", "extrap");
+ end
+
+ % Calculate the complex resonance terms
+ A = repmat(f(:).', length(wEig), 1)./ ...
+ (2 * dampRatio .* wEig .* f(:).' - ...
+ 1j * (repmat(wEig.^2, 1, length(f)) - repmat((f(:).').^2, length(wEig), 1)));
+ A = A * (rho * c^2)/prod(roomDims); % Multiply with the constant factor
+
+ % Calculate the eigenfunction values at the source and receiver positions
+ modeIdx = modeIdx * pi;
+
+ % Loop through frequencies
+ for fIdx = 1:length(f)
+ tmpA = A(:, fIdx);
+
+ % Loop through the receivers
+ for rIdx = size(rPos, 2):-1:1
+ rPsi = normFact .* prod(cos(modeIdx .* rPos(:, rIdx).'), 2);
+
+ % Loop through the sources
+ parfor sIdx = 1:size(sPos, 2)
+ waveField(rIdx, sIdx, fIdx, :) = rPsi .* tmpA .* normFact .* prod(cos(modeIdx .* sPos(:, sIdx).'), 2) * Q(sIdx, fIdx);
+ end
+ end
+ end
+
+ % ====================================================
+ % Return additional output arguments
+ % ====================================================
+ if nargout > 2
+ srcWaveField = squeeze(sum(waveField, 4));
+ end
+
+ if nargout > 3
+ broadWaveField = squeeze(sum(waveField, 3));
+ end
+
+ if nargout > 4
+ srcBroadWaveField = squeeze(sum(broadWaveField, 3));
+ end
+
+ if nargout > 5
+ sumWaveField = squeeze(sum(waveField, 2));
+ end
+
+ if nargout > 6
+ sumSrcWaveField = squeeze(sum(sumWaveField, 3));
+ end
+
+ if nargout > 7
+ sumBroadWaveField = squeeze(sum(sumWaveField, 2));
+ end
+
+ if nargout > 8
+ totalWaveField = squeeze(sum(sumBroadWaveField, 2));
+ end
+end