function u_val = poisson(pos,rectangle,f,bdry,solver,h);

% u_val = poisson(pos,rectangle,f,bdry,h,solver)
%
% solves poisson equation with dirichlet boundary conditions
% on rectangle [-a,a] x [-b,b] using a five-point stencil
%
%   pos         [x,y] a point of the rectangle
%   rectangle   [a,b]
%               2a   length in x-direction
%               2b   length in y-direction
%   f           right-hand side of -\Delta u(x) = f   (constant)
%   bdry        Dirichlet data [x0,xa,y0,yb]
%               x0  boundary value on {-a} x [-b,b]   (constant)
%               xa  boundary value on { a} x [-b,b]   (constant)
%               y0  boundary value on [-a,a] x {-b}   (constant)
%               yb  boundary value on [-a,a] x { b}   (constant)
%   solver      'Cholesky'  sparse Cholesky solver
%               'FFT'       FFT based fast solver
%   h           h discretization parameter
%
%   u_val       solution at pos = [x,y]

% the grid
n = ceil(1/h);                              
a = rectangle(1); b = rectangle(2);
nx = 2*ceil(a)*n-1; ny = 2*ceil(b)*n-1; x=1:nx; y=1:ny; 
dx = (2*a)/(nx+1); dy = (2*b)/(ny+1);

% the right hand side
r = f*ones(nx,ny); 
r(1,y) = r(1,y)+bdry(1)/dx^2; r(nx,y) = r(nx,y)+bdry(2)/dx^2;
r(x,1) = r(x,1)+bdry(3)/dy^2; r(x,ny) = r(x,ny)+bdry(4)/dy^2;

% solve it
switch solver
    case 'Cholesky' % [Dem87,§6.3.3]
        Ax = spdiags(ones(nx,1)*[-1 2 -1]/dx^2,-1:1,nx,nx);
        Ay = spdiags(ones(ny,1)*[-1 2 -1]/dy^2,-1:1,ny,ny);
        A  = kron(Ay,speye(nx)) + kron(speye(ny),Ax);
        u = A\r(:); u = reshape(u,nx,ny);
    case 'FFT'      % [Dem87,§6.7]
        u = dst2(r);    
        d = 4*(sin(x'/2*pi/(nx+1)).^2*ones(1,ny)/dx^2+ ...
                        ones(nx,1)*sin(y/2*pi/(ny+1)).^2/dy^2);
        u = u./d;       
        u = dst2(u);   
    otherwise
        error('MATLAB:badopt','%s: no such solver known',solver);
end
        
% the solution
u_val = u(round((pos(1)+a)/dx),round((pos(2)+b)/dy));

return

% subroutines for 1D and 2D fast sine transform [Dem97,p.324]

function y = dst(x) 
n = size(x,1); m = size(x,2);
y = [zeros(1,m);x]; y = imag(fft(y,2*n+2)); 
y = sqrt(2/(n+1))*y(2:n+1,:);
return

function y = dst2(x)
y = dst(dst(x)')';
return