function varargout = ita_roomacoustics_analytic_FRF_pdi(varargin)
%ITA_ROOMACOUSTICS_ANALYTIC_FRF - After Kuttruff Room Acoustics pp.66
% This function calculates the ideal FRF of an rectangular room of size
% (Lx, Ly, Lz) for source position (rs_x, rs_y, rs_z) and receiver
% positions ... PLUS: cartesian differentian of the source or receiver
%
% Syntax:
% audioObjOut = ita_roomacoustics_analytic_FRF_differentiation(geometryCoordinates, sourceCoordinates, receiverCoordinates, options)
%
% Options (default):
% 'T' (2) : description
% 'f_max' (10000) : description
% 'fft_degree' (18) : description
% 'sourceDiff' ([0 0 0]): monopole, dipole in x ([1 0 0]), etc.
% 'receiverDiff' ([0 0 0]): monopole, dipole in x ([1 0 0]), etc.
%
% Example:
% audioObjOut = ita_roomacoustics_analytic_FRF(audioObjIn)
%
% See also:
% ita_toolbox_gui, ita_read, ita_write, ita_generate
%
% Reference page in Help browser
% doc ita_roomacoustics_analytic_FRF
%
% This file is part of the application RoomAcoustics for the ITA-Toolbox. All rights reserved.
% You can find the license for this m-file in the application folder.
%
% Author: Pascal -- Email: pdi@akustik.rwth-aachen.de
% Created: 13-Jul-2011
%% Initialization and Input Parsing
sArgs = struct('pos1_data','itaCoordinates', 'pos2_data','itaCoordinates','pos3_data','itaCoordinates',...
'fft_degree',18,'T',2,'f_max',10000,'c',ita_constants('c'),'rho_0',ita_constants('rho_0'),...
'receiverDiff',[0 0 0], 'sourceDiff', [0 0 0] );
[geometry, source_pos, receiver_pos ,sArgs] = ita_parse_arguments(sArgs,varargin);
%% Inits
sArgs.f_max = sArgs.f_max;
c = double(sArgs.c);
delta_n_raw = 3*log(10)/sArgs.T;
L = geometry.cart;
r_source = source_pos.cart;
r_receiver = receiver_pos.cart; % TODO: pdi: what about multiple receiving points?
%% check loop
warning off %#ok
idx = 0;
for nx = 0:floor(2*sArgs.f_max/c * L(1))
for ny = 0:floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 ) * L(2) ))
if ~isempty(ny) && isreal(ny)
idx = idx + floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 - (ny/L(2))^2) * L(3))) + 1;
end
end
end
nModes = idx;
dummy = zeros(1,nModes);
%%
n = zeros(3,nModes);
idx = 1;
nx_max = floor(2*sArgs.f_max/c * L(1));
for nx = 0:nx_max
ny_max = floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 ) * L(2) ));
for ny = 0:ny_max
nz_max = floor( real(sqrt( (2*sArgs.f_max/c)^2 - (nx/L(1))^2 - (ny/L(2))^2) * L(3)));
idx_end = idx + nz_max;
n(1,idx:idx_end) = nx;
n(2,idx:idx_end) = ny;
n(3,idx:idx_end) = 0:nz_max;
idx = idx + nz_max + 1;
end
end
ita_verbose_info([num2str(idx - nz_max - 1) ' modes have been calculated...'],1);
%% Eigenfrequencies
f_n = c / 2 * sqrt( sum ( bsxfun(@rdivide,n,L').^2 ,1) );
%% Source Receiver Coefficients
coeff_r = potential_diff(r_receiver,L,n,sArgs.receiverDiff);
coeff_s = potential_diff(r_source,L,n,sArgs.sourceDiff);
%% new with SourceFactor [kg/s^2]
K_n = itaValue(prod(L),'m3') * 0.5.^(sum(n > 0,1));
factor = sArgs.c^2 / K_n /itaValue('Hz2');
coeff = coeff_r .* coeff_s .* double(factor).'; %unit = [];
corr = 2; % why 2??? %pdi:HUHU: not totally correct, should be j*omega after modal superposition
coeff = coeff ./ corr;
coeff(corr == 0) = 0; %avoid NaN for strange correction coefficient (zero)
if isa(delta_n_raw,'itaSuper')
delta = delta_n_raw.freq2value(f_n).';
else
delta = dummy + delta_n_raw;
end
delta = -abs(delta);
%% Set Output
if nargout == 3
varargout(1) = {f_n};
varargout(2) = {coeff};
varargout(3) = {delta};
else
pz = itaPZ;
pz.f = f_n;
pz.C = coeff;
pz.sigma = delta;
pz.exp_s = -1;
% pdi: currently it seems to be s/m !!! Jul 2011 - this is probably due
% to the douple pole handling here and the single pole handling with
% mirrowing positive poles in the itaPZ.freqresp
pz.unit = itaValue('1/m') * itaValue('1/m') ^ (sum(sArgs.receiverDiff) + sum(sArgs.sourceDiff));
varargout{1} = pz;
end
%end function
end
function res = potential_diff(pos,L,n,derivation)
% pos: coordinates (not normalized)
% L: geometry of room, length in meters
% n: eigennumber
% derivation: [0 0 0] - monopole, [1 0 0] - first dipole
res = 1;
for idx = 1:3
if mod(derivation(idx),2) % 1 3 5...
res_part = sin(pi * n(idx,:) * pos(idx)/L(idx));
else % 0 2 4
res_part = cos(pi * n(idx,:) * pos(idx)/L(idx));
end
res_part = res_part .* (pi * n(idx,:) / L(idx)).^derivation(idx) * (-1).^(1+mod(floor((derivation(idx)-1)/2),2));
res = res.*res_part;
end
end