retarded 3D Green's function uniformly moving point source

SFS_general/direction_vector.m

@@ -1,7 +1,7 @@
 function directions = direction_vector(x1,x2)
 %DIRECTION_VECTOR return unit vector(s) pointing from x1 to x2
 %
-%   Usage: n = direction_vector(x1,x2)
+%   Usage: n = direction_vector(x1,[x2])
 %
 %   Input parameters:
 %       x1  - starting point(s) [1xn] or [mxn]
@@ -11,7 +11,8 @@
 %       n   - unit vector(s) pointing in the direction(s) from x1 to x2
 %
 %   DIRECTION_VECTOR(x1,x2) calculates the unit vectors pointing from
-%   n-dimensional points x1 to the n-dimensional points x2.
+%   n-dimensional points x1 to the n-dimensional points x2. If only x1 is
+%   provided the direction vectors from zero to x1 are calculated.
 %
 %   See also: secondary_source_positions
 
@@ -46,15 +47,20 @@
 
 
 %% ===== Checking of input  parameters ===================================
-nargmin = 2;
+nargmin = 1;
 nargmax = 2;
 narginchk(nargmin,nargmax);
-if size(x1,2)~=size(x2,2)
-    error('%s: x1 and x2 need to have the same dimension.',upper(mfilename));
-end
-if size(x1,1)~=size(x2,1) && ~(size(x1,1)==1 | size(x2,1)==1)
-    error(['%s: x1 and x2 need to have the same size, or one needs to ', ...
-           'be a vector.'],upper(mfilename));
+if nargin==1
+    x2 = x1;
+    x1 = zeros(1,size(x2,2));
+else
+    if size(x1,2)~=size(x2,2)
+        error('%s: x1 and x2 need to have the same dimension.',upper(mfilename));
+    end
+    if size(x1,1)~=size(x2,1) && ~(size(x1,1)==1 | size(x2,1)==1)
+        error(['%s: x1 and x2 need to have the same size, or one needs to ', ...
+               'be a vector.'],upper(mfilename));
+    end
 end
 
 

SFS_general/retarded_time.m

@@ -0,0 +1,90 @@
+function [tau, R, Delta] = retarded_time(x, y, z, t, vs, conf)
+% RETARDED_TIME computes the retarded time for a uniformly moving point source
+%
+%   Usage: [tau, R, Delta] = retarded_time(x, y, z, t, vs, conf)
+%
+%   Input parameters:
+%       x           - x coordinate position of observer / m [nx1] or [1x1]
+%       y           - y coordinate position of observer / m [nx1] or [1x1]
+%       z           - z coordinate position of observer / m [nx1] or [1x1]
+%       t           - absolute time / s [nxm] or [1xm]
+%       vs          - velocity vector of sound source / (m/s) [1x3] or [nx3]
+%       conf        - configuration struct (see SFS_config)
+%
+%   Output parameters:
+%       tau         - retarded time between sound source and observer / s [nxm]
+%       R           - retarded distance (R = tau*c) / m [nxm]
+%       Delta       - auxilary distance / m [nxm]
+%
+%   RETARDED_TIMES(x, y, z, t, vs, conf) the retarded time for a point source
+%   uniformly moving along a line. Hereby, vs defines its direction
+%   and its velocity. The position of the sound source at t=0 is the coordinate
+%   origin.
+%
+%   See also:
+
+%*****************************************************************************
+% The MIT License (MIT)                                                      *
+%                                                                            *
+% Copyright (c) 2010-2016 SFS Toolbox Team                                   *
+%                                                                            *
+% Permission is hereby granted,  free of charge,  to any person  obtaining a *
+% copy of this software and associated documentation files (the "Software"), *
+% to deal in the Software without  restriction, including without limitation *
+% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
+% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
+% Software is furnished to do so, subject to the following conditions:       *
+%                                                                            *
+% The above copyright notice and this permission notice shall be included in *
+% all copies or substantial portions of the Software.                        *
+%                                                                            *
+% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
+% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
+% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
+% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
+% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
+% DEALINGS IN THE SOFTWARE.                                                  *
+%                                                                            *
+% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
+% methods like wave field synthesis or higher order ambisonics.              *
+%                                                                            *
+% http://sfstoolbox.org                                 [email protected] *
+%*****************************************************************************
+
+%% ===== Checking of input  parameters ==================================
+% nargmin = 6;
+% nargmax = 6;
+% narginchk(nargmin,nargmax);
+% isargmatrix(x,y,z,t,vs);
+% isargstruct(conf);
+
+
+%% ===== Configuration ==================================================
+c = conf.c;
+
+
+%% ===== Computation ====================================================
+
+v = vector_norm(vs,2);  % velocity of sound source [n x 1] or [1 x 1]
+nvs = bsxfun(@rdivide,vs,v);  % direction of movement [n x 3] or [1 x 3]
+M = v/c;  % Mach number [n x 1] or [1 x 1]
+
+% vector between observer and point source d(t) = x - vs*t
+% [n x m] or [1 x m]
+dx = bsxfun(@minus, x, bsxfun(@times,t, vs(:,1)));
+dy = bsxfun(@minus, y, bsxfun(@times,t, vs(:,2)));
+dz = bsxfun(@minus, z, bsxfun(@times,t, vs(:,3)));
+% instant distance between sound source and observer position
+r = sqrt( dx.^2 + dy.^2 + dz.^2 );  % [n x m]
+% component of d(t) parallel to the direction of movement: d(t)^T * nvs
+% [n x m]
+xpara = nvs(:,1).*dx + nvs(:,2).*dy + nvs(:,3).*dz;
+% auxiliary distance  
+Delta = sqrt( M.^2.*xpara.^2 + (1-M.^2).*r.^2 );  % [n x 1]
+% retarded distance
+R = (M.*xpara + Delta)./(1-M.^2);
+% retarded time
+tau = R./c;
+
+end

SFS_helper/isargsecondarysource.m

@@ -54,9 +54,9 @@ function isargsecondarysource(varargin)
             'and weights.'], ...
             inputname(ii));
     end
-    if ~all(abs(vector_norm(x0(:,4:6),2)-1)<1e-10)
-        error(['The norm of the direction of %s is not 1. ', ...
-            'You can use the direction_vector function to get ', ...
-            'correct direction values.'],inputname(ii));
-    end
+%     if ~all(abs(vector_norm(x0(:,4:6),2)-1)<1e-10)
+%         error(['The norm of the direction of %s is not 1. ', ...
+%             'You can use the direction_vector function to get ', ...
+%             'correct direction values.'],inputname(ii));
+%     end
 end

SFS_monochromatic/greens_function_mono.m

@@ -25,6 +25,7 @@
 %   References:
 %       H. Wierstorf, J. Ahrens, F. Winter, F. Schultz, S. Spors (2015) -
 %       "Theory of Sound Field Synthesis"
+%       J. Ahrens (2012) - "Analytic Methods of Sound Field Synthesis"
 %
 %   See also: sound_field_mono
 
@@ -99,6 +100,32 @@
     %
     G = 1/(4*pi) .* (1i*omega/c + 1./r) .* scalar./r.^2 .* exp(-1i*omega/c.*r);
 
+elseif strcmp('mps',src)
+    % Source model for a uniformly moving point source: retarded 3D Green's 
+    % function.
+    %
+    %                 1   e^(+-i w/c R)
+    % G(x-xs,vs,w) = --- --------------
+    %                4pi      R_1
+    %
+    % vs = v * ns   % velocity vector
+    % M = v/c  % Mach number
+    % xparallel = (x-xs)*ns  % component in direction of movement
+    % R_1 = sqrt( M^2*xparallel.^2 + (1-M^2)*|x-xs|.^2 )
+    % R = (M*xparallel + R1)./(1-M*2);
+    %    
+    % see: Ahrens (2012), eq.(5.60)
+    %    
+    [~,x1]  = xyz_axes_selection(x,y,z); % get first non-singleton axis
+    % shift and vectorize coordinates
+    x = x(:)-xs(1);
+    y = y(:)-xs(2);
+    z = z(:)-xs(3);
+    % retarded time
+    [tau, ~, R1] = retarded_time(x,y,z,0,xs(4:6),conf);
+    %
+    G = zeros(size(x1));
+    G(:) = 1/(4*pi) * exp(-1i*omega.*tau)./R1;
 elseif strcmp('ls',src)
     % Source model for a line source: 2D Green's function.
     %