function chladni_figures ( n, mu, k )

%*****************************************************************************80
%
%% chladni_figures() computes Chladni figures for a vibrating plate.
%
%  Discussion:
%
%    The printed program uses n = 100, mu = 0.225, k = 90.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2022
%
%  Author:
%
%    Original MATLAB version by Martin Gander, Felix Kwok.
%    This version by John Burkardt
%
%  Reference:
%
%    Martin Gander, Felix Kwok,
%    Chladni figures and the Tacoma bridge,
%    SIAM Review,
%    Volume 54, Number 3, pages 573-596, September 2012.
%
%  Input:
%
%    integer n: the number of nodes in the x and y directions,
%    There will be two additional ghost nodes, for a total of n+2.
%
%    real mu: a material constant, between 0 and 1.
%
%    integer k: the number of eigenmodes to compute.
%
  h = 2.0 / ( n - 3 );
%
%  Number the points in the grid.
%  'S' means the region is a square.
%  n+2 is the number of nodes.
%
  G = numgrid ( 'S', n + 2 );
%
%  Construct the sparse 5-point Laplacian on the given grid.
%
  D = delsq ( G );
%
%  Define the boundary points.
%
  bl = G(3:n,3);
  br = G(3:n,n);
  bt = G(3,3:n)';
  bb = G(n,3:n)';
%
%  Define ghost points.
%
  gl = G(3:n,2);
  gr = G(3:n,n+1);
  gt = G(2,3:n)';
  gb = G(n+1,3:n)';
%
%  Initialize N and L to the discrete Laplacian.
%
  L = D;
  N = D;
%
%  Correct the outer Laplacian for the boundary conditions.
%
  N(bl,bl) = N(bl,bl) / 2.0;
  N(br,br) = N(br,br) / 2.0;
  N(bt,bt) = N(bt,bt) / 2.0;
  N(bb,bb) = N(bb,bb) / 2.0;
%
%  Modify the stencil in L at the ghost points to approximate
%  ddu/dndt, which gives the correct boundary conditions.
%
  L([gl;gr;gt;gb],:) = 0.0;
%
%  Left:
%
  for i = gl(1:end-1)'
    L([i,i+1],[i,i+1,i+2*n,i+2*n+1] ) = ...
    L([i,i+1],[i,i+1,i+2*n,i+2*n+1] ) ...
    + 0.5 * ( mu - 1.0 ) * [1.0,-1.0,-1.0,1.0;-1.0,1.0,1.0,-1.0];
  end
%
%  Right
%
  for i = gr(1:end-1)'
    L([i,i+1],[i,i+1,i-2*n,i-2*n+1] ) = ...
    L([i,i+1],[i,i+1,i-2*n,i-2*n+1] ) ...
    + 0.5 * ( mu - 1.0 ) * [1.0,-1.0,-1.0,1.0;-1.0,1.0,1.0,-1.0];
  end
%
%  Top
%
  for i = gt(1:end-1)'
    L([i,i+n],[i+n,i,i+n+2,i+2] ) = ...
    L([i,i+n],[i+n,i,i+n+2,i+2] ) ...
    - 0.5 * ( mu - 1.0 ) * [1.0,-1.0,-1.0,1.0;-1.0,1.0,1.0,-1.0];
  end
%
%  Bottom
%
  for i = gb(1:end-1)'
    L([i,i+n],[i+n,i,i+n-2,i-2] ) = ...
    L([i,i+n],[i+n,i,i+n-2,i-2] ) ...
    - 0.5 * ( mu - 1.0 ) * [1.0,-1.0,-1.0,1.0;-1.0,1.0,1.0,-1.0];
  end
%
%  Compose N and L to get the fourth order operator.
%
  A = N * L;
%
%  Use boundary conditions to eliminate ghost points.
%
  A([gl;gr;gt;gb],:) = 0.0;
%
%  Left.
%
  for i = gl'
    A(i,[i+n,i,i+n-1,i+n+1,i+2*n]) = [ 2*(1+mu), -1, -mu, -mu, -1 ];
  end
%
%  Right.
%
  for i = gr'
    A(i,[i-n,i,i-n-1,i-n+1,i-2*n]) = [ 2*(1+mu), -1, -mu, -mu, -1 ];
  end
%
%  Top.
%
  for i = gt'
    A(i,[i+1,i,i+1+n,i+1-n,i+2]) = [ 2*(1+mu), -1, -mu, -mu, -1 ];
  end
%
%  Bottom.
%
  for i = gb'
    A(i,[i-1,i,i-1+n,i-1-n,i-2]) = [ 2*(1+mu), -1, -mu, -mu, -1 ];
  end
%
%  Eliminate the ghost points.
%
  phys = G(3:n,3:n);
  phys = phys(:);
  ghost = [ gl; gr; gt; gb ];
%
%  This is an obscure MATLAB operation that deserves a comment!
%
  A0 = A(phys,phys) - A(phys,ghost) / A(ghost,ghost ) * A(ghost,phys);
%
%  For right hand side, take into account half and quarter cells.
%
  B = speye ( n^2 );
  B(bl,bl) = B(bl,bl) / 2.0;
  B(br,br) = B(br,br) / 2.0;
  B(bt,bt) = B(bt,bt) / 2.0;
  B(bb,bb) = B(bb,bb) / 2.0;
  B0 = B(phys,phys);
%
%  Solve the generalized eigenvalue problem.
%
  [ V, Lambda ] = eigs ( A0 / h^4, B0, k, 'SM' );
  [ y, p ] = sort ( diag ( Lambda ) );
  x = [ -1:2/(n-3):1];
%
%  Plot each eigenmode.
%
  for i = 1 : k
    contourf ( x, x, reshape ( real ( V(:,p(i)) ), n-2, n-2 ), [0,0], 'k-' );
    axis ( 'equal' );
    s = sprintf ( '\\lambda_{%d}', i );
    title ( s );
    filename = sprintf ( 'chladni_figure_%02d.png', i );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  end

  return
end

