function f = membraneForces(z,nx,ny,T,dx)

% Forces on a 2D membrane of beads, given by 
%   f(i,j) = k*( z(i+1,j) + z(i-1,j) + z(i,j+1) + z(i,j-1) - 4*z(i,j) ),
% using periodic boundary conditions (z(nx+1,j)=z(1,j), z(i,ny+1)=z(i,1))
% Input:
%   z(1,nx*ny) : vertical displacements of nx*ny beads (m)
%   nx         : number of beads in x direction
%   ny         : number of beads in y direction
%   T      : string tension (N)
%   dx     : distance between beads (m)
% Output:
%   f(1,nx*ny) : forces on nx*ny beads (N)
% JMS, Nov.2019

z = reshape(z,nx,ny);
f = T/dx*( circshift(z,[+1,0]) + ...
           circshift(z,[-1,0]) + ...
           circshift(z,[0,+1]) + ...
           circshift(z,[0,-1]) - 4*z );
f = reshape(f,1,nx*ny);

end
