function drawdown=hantushn(Q,r,t,St,c,Sf,kD,N)
% drawdown=hantushn(Q,r,t,St,c,Sf,kD[,N])
% Analytical N-layer solution of groundwater flow in semi-confined multiple aquifer system
% Back-transformation from Laplace space according to Stehfest
% Solution by Jochem Ftirz and Rob Tijsen (2008) probably taken from another reference
% used by the interpretation of their pumping test De Zilk, Pastoorslaan Sept 3-14 2007
% Q = extraction vector (positive extractions will create positive dd')
% r = distance vector
% t = time vector
% drawdown will be in [length(Aquif,length(r),length(t)] format
% If instead r=r_iz_t list, the drawdown will be in list format
% provinding the drawdown for every coordinate pair
% in that case t is dummy !
% St = Storage coefficient of aquitards (first and last layers are aquitards)
% c = aquitard resistance vector one more than aquifers
% Sf = Storage coefficient of aquifers (second layer in a aquifer)
% kD = Transmissivity vector
% if length(c)=length(kD)+1, then also length(St)=length(Sf)+1
% in that case the lowest layer is an aquitard below which the head
% of zero is maintained, which may be tricky. If however
% length(kD)==length(c), then the lowest layer is an aquifer with
% no more layers below, and so the bottom is impvervious
% N  = Stehfest's parameter, default 10
% TO 090329, 0390330

fprintf('Running Hantushn\n');
%% Selftest ??
if nargin<7
    drawdown=selftest;
    return;
elseif nargin<8
    N=10;
end  % Stehfest's parameter

%% Check if r and t together form a list, in this case they mustb be column
% vectors of equal length
islist=(size(r,2)==3); % interpreting r as [r iz t] triples convenient for pumping tests
                     % further values are not used
%% Housekeepiong
Q=Q(:); t=t(:); St=St(:); c=c(:); Sf=Sf(:); kD=kD(:);
if islist,r_iz_t=r; else r=r(:); end

%% Stehfest v-coefficients, need to be calculated only once
v=NaN(N,1);
for i = 1:N
    dum=0;
    for k=floor((i+1)/2):min(i,N/2)
        dum=dum+k^ (N/2)*factorial(2*k)/(factorial(N/2-k)*factorial(k)*factorial(k-1)*factorial(i-k)*factorial(2*k-i));
    end
    v(i)=(-1)^(i+N/2)*dum;
end

%% Compute drawdown
switch islist
    case 1  % [r z] is just a list of coordinate pairs
        drawdown=zeros(size(r_iz_t,1),1);  % length(r)==length(t)
        for i=1:length(r_iz_t)
          dd=ddOnePoint(Q,r_iz_t(i,1),r_iz_t(i,3),St,c,Sf,kD,v);
          drawdown(i)=dd(r_iz_t(i,2));
        end
    case 0 % r z define an array
        drawdown=zeros(length(kD),length(r),length(t));
        for ir=1:length(r)
            for it=1:length(t)
                drawdown(:,ir,it)=ddOnePoint(Q,r(ir),t(it),St,c,Sf,kD,v);
            end
        end
end

function dd=ddOnePoint(Q,r,t,St,c,Sf,kD,v)
s=zeros(size(kD));
for iStehfest=1:length(v)
    p=iStehfest*log(2)/t;
    d=p*Sf./kD;
    b=sqrt(p*St.*c);
    if length(kD)==1
        eii=(b(1).*coth(b(1)))./(c(1).*kD);
        eij=(b(2).*coth(b(2)))./(c(2).*kD);
        A=eii+eij+d;
    else
        if length(c)>length(kD)
            bcothb=b.*coth(b); bcothb(isnan(bcothb))=1; % zero aquitard storage
            bsinhb=b./sinh(b); bsinhb(isnan(bsinhb))=1; % zero aquitard storage
            eii= bcothb(1:end-1)./(c(1:end-1).*kD);
            eij= bcothb(2:end)  ./(c(2:end)  .*kD);
            fii=bsinhb(2:end-1)./(c(2:end-1).*kD(2:end));
            fij=bsinhb(2:end-1)./(c(2:end-1).*kD(1:end-1));
            A=diag(eii+eij+d)-diag(fii,-1)-diag(fij,+1);
        else
            bcothb=ones(size(b)); %% b.*coth(b); bcothb(isnan(bcothb))=1; % zero aquitard storage
            bsinhb=ones(size(b)); %% b./sinh(b); bsinhb(isnan(bsinhb))=1; % zero aquitard storage
            eii= bcothb./(c.*kD);
            eij= [bcothb(2:end)./(c(2:end).*kD(1:end-1));0];
            fii=bsinhb(2:end)./(c(2:end).*kD(2:end));
            fij=bsinhb(2:end)./(c(2:end).*kD(1:end-1));
            A=diag(eii+eij+d)-diag(fii,-1)-diag(fij,+1);
        end
    end
    [V,D]=eig(sqrtm(A));
    s=s+(v(iStehfest)/(2*pi*p))*V*diag(besselk(0,r*diag(D)))*V^(-1)*(Q./kD);
end
dd=s*log(2)/t;

function dd=selftest

c =[1000 1500 1000 4000 2000];
St=[   3     5   3    2 1]*1e-3;
kD=[2000  1500 500 2000];
Sf=[  10     4   1    3]*1e-4;
Q =[0, 10000, 0, 0];

r =logspace(1,log10(6000),40);
t =[1e-3 1e-2 1e-1 1 10];

dd=hantushn(-Q,r,t,St,c,Sf,kD);

for it=1:length(t)
    semilogx(r,dd(:,:,it)); hold on
end