function sol = simplePoisson(xGrid,yGrid,data)

% Copyright (C) 2017 Leon van Dommelen.  This software is
% made available under the GNU Public License version 3.
%
% Function that finds the solution to the two-dimensional
% Poisson equation on a rectangular grid when the
% interior "forcing" and solution on the boundaries is
% given.
%
% Input:
%
%   xGrid: an m x n array containing the values of the
%          first independent variable at the m x n grid
%          points.  All xGrid(:,j) values in the same
%          column j must be the same.  And the xGrid(:,j)
%          values must increase with increasing column
%          number j.
%
%   yGrid: an m x n array containing the values of the
%          second independent variable at the m x n grid
%          points.  All yGrid(i,:) values in the same row
%          i must be the same.  And the yGrid(i,:) values
%          must increase with increasing row number i.
%
%   data: an m x n array that contains data about the
%         desired solution.  The value of data(i,j) at
%         any interior point of the grid should be the
%         "forcing" at that point (i.e. whatever it is at
%         the point that wants to make the solution
%         nontrivial).  The value of data(i,j) at any
%         boundary point, (i=1, i=m, j=1, or j=n), should
%         be the desired value of the solution at that
%         point.
%
% Output:
%
%   sol: an m x n array containing the values of the
%        approximate solution at the m x n points.
%
% Hint: arrays xGrid and yGrid must be in the form that
% you would get from
%
%   [xGrid yGrid] = meshgrid(xValues,yValues)
%
% where xValues is a one-dimensional array of n
% increasing x-values, and yValues is a one-dimensional
% array of m increasing y-values.
%
% Warning: You would not want to use this function for
% very large values of m x n.  Its simple algorithm
% wastes lots of storage and computer time for big m x n.
% But at the time of writing, you can probably make m x n
% quite big, like say 100 x 100, without much trouble.

% find the size of the arrays
[m n] = size(xGrid)

% check that we have interior points
if m<3 | n<3
    fprintf('*** simplePoisson: ')
    error('There are no interior points!')
end

% check that all three arrays have the same size
[mm nn] = size(yGrid);
if mm~=m | nn~= n
    fprintf('*** simplePoisson: ')
    error('The input arrays are not all the same size!')
end
[mm nn] = size(data);
if mm~=m | nn~= n
    fprintf('*** simplePoisson: ')
    error('The input arrays are not all the same size!')
end

% one check that xGrid is constant in columns and yGrid in rows
if xGrid(1,2)~=xGrid(3,2) | yGrid(2,1)~=yGrid(2,3)
    fprintf('*** simplePoisson: ')
    error('The xGrid and/or yGrid inputs are incorrect!')
end

% compute the number of unknown solution values
N=m*n
% zero the matrix of the system of equations to solve
A=zeros(N);
% zero the vector of right hand sides
b=zeros(N,1);

% form the interior equations to solve
for i = 2:m-1
    for j = 2:n-1
        % compute the equation number
        I = i + (j-1)*m;
        % compute the mesh increments first
        dx1=xGrid(i,j)-xGrid(i,j-1);
        dx2=xGrid(i,j+1)-xGrid(i,j);
        dy1=yGrid(i,j)-yGrid(i-1,j);
        dy2=yGrid(i+1,j)-yGrid(i,j);
        if dx1<=0 | dx2<=0 | dy1<=0 | dy2<=0
            fprintf('*** simplePoisson: ')
            error('xGrid or yGrid values are not increasing!')
        end
        % find the average increments
        dxAve=.5*(dx1+dx2);
        dyAve=.5*(dy1+dy2);
        % diagonal coefficient of the equation
        A(I,I)=-1/(dx2*dxAve)-1/(dx1*dxAve)...
               -1/(dy2*dyAve)-1/(dy1*dyAve);
        % off-diagonal coefficients
        A(I,I+m)=1/(dx2*dxAve);
        A(I,I-m)=1/(dx1*dxAve);
        A(I,I+1)=1/(dy2*dyAve);
        A(I,I-1)=1/(dy1*dyAve);
        % right hand side of the equation
        b(I)=data(i,j);
    end
end

% add the x-boundary conditions at i = 1 and m
for j = 1:n
    % x-boundary condition at i = 1
    I = 1 + (j-1)*m;
    A(I,I)=1;
    b(I)=data(1,j);
    % x-boundary condition at i = m
    I = m + (j-1)*m;
    A(I,I)=1;
    b(I)=data(m,j);
end

% add the y-boundary conditions at j = 1 and n
for i = 2:m-1
    % y-boundary condition at j = 1
    I = i;
    A(I,I)=1;
    b(I)=data(i,1);
    % y-boundary condition at j = n
    I = i + (n-1)*m;
    A(I,I)=1;
    b(I)=data(i,n);
end

% check the condition number of the matrix
condA=cond(A);
relErrorMatlab=condA*eps(1)
if relErrorMatlab >= 0.1
    fprintf('*** simplePoisson: ')
    error('The equations are singular.  Check xGrid and yGrid!')
end

% find the solution as a vector
solVec=A\b;
% reformat it to the same form as the input arrays
sol=reshape(solVec,[m n]);

end
