function status = test_sphbesselh_zeros(modus)
%TEST_SPHBESSELH_ZEROS tests the computation of the roots for spherical Bessel
%functions
%
%   Usage: status = test_sphbesselh_zeros(modus)
%
%   Input parameters:
%       modus   - 0: numerical
%                 1: visual
%
%   Output parameters:
%       status  - true or false

%*****************************************************************************
% 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, subiect 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);


%% ===== Main ============================================================
markers = {'.', '^', '*', 'o'};
linecolors = [
         0 0.4470 0.7410
    0.8500 0.3250 0.0980
    0.9290 0.6940 0.1250
    0.4940 0.1840 0.5560
    0.4660 0.6740 0.1880
    0.3010 0.7450 0.9330
    0.6350 0.0780 0.1840
];
subplotdims = [2,3];

Nsym = length(markers);
Nsubplots = prod(subplotdims);
Ncolors = size(linecolors,1);
step = 2;

% Check for warnings and plot zeros
if modus
    figure;
    for order=1:step:300
        % Some indices
        gdx = (order-1)./step;  % graph index
        sdx = mod(gdx,Nsubplots)+1;  % subfig index
        mdx = mod(floor(gdx./(Nsubplots*Ncolors)), Nsym)+1;  % marker index
        cdx = mod(floor(gdx./Nsubplots), Ncolors)+1;  % color index
        % Plotting
        [z,~] = sphbesselh_zeros(order);
        subplot(subplotdims(1),subplotdims(2),sdx);
        hold on;
        plot(real(z), imag(z), ...
            'Color', linecolors(cdx,:), ...
            'Marker', markers{mdx}, ...
            'DisplayName', sprintf('%d',order));
        hold off;
    end
    % Add titles, legends, etc.
    for sdx=1:Nsubplots
        subplot(subplotdims(1),subplotdims(2),sdx);
        xlim([min(real(z)),0]);
        ylim([min(imag(z)),max(imag(z))]);
        title_str = sprintf(['Zeros of spherical Bessel functions of ', ...
                             'orders %i:12:%i'],2*sdx-1,2*sdx-1+288);
        title(title_str);
    end
else
    for order=[1 10 100]
        z = sphbesselh_zeros(order);
        z_ref = ref_sphbesselh_zeros(order);
        if any(abs(z-z_ref)>1e8)
            return
        end
    end
end

status = true;

end

function z = ref_sphbesselh_zeros(N)
    % Return reference values for some orders as calculated by
    % the following Python code:
    %
    % from scipy import signal
    % for N in [1, 10, 100]:
    %   r = signal.filter_design._bessel_zeros(N)
    %   print(1/r)
    switch N
    case 1
        z = [-1-0i];
    case 10
        z = [-3.10891623+8.23269946i; -4.88621957+6.22498548i; ...
             -5.96752833+4.38494719i; -6.61529097+2.61156792i; ...
             -6.92204491+0.8676652i;  -6.92204491-0.8676652i; ...
             -6.61529097-2.61156792i; -5.96752833-4.38494719i; ...
             -4.88621957-6.22498548i; -3.10891623-8.23269946i];
    case 100
        z = [ -7.27958017+96.07487622i; -16.56840794+89.69305883i; ...
             -12.47552337+92.61711308i; -20.06460298+87.0388997i; ...
             -23.16077972+84.55529001i; -25.95970766+82.1917017i; ...
             -28.52378431+79.91820025i; -30.89456430+77.71529775i; ...
             -33.10152255+75.5694611i;  -35.16654671+73.47083873i; ...
             -37.10646671+71.41199594i; -38.93458277+69.38715996i; ...
             -40.66164008+67.39174367i; -43.84648113+63.47496561i; ...
             -46.71605927+59.63895796i; -45.31789403+61.54799302i; ...
             -48.04558910+57.74602015i; -49.31049835+55.86759414i; ...
             -50.51430689+54.00230256i; -51.66011999+52.14893987i; ...
             -42.29647848+65.42203167i; -52.75069221+50.30644364i; ...
             -53.78847888+48.47387175i; -54.77567785+46.65038396i; ...
             -55.71426383+44.83522681i; -56.60601678+43.02772133i; ...
             -57.45254558+41.22725279i; -58.25530796+39.43326208i; ...
             -59.01562729+37.64523857i; -59.73470691+35.86271397i; ...
             -60.41364233+34.08525709i; -61.05343172+32.31246944i; ...
             -61.65498500+30.54398132i; -62.21913160+28.77944847i; ...
             -62.74662735+27.01854918i; -63.23816028+25.26098169i; ...
             -63.69435590+23.50646194i; -64.11578160+21.75472159i; ...
             -64.50295064+20.00550618i; -64.85632558+18.25857355i; ...
             -65.17632127+16.5136924i;  -65.46330747+14.77064097i; ...
             -65.71761112+13.0292058i;  -65.93951828+11.28918071i; ...
             -66.12927583+9.5503657i;  -66.28709289+7.81256601i; ...
             -66.41314197+6.07559126i; -66.50755997+4.33925457i; ...
             -66.57044891+2.60337173i; -66.60187650+0.86776044i; ...
             -66.60187650-0.86776044i; -66.57044891-2.60337173i; ...
             -66.50755997-4.33925457i; -66.41314197-6.07559126i; ...
             -66.28709289-7.81256601i; -66.12927583-9.5503657i; ...
             -65.93951828-11.28918071i; -65.71761112-13.0292058i; ...
             -65.46330747-14.77064097i; -65.17632127-16.5136924i; ...
             -64.85632558-18.25857355i; -64.50295064-20.00550618i; ...
             -64.11578160-21.75472159i; -63.69435590-23.50646194i; ...
             -63.23816028-25.26098169i; -62.74662735-27.01854918i; ...
             -62.21913160-28.77944847i; -61.65498500-30.54398132i; ...
             -61.05343172-32.31246944i; -60.41364233-34.08525709i; ...
             -59.73470691-35.86271397i; -59.01562729-37.64523857i; ...
             -58.25530796-39.43326208i; -57.45254558-41.22725279i; ...
             -56.60601678-43.02772133i; -55.71426383-44.83522681i; ...
             -54.77567785-46.65038396i; -53.78847888-48.47387175i; ...
             -52.75069221-50.30644364i; -42.29647848-65.42203167i; ...
             -51.66011999-52.14893987i; -50.51430689-54.00230256i; ...
             -49.31049835-55.86759414i; -48.04558910-57.74602015i; ...
             -45.31789403-61.54799302i; -46.71605927-59.63895796i; ...
             -43.84648113-63.47496561i; -40.66164008-67.39174367i; ...
             -38.93458277-69.38715996i; -37.10646671-71.41199594i; ...
             -35.16654671-73.47083873i; -33.10152255-75.5694611i; ...
             -30.89456430-77.71529775i; -28.52378431-79.91820025i; ...
             -25.95970766-82.1917017i;  -23.16077972-84.55529001i; ...
             -20.06460298-87.0388997i;  -12.47552337-92.61711308i; ...
             -16.56840794-89.69305883i;  -7.27958017-96.07487622i];
    end
end
