function [Phi,Qt,Qx,Qy,Qs]=fdm2t(x,y,t,kx,ky,S,IH,FH,FQ,varargin)
% [Phi,Q,Qx,Qy,Qs]=fdm2d(x,y,t,kx,ky|[],S,IH,FH,FQ [,radial])
% 2D block-centred transient finite difference model
% INPUT:
%  x(Nx+1)         x-coordinate of mesh/grid
%  y(Ny+1)         y-coordinate of mesh/grid
%  t(Nt+1)         time for output. 0 will be added if necessary
%  kx(Ny,Nx)       conductivity in x-direction
%  ky(Ny,Nx)       same in y direction, ky=kx if ky==[]
%  S(Ny,Nx)        primary storage (S+Sy)
%  IH(Ny,Nx)       initial head
%  FH(Ny,Nx)       fixed heads (NaN for ordinary points), Q=fixed nodal flows
%  Radial          Arbitrary input caused model to assume axial flow with
%                 x=r and FQ are ring flows while als flows have dimension
%                 L3/T instead of L2/t. For clearness use 'R' or 'Radial'
%                 as input at the function call.
% OUTPUT
%  Phi(Ny,Nx,Nt+1) computed heads with Phi(Ny,Nx,1) initial heads for t=0
%, Qt(Ny,Nx,Nt)    computed total cell balance during timestep it
%  Qx(Ny,Nx-1,Nt)  hor.  cell face flow in x-direction positive along increasing col number
%  Qy(Ny-1,Nx,Nt)  vert. cell face flow in y-directin  positive along increasing row number
%  Qs(Ny,Nx,Nt)    storage change of node during timestep it
%
% TO 991017  TO 000530 001026 070414 080301


% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later
% TO 080226 (added inactive cells and true fixed heads, t(1) is now always made zero

% Copyright 2009 Theo Olsthoorn, without any warranty
% under free software foundnation GNU version 3 or later licence

theta=0.67;  % implicitness

if isempty(ky), ky=kx; end

if t(1)~=0, t=unique([0 abs(t)]); else t=unique(abs(t)); end     % Initial time is always zero

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

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
if nargin>9
    fprintf('Fdm2t 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,:));
    Cs=pi*dy*(x(2:end).^2-x(1:end-1).^2).*S;
else
    fprintf('Fdm2t 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,:));
    Cs=dy*dx.*S;
end
Cs=Cs(:);  % storage conductacne when devided by dt*theta

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 =Nodes( kx>0 | ky>0 | S>0);               % active cells
I    =Nodes((kx>0 | ky>0 | S>0) &  isnan(FH)); % active cells but not fixed heads = cells with heads to be computed
Ifh  =Nodes((kx>0 | ky>0 | S>0) & ~isnan(FH)); % active cells with fixed heads

Phi=NaN(Ny*Nx,Nt+1);  % allocate space to store the entire head matrix
Qt =NaN(Ny*Nx,Nt);    % allocate memory for Qt
Qs =NaN(Ny*Nx,Nt);    % allocate memory for Qs
Qx =NaN(Ny,Nx-1,Nt);  % allocate memory for Qx
Qy =NaN(Ny-1,Nx,Nt);  % allocate memory for Qy


Phi(IAct,1)=IH(IAct);  % store initial head at Phi(:,:,1)
FQ=FQ(:); 
Fi =FH(:);            % head computed in this time step
for it=1:Nt
    Fi (I)         =spdiags(Adiag(I)+Cs(I)/dt(it)/theta,0,A(I,I))\(FQ(I)-A(I,Ifh)*FH(Ifh)+Cs(I).*Phi(I,it)/dt(it)/theta); % solve
    Phi(IAct ,it+1)=Fi(IAct)/theta-(1-theta)/theta*Phi(IAct,it);
    Qt (IAct ,it  )=spdiags(Adiag(IAct),0,A(IAct,IAct))*Fi(IAct);
    Qs (IAct ,it  )=-Cs(IAct).*(Phi(IAct,it+1)-Phi(IAct,it))/dt(it);   % Storage in time step m3 for cell
    Qx (:,: ,it   )=-Cx.*diff(reshape(Fi,Ny,Nx),1,2);	% Flow across horizontal cell faces m3/d per cell
    Qy (:,: ,it   )=-Cy.*diff(reshape(Fi,Ny,Nx),1,1);   % Flow across vertical cell faces, m3/d per cell
end
Phi=reshape(Phi,Ny,Nx,Nt+1);                   % NaN if inactive
Qt =reshape(Qt ,Ny,Nx,Nt  ); Qt(isnan(Qt))=0;  % 0 if inactive
Qs =reshape(Qs ,Ny,Nx,Nt  ); Qs(isnan(Qs))=0;  % 0 if inactive
Qx(isnan(Qx))=0; Qy(isnan(Qy))=0;              % 0 if inactive
