function U = fhs2d(F,k,Dx,Dy)
%FHS2D  Fast 2D Helmholtz equation solver.
%  U = FHS2D(F,K,DX,DY) solves the Helmholtz equation in two-space 
%  dimensions with source term F(X,Y) defined on a rectangular domain 
%  assuming Dirichlet boundary conditions and for constant wavenumber K.  
%  The boundary terms are specified in the perimeter of the NY x NX matrix 
%  F generated using meshgrid. Thus EAST, WEST, SOUTH, and NORTH boundary 
%  conditions are to be specified, respectively, in F(:,1) F(:,END), F(1,:),
%  and F(END,:). Meshgrid widths in X and Y directions are specified in DX 
%  and DY, respectively. Inversion of the discrete five-point Laplacian is
%  perfomed using the fast sine transform. For best performance speed 
%  NX - 1 and NY - 1 should be chosen as powers of 2.
%
%  See also HS2D, PS2D, FPS2D.

%  Francisco J. Beron-Vera, 2006/11/15
%  $Revision: 1.1 $  $Date: 2006/11/22 21:11:20 $

% Get dimensions.
Nx = size(F,2)-2;
Ny = size(F,1)-2;

% Fold boundary into source term.
U = F;
U(2:end-1,2:end-1) = 0;
U(2:end-1,2:end-1) = 0;
F = F - 4*del2(U,Dx,Dy);
F([1 end],:) = [];
F(:,[1 end]) = [];

% Eigenvalues of discrete Laplacian (nab2).
[nx ny] = meshgrid(1:Nx,1:Ny);
L = 2/Dx^2*(cos(nx*pi/(Nx+1))-1) + 2/Dy^2*(cos(ny*pi/(Ny+1))-1);

% Solve nab2*U + k^2*U = F using fast sine transform.
U(2:end-1,2:end-1) = idst(idst(dst(dst(F)')'./(L+k^2))')';