diff --git a/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
new file mode 100644
index 0000000000000000000000000000000000000000..909e7808f4a9562addb29fb337adf3682a49feb1
--- /dev/null
+++ b/Sound Fields/MATLAB/Functions/ptSrcFieldTD.m
@@ -0,0 +1,148 @@
+%% Signal generated by point sources in the time-domain
+% --------------------------------------------------
+% Author: Achilles Kappis
+% e-mail: axilleaz@protonmail.com
+%
+% Date: 13/09/2024 (DD/MM/YYYY)
+%
+% Copyright: MIT
+% --------------------------------------------------
+% Functionality: Calculate the signals generated by ideal point sources in
+% the time-domain.
+% --------------------------------------------------
+% Input
+%
+% sPos [numeric]: The Cartesian coordinates of the source positions. This
+% must be an Nx3 matrix where N is the number of sources.
+%
+% rPos [numeric]: The Cartesian coordinates of the receiver positions. This
+% must be an Mx3 matrix where M is the number of receivers.
+%
+% fs [numeric]: The sampling frequency. This is required in order to
+% calculate the time-of-flight in samples and must be a
+% positive real scalar.
+%
+% Q [numeric] (Optional): The source signals. This can be either an IxNxJ
+% array, where I is the length of the source
+% signals in samples and J is the number of
+% trials/sound field realisations, or a real
+% positive integer denoting the length of the
+% source signals. In case the length of the source
+% signals is given, the generated signals are
+% zero-mean, (approximately) uniformly distributed
+% and normalised to unity. [Default: 128].
+%
+% nTrials [numeric] (Optional): This is the number of trials/sound field
+% realisations that will be generated. If "Q"
+% is provided, this argument is ignored. It
+% must be a positive real scalar.
+% [Default: 1]
+%
+% c [numeric] (Optional): The speed of sound. [Default: 343].
+%
+% --------------------------------------------------
+% Output
+%
+% rSig [numeric]: The signals at the receiver positions for all
+% trials/sournd field realisations. This is an IxMxJ array.
+%
+% rSigMean [numeric]: The signals at the receiver positions averaged over
+% all trials/sound field realisations. This is an IxM
+% matrix.
+%
+% rSigMtx [numeric]: The signals at each receiver position due to each
+% source for each trial/sound field realisation. This is
+% an IxMxNxJ array.
+%
+% rSigMtxMean [numeric]: The signals at each receiver position due to each
+% source averaged over the trials/sound field
+% realisations. This is an IxMxN array.
+%
+% --------------------------------------------------
+% Notes
+%
+% - It is the responsibility of the user to pick (or provide) adequately
+% long source signals
+%
+% --------------------------------------------------
+function [rSig, rSigMean, rSigMtx, rSigMtxMean] = ptSrcFieldTD(sPos, rPos, fs, Q, nTrials, c)
+ % ====================================================
+ % Check for number of arguments
+ % ====================================================
+ narginchk(3, 6);
+ nargoutchk(0, 4);
+
+ % ====================================================
+ % Validate input arguments
+ % ====================================================
+ % Validate mandatory arguments
+ validateattributes(sPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the source positions", 1);
+ validateattributes(rPos, "numeric", {'2d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', 3}, mfilename, "Cartesian coordinates of the receiver positions", 2);
+ validateattributes(fs, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive'}, mfilename, "The sampling frequency", 3);
+
+ % Validate optional arguments
+ if nargin > 3 && ~isempty(Q)
+ if isscalar(Q)
+ validateattributes(Q, "numeric", {'scalar', 'real', 'nonnan', 'nonempty', 'finite', 'positive', 'integer'}, mfilename, "Length of source signals in samples", 4);
+ else
+ validateattributes(Q, "numeric", {'3d', 'real', 'nonnan', 'nonempty', 'finite', 'ncols', size(sPos, 1)}, mfilename, "Source signals", 4);
+ nTrials = size(Q, 3);
+ end
+ else
+ Q = 128;
+ end
+
+ if nargin > 4 && ~isempty(nTrials) && isscalar(Q)
+ validateattributes(nTrials, "numeric", {'scalar', 'positive', 'nonnan', 'nonempty', 'finite'}, mfilename, "Number of trials/sound field realisations", 5);
+ elseif nargin > 4 && isempty(nTrials)
+ nTrials = 1;
+ end
+
+ if nargin > 5 && ~isempty(c)
+ validateattributes(c, "numeric", {'scalar', 'nonempty', 'nonnan', 'finite', 'real', 'positive'}, mfilename, "The speed of sound", 6);
+ else
+ c = 343;
+ end
+
+
+ % ====================================================
+ % Pre-process data
+ % ====================================================
+ % Generate source signals
+ if isscalar(Q)
+ Q = rand(Q, size(sPos, 1), nTrials);
+ Q = Q - mean(Q);
+ Q = Q./max(abs(Q));
+ end
+
+ % Calculate source-receiver distances
+ dist = twoPtDist(sPos, rPos); % Distances
+ del = dist/c; % Delays in seconds
+
+ % ====================================================
+ % Calculate signals
+ % ====================================================
+ % Go through the trials/sound field realisations
+ for jIdx = size(Q, 3):-1:1
+ % Go through the sources (calculate for all receivers in one go)
+ for sIdx = size(dist, 1):-1:1
+ rSigMtx(:, :, sIdx, jIdx) = (1./dist(sIdx, :)) .* delayseq(Q(:, sIdx, jIdx), del(sIdx, :), fs);
+ end
+ end
+
+ % Sum source signals at each receiver position
+ rSig = squeeze(sum(rSigMtx, 3));
+
+ % ====================================================
+ % Calculate output arguments
+ % ====================================================
+ % Mean receiver signal over all trials/sound field realisations
+ if nargout > 1
+ rSigMean = squeeze(mean(rSig, 3));
+ end
+
+ % Mean receiver signal due to each source over all trials/sound field realisations
+ if nargout > 3
+ rSigMtxMean = squeeze(mean(rSigMtx, 4));
+ end
+end
\ No newline at end of file