function status = test_modal_weighting(modus)
%TEST_MODAL_WEIGHTING evaluates the modal weighting implementation
%
%   Usage: status = test_modal_weighting(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual
%
%   Output parameters:
%       status  - true or false
%
%   TEST_MODAL_WEIGHTING(modus) tests the implementation of the modal weighting
%   window for different local sound field synthesis techniques.

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% 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.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************


status = false;


%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);


%% ===== Configuration ==================================================
window_type = {...
    'rect'; ...
    'max-rE'; ...
    'kaiser'; ...
    'tukey'; ...
    };
% For plotting
conf = SFS_config;
conf.dimension = '2.5D';
conf.secondary_sources.number = 16;
conf.secondary_sources.geometry = 'circular';
conf.secondary_sources.size = 3;
conf.resolution = 300;
conf.plot.usedb = true;
conf.t0 = 'source';
X = [-2 2];
Y = [-2 2];
Z = 0;
src = 'pw';


%% ===== Computation and Plotting ========================================
if modus % graphical mode
    for window = window_type'
        conf.modal_window = window{:};
        update_title = false;
        base_title = ['Window: ' window{:}];
        switch conf.modal_window
        case {'rect', 'max-rE'}
            parameter = 1;
        otherwise
            parameter = [0 0.5 1];
            update_title = true;
        end
        for pp = 1:length(parameter)
            conf.modal_window_parameter = parameter(pp);
            sound_field_imp_nfchoa(X,Y,Z,[0 -1 0],'pw',0,conf);
            if update_title
                figure_title = [base_title ...
                                ', Parameter: ' ...
                                num2str(parameter(pp))];
            else
                figure_title = base_title;
            end
            title(figure_title);
        end
    end
else % numerical mode
    for window = window_type'
        conf.modal_window = window{:};
        conf.modal_window_parameter = 0.5;
        for order = [1 2]
            for ndtft = [2*order+1 3*order+1]
                [win,Win,Phi] = ...
                    modal_weighting(order,ndtft,conf);
                [win_ref,Win_ref,Phi_ref] = ...
                    ref_modal_weighting(order,ndtft,conf);
                if any([abs(win-win_ref) ...
                        abs(Win-Win_ref) ...
                        abs(Phi-Phi_ref)])>1e14
                    return
                end
            end
        end
    end
end


status = true;

end


function [win,Win,Phi] = ref_modal_weighting(order,ndtft,conf)
    % Returns reference values for some orders as calculated by the Matlab
    % implementation of commit ...
    switch conf.modal_window
    case 'rect'
        if order==1 && ndtft==2*order+1
            win = [1 1];
            Win = [1 0 0];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 1];
            Win = [0.75 0.25 -0.25 0.25];
            Phi = [0 1.570796326794897 3.141592653589793 4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 1 1];
            Win = [1 0 0 0 0];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 1 1];
            Win = [0.714285714285714 ...
                   0.257419676543548 ...
                  -0.178139943388210 ...
                   0.063577409701804 ...
                   0.063577409701804 ...
                  -0.178139943388210 ...
                   0.257419676543548];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case 'max-rE'
        if order==1 && ndtft==2*order+1
            win = [1 0.707106781186548];
            Win = [0.804737854124365 0.097631072937817 0.097631072937817];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.707106781186548];
            Win = [0.603553390593274 ...
                   0.250000000000000 ...
                  -0.103553390593274 ...
                   0.250000000000000];
            Phi = [0 ...
                   1.570796326794897 ...
                   3.141592653589793 ...
                   4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 0.866025403784439 0.500000000000000];
            Win = [0.746410161513775 ...
                   0.145243228056937 ...
                  -0.018448308813825 ...
                  -0.018448308813825 ...
                   0.145243228056937];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 0.866025403784439 0.500000000000000];
            Win = [0.533150115366983 ...
                   0.265342154409152 ...
                  -0.040912347323205 ...
                   0.008995135230562 ...
                   0.008995135230562 ...
                  -0.040912347323205 ...
                   0.265342154409152];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case {'kaiser', 'kaiser-bessel'}
        % conf.modal_window_parameter == 0.5
        if order==1 && ndtft==2*order+1
            win = [1 0.581816879161903];
            Win = [0.721211252774602 0.139394373612699 0.139394373612699];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.581816879161903];
            Win = [0.540908439580952 ...
                   0.250000000000000 ...
                  -0.040908439580952 ...
                   0.250000000000000];
            Phi = [0 ...
                   1.570796326794897 ...
                   3.141592653589793 ...
                   4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 0.883766863470753 0.581816879161903];
            Win = [0.786233497053062 ...
                   0.120959694808693 ...
                  -0.014076443335224 ...
                  -0.014076443335224 ...
                   0.120959694808693];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 0.883766863470753 0.581816879161903];
            Win = [0.561595355037902 ...
                   0.263300911786297 ...
                  -0.063101577944049 ...
                   0.019002988638802 ...
                   0.019002988638802 ...
                  -0.063101577944049 ...
                   0.263300911786297];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    case 'tukey'
        % conf.modal_window_parameter == 0.5
        if order==1 && ndtft==2*order+1
            win = [1 0.75];
            Win = [0.833333333333333 0.083333333333333 0.083333333333333];
            Phi = [0 2.094395102393196 4.188790204786391];
        elseif order==1 && ndtft==3*order+1
            win = [1 0.75];
            Win = [0.625 0.250 -0.125 0.250];
            Phi = [0 1.570796326794897 3.141592653589793 4.712388980384690];
        elseif order==2 && ndtft==2*order+1
            win = [1 1 0.5];
            Win = [0.800000000000000 ...
                   0.161803398874989 ...
                  -0.061803398874989 ...
                  -0.061803398874989 ...
                   0.161803398874989];
            Phi = [0 ...
                   1.256637061435917 ...
                   2.513274122871834 ...
                   3.769911184307752 ...
                   5.026548245743669];
        elseif order==2 && ndtft==3*order+1
            win = [1 1 0.5];
            Win = [0.571428571428571 ...
                   0.289208381394450 ...
                  -0.049430105116435 ...
                  -0.025492561992301 ...
                  -0.025492561992301 ...
                  -0.049430105116435 ...
                   0.289208381394450];
            Phi = [0 ...
                   0.897597901025655 ...
                   1.795195802051310 ...
                   2.692793703076966 ...
                   3.590391604102621 ...
                   4.487989505128276 ...
                   5.385587406153931];
        end
    end
end
