function kettledrum()         % S.W. Rienstra 2021-02-25
   clear all;
   warning off MATLAB:MKDIR:DirectoryExists
   format long
   format compact

   %% parameters
   a     = 0.656/2    ; % radius of cavity
   L     = 0.4142     ; % length of cavity
   d     = 0.028/2    ; % radius hole
   rhoa  = 1.21       ; % density air
   ca    = 344        ; % sound speed air
   T     = 3990       ; % membrane tension
   sgm   = 0.2653     ; % membrane density
   cM    = sqrt(T/sgm); % sound speed membrane
   V     = pi*a^2*L   ; % volume cavity
   mmax  = 7      ;     % number of m-orders
   frqmx = 1000   ;     % maximum frequency considered for border-case resonance frequencies: without cavity and without membrane
   Bmax  = 1200;        % estimated maximum values of Bessel zeros x and y
   xmax  = 500;         % number of x and y
   Cmax  = 50;          % summation of C(m,n)
   Smax  = 50;          % summation of S(n)
   Nmax  = 15;          % size of matrix. Consider only first few eigenvalues. Otherwise increase Nmax
   jmax  = 50;          % number of k-iterations. Maybe too high. Check output
   eta   = 0.5;         % moderating parameter. If convergent, eta=0 gives superlinear convergence, but especially for m=0 not always.
   idealmodes    = 0;
   integraaltest = 0;
   testoutput    = 0;
   shotofhail    = 0;
   %
   if shotofhail
      for m=0:mmax
        m1 = m+1;
        kstart(m1,:)   = [100:10:1000]*2*pi*a/ca;  % do this first, and then select suitable starting values
      end
   else
      kstart(0+1,:)   = [125, 253, 416,  477, 614, 776, 808,  986]*2*pi*a/ca;
      kstart(1+1,:)   = [151, 311, 351,  507, 567, 707, 879, 1080]*2*pi*a/ca;
      kstart(2+1,:)   = [228, 412, 600,  795, 969,   0,   0,    0]*2*pi*a/ca;
      kstart(3+1,:)   = [300, 493, 678, 1060,   0,   0,   0,    0]*2*pi*a/ca;
      kstart(4+1,:)   = [370, 570, 763,  952,   0,   0,   0,    0]*2*pi*a/ca;
      kstart(5+1,:)   = [439, 646, 843, 1033,   0,   0,   0,    0]*2*pi*a/ca;
      kstart(6+1,:)   = [508, 720, 920,    0,   0,   0,   0,    0]*2*pi*a/ca;
      kstart(7+1,:)   = [575, 794, 997,    0,   0,   0,   0,    0]*2*pi*a/ca;
   end
   for m=0:mmax
       m1 = m+1;
       Nstart(m1) = find(kstart(m1,:)~=0,1,'last');
   end
   %this is probably not necessary:
   %x     = ones(mmax+1,xmax)*NaN;
   %y     = ones(mmax+1,xmax)*NaN;
   %gmn   = ones(mmax+1,xmax)*NaN;
   %% make x_mn
   for m=0:mmax
     m1       = m+1;
     j0       = [0.01:0.1:(Bmax+m/2-1/4)*pi/2];                                       % fine grid initial values
     jm0      = j0(find(abs(diff(J(m,j0)>=0))));                                      % find the sign changes >0 and add to initial values jm0
     jm       = fsolve(@(x) J(m,x),jm0,optimoptions('fsolve','Display','off'));
     x(m1,1:xmax) = jm(1:xmax);
   end
   %% make y_mn
   for m=0:mmax
     m1       = m+1;
     dj0      = [0.01:0.1:(Bmax+m/2-3/4)*pi/2];                                       % fine grid initial values
     djm0     = dj0(find(abs(diff(dJ(m,dj0)>=0))));                                   % find the sign changes >0 and add to initial values djm0
     djm      = fsolve(@(x) dJ(m,x),djm0,optimoptions('fsolve','Display','off'));
     if m==0
           djm = [0,djm];                        % add 0 as first non-trivial zero of J_m' if m=0.
     end
     y(m1,1:xmax) = djm(1:xmax);
   end
   %% make om0_mn  (membrane resonances, i.e. without air loading)
   %% and  omc_mnl (cavity resonances, i.e. with infintely stiff membrane)
   for m=0:mmax
     m1       = m+1;
     N0(m1)   = find(x(m1,:)*cM/a>2*pi*frqmx,1);
     if ~isempty(N0(m1))
        om0(m1,1:N0(m1)) = x(m1,1:N0(m1))*cM/a;
     end
     Nc(m1)   = find(y(m1,:)*ca/a>2*pi*frqmx,1);
     if ~isempty(Nc(m1))
        for n=1:Nc(m1)
           if (2*L*frqmx/ca)^2 > (y(m1,n)*L/pi/a)^2
             for l=0:floor(sqrt((2*L*frqmx/ca)^2-(y(m1,n)*L/pi/a)^2))
                l1 = l+1;
                omc(m1,n,l1) = ca/a*sqrt(y(m1,n).^2 + (l*pi*a/L)^2);
             end
           end
        end
     end
   end
   if idealmodes
      om0
      omc
      return
   end
   % test of critical complex contour integral
   if integraaltest
      k    = 2.0005+0.5*i;
      m    = 0;
      m1   = m+1;
      n    = 2;
      n1   = 2;
      splt = [0:0.001:10];
      fplt = ftest(splt,k,m,n,n1,x);
      figure(1)
      plot(splt,real(fplt),splt,imag(fplt))
      disp([sprintf('  %0.8g,      ',ftest(10,k,m,n,n1,x)), sprintf('%0.8g',ftest(100,k,m,n,n1,x))])
      %sum(fplt)*0.001  % crude check of integral
      disp([integral(@(s) ftest(s,k,m,n,n1,x),0,10) + integral(@(s) ftest(s,k,m,n,n1,x),10,100)])  % 100 is usually enough for O(1E-12)
      if n==n1
         disp([findImn(2,k,m,n,x)])
      else
         disp([(findImn(1,k,m,n,x)-findImn(1,k,m,n1,x))/(x(m1,n)^2-x(m1,n1)^2)])
      end
      return
   end
   % here the iteration really starts.
   for m=0:mmax             % outer iteration over m
     m1 = m+1;
     for mm = 1:Nstart(m1)  % inner iteration over the resonace frequencies per m
       k  = kstart(m1,mm);
       if testoutput
         disp(sprintf('m=[%1i,%1i],          kstart=%15.12f,     frq=%10.5f',m,mm,k,k*ca/(2*pi*a)));
       end
       for j=1:jmax          % iteratie over k, per frequency
           for n =1:xmax
               gmn(m1,n) = gm(y(m1,n),k);
           end
           if m==0
              for n = 1:Nmax
                  S(n) = 0;
                  for nn=1:Smax
                      S(n) = S(n) + 1/((x(1,n)^2-y(1,nn)^2)*gmn(1,nn)*sin(gmn(1,nn)*L/a)*J(0,y(1,nn)));
                  end
              end
           end
           for n = 1:Nmax
               Cm1(m1,n) = 0;
               Cm2(m1,n) = 0;
               for nn=1:Cmax
                  if m==0
                     Cm1( 1,n) = Cm1( 1,n) + cot(gmn( 1,nn)*L/a)/gmn( 1,nn)/(y( 1,nn)^2-x( 1,n)^2);
                     Cm2( 1,n) = Cm2( 1,n) + cot(gmn( 1,nn)*L/a)/gmn( 1,nn)/(y( 1,nn)^2-x( 1,n)^2)^2;
                  else
                     Cm1(m1,n) = Cm1(m1,n) + cot(gmn(m1,nn)*L/a)/gmn(m1,nn)/(y(m1,nn)^2-x(m1,n)^2)/(1-(m/y(m1,nn))^2);
                     Cm2(m1,n) = Cm2(m1,n) + cot(gmn(m1,nn)*L/a)/gmn(m1,nn)/(y(m1,nn)^2-x(m1,n)^2)^2/(1-(m/y(m1,nn))^2);
                  end
               end
           end
           for n = 1:Nmax
               Im1(m1,n) = findImn(1,k,m,n,x);
               Im2(m1,n) = findImn(2,k,m,n,x);
           end
           for n = 1:Nmax        % this could be more efficient. use symmetry of A
             for n1 = 1:Nmax
                 if n1==n
                    if m==0
                       Am(n,n) = 1/x(1,n)^2 - 4*rhoa*a/sgm*(Cm2(1,n)-i/2*Im2(1,n) - 2*d/pi/a*S(n)^2);
                    else
                       Am(n,n) = 1/x(m1,n)^2 - 4*rhoa*a/sgm*(Cm2(m1,n)-i/2*Im2(m1,n));
                    end
                 else
                    if m==0
                       Am(n,n1) = - 4*rhoa*a/sgm*x(1,n1)/x(1,n)*((Cm1(1,n)-Cm1(1,n1) -i/2*(Im1(1,n)-Im1(1,n1)))/(x(1,n)^2-x(1,n1)^2) - 2*d/pi/a*S(n)*S(n1));
                    else
                       Am(n,n1) = - 4*rhoa*a/sgm*x(m1,n1)/x(m1,n)*((Cm1(m1,n)-Cm1(m1,n1) -i/2*(Im1(m1,n)-Im1(m1,n1)))/(x(m1,n)^2-x(m1,n1)^2));
                    end
                 end
             end
           end
           E     = eig(Am);
           [q,p] = min(abs(k-cM/ca./sqrt(E)));     % don't forget the dot
           k     = eta*k+(1-eta)*cM/ca/sqrt(E(p));
           if (testoutput || j==jmax)
             disp(sprintf('m=[%1i,%1i], q=%10.3E, k=%8.5f %+8.5fi,  frq=%10.5f, t60=%10.5f',m,mm, q, real(k), imag(k), real(k*ca/(2*pi*a)),-3/imag(k*ca/a)*log(10)))
           end
       end  % j-loop
       if testoutput
         disp(sprintf(' '))
       end
     end  % mm-loop
   end  % m-loop
   %
   %    End of program
   %
end
%% external auxiliary functions
function y = gm(zz,kk)
      y = i.*sqrt(i*(zz-kk)).*sqrt(-i*(zz+kk));
end
function y = findImn(r,k,m,n,x)
    m1 = m+1;
    Th0 = asin(i*real(k)/k);
    Th1 = asin(real(k)/k);
    Th2 = pi/2;
    switch r
      case 1
             y = 2*i/pi*integral(@(s) f10(s),0,1/real(k));      % integral(@(s) f(s),0,1,'AbsTol',1e-6,'RelTol',1e-3)    voor als je nauwkeurigheid wilt sturen
             if imag(k)<0
               y = y - 2/pi*(Th1-Th0)*integral(@(s) f11(s),0,1);
               y = y - 2/pi*(Th2-Th1)*integral(@(s) f12(s),0,1);
             else
               y = y - 2/pi*(Th2-Th0)*integral(@(s) f13(s),0,1);
             end
      case 2
             y = 1/(2*gm(x(m1,n),k)*x(m1,n)^2);
             y = y -i/pi*integral(@(s) f20(s),0,1/(real(k)^2));  % integral(@(s) f(s),0,1,'AbsTol',1e-6,'RelTol',1e-3)  voor als je nauwkeurigheid wilt sturen
             if imag(k)<0
               y = y + 2/pi*(Th1-Th0)*integral(@(s) f21(s),0,1);
               y = y + 2/pi*(Th2-Th1)*integral(@(s) f22(s),0,1);
             else
               y = y + 2/pi*(Th2-Th0)*integral(@(s) f23(s),0,1);
             end
    end

    function y = f10(s)
             s1 = [s~=0].*s;
             s2 = s1.*s1;
             y = [s==0].*0 + [s~=0].*IK(m,1./s1,1)./sqrt(1+k^2.*s2)./(1+x(m1,n)^2.*s2);
    end
    function y = f20(s)
             s1 = [s~=0].*s;
             s2 = sqrt(s1);
             y = [s==0].*0 + [s~=0].*s2.*IK(m,1./s2,1)./sqrt(1+k^2.*s1)./((1+x(m1,n)^2.*s1).^2);
    end
    function y = f11(h)
             G1h = G1(h);
             y = G1h.*IK(m,G1h,0)./(G1h.^2+x(m1,n)^2);
    end
    function y = f12(h)
             G2h = G2(h);
             y = G2h.*IK(m,G2h,0)./(G2h.^2+x(m1,n)^2);
    end
    function y = f13(h)
             G3h = G3(h);
             y = G3h.*IK(m,G3h,0)./(G3h.^2+x(m1,n)^2);
    end
    function y = f21(h)
             G1h = G1(h);
             y = G1h.*IK(m,G1h,0)./((G1h.^2+x(m1,n)^2).^2);
    end
    function y = f22(h)
             G2h = G2(h);
             y = G2h.*IK(m,G2h,0)./((G2h.^2+x(m1,n)^2).^2);
    end
    function y = f23(h)
             G3h = G3(h);
             y = G3h.*IK(m,G3h,0)./((G3h.^2+x(m1,n)^2).^2);
    end
    function y = G1(h)
             y = -i*k*sin(Th0+(Th1-Th0)*h);
    end
    function y = G2(h)
             y = -i*k*sin(Th1+(Th2-Th1)*h);
    end
    function y = G3(h)
             y = -i*k*sin(Th0+(Th2-Th0)*h);
    end
end
function y  = ftest(t,k,m,n,n1,x)
         m1 = m+1;
         y0 = J(m,t);
         y1 = y0./(t-x(m1,n));
         y2 = y0./(t-x(m1,n1));
         y  = t.*y1.*y2./(t + x(m1,n))./(t+x(m1,n1))./gm(t,k);
end
function y = dJ(n,x)
  if n==0
     y = -besselj(1,x);
  else
     y = (besselj(n-1,x)-besselj(n+1,x))/2;
  end
end
function y = J(n,x)
  y = besselj(n,x);
end
function y = IK(n,x,scale)                     % Note: I and K scale only the same if x is real positive
  y = besseli(n,x,scale).*besselk(n,x,scale);
end
