function [Phi,Q,Qx,Qy,Qz]=fdm3(x,y,Z,kx,ky,kz,FH,FQ) %% future: ,hdr,Cdr,hrv,hbr,Crv,hgh,Cgh)
% [Phi,Q,Qx,Qy,Qz]=fdm3(x,y,Z,x,kx,ky,kz,FH,FQ)
% 3D block-centred steady-state finite difference model
% x,y,Z mesh coordinates
% Z may be a vector, a 3D vertical vector or a full size 3D array of size Ny,Nx,Nz+1
% kx,ky,kz conductivity arrays
% FH=fixed heads (NaN for ordinary points), Q=fixed nodal flows
% Phi,Q computed heads and cell balances
% Qx,Qy,Qz is cell face flow, positive along positive axis direction
% Future:  hdr=elevation drains, Cdr=resistance drains
%          hrv=elevation river , Crv=resistance river, hbr=elevation bottom river
%          hgh=elevation general head, Cgh is resistance general head
% TO 991017  TO 000530 001026 070414 070513
% TO 080226 implemented inactive cells and true fixed heads
% TO 090216 small edits


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

x=x(:)'; dx=abs(diff(x)); Nx=length(dx);
y=y(:);  dy=abs(diff(y)); Ny=length(dy);

% z may be a 3D array with the z values of the bottom and top of all cells
sz=size(Z);
if length(sz)<3 % than it is a vector which will be orientated vertially
    if ~(any(sz)==1),
        error('Z must be vector or a full size 3D matrix\n');
    end
    Z=permute(flipud(Z(:)),[1,1,length(Z(:))]); dz=abs(diff(Z)); Nz=length(dz);
    [Dx,Dy,Dz]=meshgrid(dx,dy,dz);
else  % z is specified as a 3D array, then ...
    Nz=sz(end)-1;
    if sz(1)==1 && sz(2)==1  % only a vertical column of z was given
        Z=Z(ones(Ny,1),ones(1,Nx),1:Nz+1); %overwrite into 3D array
    end
    Dz=abs(diff(Z,1,3));
    Dx=dx(ones(Ny,1),1:Nx,ones(Nz,1));
    Dy=dy(1:Ny,ones(1,Nx),ones(Nz,1));
end

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

% resistances and conducctances
warning('off','all');  % to allow zero conductivities for inactive cells
RX=0.5*Dx./(Dy.*Dz.*kx); Cx=1./(RX(:,1:end-1,:)+RX(:,2:end,:));
RY=0.5*Dy./(Dz.*Dx.*ky); Cy=1./(RY(1:end-1,:,:)+RY(2:end,:,:));
RZ=0.5*Dz./(Dx.*Dy.*kz); Cz=1./(RZ(:,:,1:end-1)+RZ(:,:,2:end));
warning('on','all');

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

IAct =Nodes( kx>0 | ky>0 | kz>0);               % active cells
I    =Nodes((kx>0 | ky>0 | kz>0) &  isnan(FH)); % active cells but not fixed heads = cells with heads to be computed
Ifh  =Nodes((kx>0 | ky>0 | kz>0) & ~isnan(FH)); % active cells with fixed heads

Phi=FH(:);
FQ =FQ(:);
Q  =FQ;

Phi( I)=spdiags(Adiag(I   ),0,A(I   ,I   ))\(FQ(I)-A(I,Ifh)*Phi(Ifh)); % solve
Q(IAct)=spdiags(Adiag(IAct),0,A(IAct,IAct))*Phi(IAct );		           % reshape
  
Phi=reshape(Phi,Ny,Nx,Nz);
Q  =reshape(Q  ,Ny,Nx,Nz);
Qx=-Cx.*diff(Phi,1,2)*sign(x(end)-x(1)); Qx(isnan(Qx))=0;  % Flow across horizontal cell faces
Qy=-Cy.*diff(Phi,1,1)*sign(y(end)-y(1)); Qy(isnan(Qy))=0;  % Flow across vertical cell faces
Qz=-Cz.*diff(Phi,1,3)*sign(Z(end)-Z(1)); Qz(isnan(Qz))=0;  % Flow across vertical cell faces
