function [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu)
% PLATE_MATRICES  Sparse Laplacian, stiffness and mass matrices for the plate.
%
%   [L, K, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, bc, nu) assembles
%   the spatial operators of the finite-difference schemes as sparse matrices
%   acting on the vector w(idx) of unknowns, where idx indexes the ny x nx field
%   with y running fastest (column-major, as w(:) does) and the unknowns form an
%   my x mx array.  The semi-discrete equation is
%
%     M w_tt + 2 sigma0 M w_t - 2 sigma1 L w_t + kappa^2 K w = F(t) M s / mu
%
%   Simply supported and clamped edges (w = 0 on the edge nodes).  The unknowns
%   are the mx = nx-2 by my = ny-2 interior nodes, M is the identity, L the
%   five-point Laplacian and K the thirteen-point biharmonic, with the edges
%   handled exactly as in the ghost-node scheme of plate_fd / plate_pad: D2*D2
%   encodes the odd reflection (w_nn = 0), and for clamped edges the even
%   reflection changes the end coefficient of the fourth difference from 5/dx^4
%   to 7/dx^4.  The mixed term 2 w_xxyy never reaches a ghost node.  nu is not
%   used.  A test in the lab checks K against the explicit stencil.
%
%   Free edges.  Every node is an unknown (mx = nx, my = ny).  The free-edge
%   conditions, no bending moment and no Kirchhoff shear on an edge and no
%   twisting moment at a corner, involve third derivatives and do not fit in a
%   ghost layer.  K is therefore assembled from the bending energy
%
%     U = D/2 * integral( w_xx^2 + w_yy^2 + 2 nu w_xx w_yy + 2 (1-nu) w_xy^2 ) dA
%
%   using second differences at the nodes, the cross difference at the cell
%   centres and trapezoidal quadrature weights (the finite-difference energy
%   method; Bushnell, Computers & Structures 1973).  The free-edge conditions
%   are then natural boundary conditions of the discrete energy: the scheme is
%   second order inside and first order at the boundary.  M holds the trapezoidal
%   node weights (1 inside, 1/2 on an edge, 1/4 at a corner) and L is the
%   Laplacian energy operator with the natural zero-normal-derivative condition.
%   K is symmetric positive semidefinite and its null space is exactly the three
%   rigid-body motions 1, x and y, which plate_rigid_modes removes from the
%   state.  All three matrices are divided by dx*dy, so inside the plate M = 1
%   and K and L reduce to the usual stencils, and the largest eigenvalue of
%   M \ K stays below the interior bound (4/dx^2 + 4/dy^2)^2 that plate_fd_plan
%   uses for the time step.
%
%   Everything is built from Kronecker products so the matrices stay sparse;
%   kron(Ax, By) acts on w(:) as By * W * Ax'.

if nargin < 6
    nu = NaN;
end
switch bc
    case {'simply-supported', 'clamped'}
        mx = nx - 2;
        my = ny - 2;
        ex = ones(mx, 1);
        ey = ones(my, 1);
        D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
        D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
        Ix = speye(mx);
        Iy = speye(my);
        D4x = D2x * D2x;
        D4y = D2y * D2y;
        if strcmp(bc, 'clamped')
            D4x(1, 1) = 7 / dx^4;   D4x(mx, mx) = 7 / dx^4;
            D4y(1, 1) = 7 / dy^4;   D4y(my, my) = 7 / dy^4;
        end
        L = kron(Ix, D2y) + kron(D2x, Iy);
        K = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
        M = speye(mx * my);
        mask = false(ny, nx);
        mask(2:end - 1, 2:end - 1) = true;
        idx = find(mask);
    case 'free'
        if ~isfinite(nu)
            error('plate_matrices: free edges need Poisson''s ratio nu');
        end
        mx = nx;
        my = ny;
        % One-dimensional pieces, as triplets so the layout is unambiguous:
        % second differences at the interior nodes ((n-2) x n), first differences
        % at the cell midpoints ((n-1) x n), selectors of the interior nodes, and
        % trapezoidal node weights.
        kx = (1:nx - 2)'; ky = (1:ny - 2)';
        d2x = sparse([kx; kx; kx], [kx; kx + 1; kx + 2], [ones(nx - 2, 1); -2 * ones(nx - 2, 1); ones(nx - 2, 1)], nx - 2, nx) / dx^2;
        d2y = sparse([ky; ky; ky], [ky; ky + 1; ky + 2], [ones(ny - 2, 1); -2 * ones(ny - 2, 1); ones(ny - 2, 1)], ny - 2, ny) / dy^2;
        jx = (1:nx - 1)'; jy = (1:ny - 1)';
        d1x = sparse([jx; jx], [jx; jx + 1], [-ones(nx - 1, 1); ones(nx - 1, 1)], nx - 1, nx) / dx;
        d1y = sparse([jy; jy], [jy; jy + 1], [-ones(ny - 1, 1); ones(ny - 1, 1)], ny - 1, ny) / dy;
        Sx = sparse(kx, kx + 1, ones(nx - 2, 1), nx - 2, nx);
        Sy = sparse(ky, ky + 1, ones(ny - 2, 1), ny - 2, ny);
        tx = ones(nx, 1); tx([1 end]) = 0.5; Tx = spdiags(tx, 0, nx, nx);
        ty = ones(ny, 1); ty([1 end]) = 0.5; Ty = spdiags(ty, 0, ny, ny);
        % The four terms of the energy density, each as D' * (weights) * D
        Kxx = kron(d2x' * d2x, Ty);                                                % w_xx^2 at nodes 2..nx-1, every row
        Kyy = kron(Tx, d2y' * d2y);                                                % w_yy^2 at nodes 2..ny-1, every column
        Kc = nu * (kron(d2x' * Sx, Sy' * d2y) + kron(Sx' * d2x, d2y' * Sy));      % 2 nu w_xx w_yy where both exist
        Kxy = 2 * (1 - nu) * kron(d1x' * d1x, d1y' * d1y);                         % 2 (1-nu) w_xy^2 at the cell centres
        K = Kxx + Kyy + Kc + Kxy;
        L = -(kron(d1x' * d1x, Ty) + kron(Tx, d1y' * d1y));
        M = kron(Tx, Ty);
        idx = (1:nx * ny)';
    otherwise
        error('plate_matrices: unknown boundary condition "%s"', bc);
end
end