function w = helmsolve(kappa,f)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function w = helmsolve(kappa,f)
% This function computes an approximate solution to
% the Helmholtz equation
% 
%     - u_{xx} - u_{yy} - kappa^2 u = g
% 
% on the unit square, with u=0 on the boundary of the
% square.  
% 
% The discretization is a 2nd order accurate finite
% difference method, with mesh points (x_j,y_k),
%   x_j = jh,  y_k = kh,  j,k=1,...,n, and h = 1/(n+1).
% 
% The solution algorithm is a "fast Poisson solver",
% using the discrete sine transform to diagonalize
% the matrix.  It will be unstable if kappa is too
% close to an eigenvalue of the discrete Laplacian.
%
% The program is designed to be rather efficient in
% time and storage (although the nxn array lam
% could be stored in place of f if storage is an issue).
% 
% On input:
%    kappa  is a scalar equal to the parameter in the
%           differential equation
%    f      is an n x n array, with f(j,k) = g(x_j,y_k).
% On ouput:
%    w      is an n x n array, 
%           with w(j,k) approx. u(x_j,y_k).
% 
% helmsolve.m     Dianne O'Leary  05/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,m] = size(f)

% Take the Fourier transform in the y direction.

w = dst(f);

% Take the Fourier transform in the x direction.

w = dst(w')';

% Divide by the eigenvalues.

j = [1:n]';
cosj = cos(j*pi/(n+1));
lam = (n+1)^2*(4 -2*(cosj*ones(1,n) + ones(n,1)*cosj')) - kappa^2;
w = w ./ lam;

% Take the inverse Fourier transform in the x direction.

w = idst(w')';

% Take the inverse Fourier transform in the y direction.

w = idst(w);