function [Phi,Q,Psi,Qx,Qy]=fdm2(x,y,Kx,Ky,FH,FQ,varargin)
% [Phi,Q,Psi,Qx,Qy]=fdm2d(x,y,Kx,[Ky],FH,FQ [,radial])
% 2D block-centred steady-state finite difference model
% x,y mesh coordinates, Kx,Ky conductivities, if Ky==[], Ky=Kx
% FH=fixed heads (NaN for ordinary points), Q=fixed nodal flows
% Phi,Q computed heads and cell balances
% Qx(Ny-1,Nx-2) horizontal cell face flow positive in positive x direction
% Qy(Ny-2,Nx-1) vertial    cell face flow, postive in positive y direction
% Psi(Ny,Nx-2)  stream function
% TO 991017  TO 000530 001026 070414 090314


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

if isempty(Ky), Ky=Kx; end

x=x(:)'; Nx=length(x)-1; dx=diff(x); xm=0.5*(x(1:end-1)+x(2:end));
y=y(:);  Ny=length(y)-1; dy=abs(diff(y));

Nodes = reshape(1:Nx*Ny,Ny,Nx);               % Node numbering
IE=Nodes(:,2:end);   IW=Nodes(:,1:end-1);
IS=Nodes(2:end,:);   IN=Nodes(1:end-1,:);

% resistances and conducctances
switch nargin
    case 7
        fprintf('Fdm2 in radial mode.\n')
        RX=(1./dy)*log(x(2:end-1)./xm(1:end-1))./(2*pi*Kx(:,1:end-1))+...
           (1./dy)*log(xm(2:end)./x(2:end-1)) ./(2*pi*Kx(:,2:end));
        RY=0.5/pi*dy*(1./(x(2:end).^2-x(1:end-1).^2))./Ky;
        Cx=1./RX; Cy=1./(RY(1:end-1,:)+RY(2:end,:));
    otherwise
        fprintf('Fdm2 in flat mode.\n')
        RX=0.5*(1./dy)*dx./Kx; Cx=1./(RX(:,1:end-1)+RX(:,2:end));
        RY=0.5*dy*(1./dx)./Ky; Cy=1./(RY(1:end-1,:)+RY(2:end,:));
end

A=sparse([IE(:); IW(:); IN(:); IS(:)],...
         [IW(:); IE(:); IS(:); IN(:)],...
        -[Cx(:); Cx(:); Cy(:); Cy(:)],...
         Ny*Nx,Ny*Nx,5*Ny*Nx);   % System matrix
Adiag= -sum(A,2);                % Main diagonal

IAct=(Kx>0 | Ky>0);
I   = find(IAct &  isnan(FH));
Ifh = find(IAct & ~isnan(FH));
IAct= find(IAct);

Phi=FH;
Q  =zeros(Ny,Nx);

Phi(I)=spdiags(Adiag(I),0,A(I,I))\(FQ(I)-A(I,Ifh)*FH(Ifh)); % solve

Q(IAct)=spdiags(Adiag(IAct),0,A(IAct,IAct))*Phi(IAct);		% reshape

Qx=-Cx.*diff(Phi,1,2)*sign(x(end)-x(1));   % Flow across vertical   cell faces
Qy=-Cy.*diff(Phi,1,1)*sign(y(end)-y(1));   % Flow across horizontal cell faces

Psi= [flipud(cumsum(flipud(Qx)));zeros(size(Qx(1,:)))];
