function x = S(A,i,j,t,index)
% S = S(A,i,j,t,index)
%
% This function calculates the quantity S(t), defined in the section 
% "decomposing near-term sensitivity analysis" of the paper by J. Yearsley, 
% "Transient population dynamics and short-term sensitivity analysis of matrix
% population models", Ecological Modelling.
% S(t) can be defined as the contribution to the near-term sensitivity which requires no 
% additional (compared to an asymptotic sensitivity analysis) data collection. 
% S(t) describes the near-term sensitivity when a population is initially at a stable
% stage distribution, and it can be calculated from the Leslie matrix alone.
%
% Input arguments:
%   A is the Leslie matrix
%   i,j are the indices describing the element of the Leslie matrix for
%       which a sensitivity is being calculated.
%   t is the the time after an initial perturbation for which the sensitivity is
%       being calculated
%
% Optional argument:
%   index defines whether a component of S is being computed. S can be
%   decomposed into components from each eigenvector.
%     If index=0 then calculate S(t) in full
%     If index=i (i>0) then calculate ith term corresponding to the
%       contribution from the ith eigenvector
%
% Output:
%   The value for S(t), which is a vector with the same dimensions as the
%   population structure vector.

% Jon Yearsley 22/7/2003
% j.yearsley@macaulay.ac.uk

if nargin<5,
    index=0;
end

% Right eigenvectors and eigenvalues of Leslie matrix A
[w,lambda] = eig(A);
[~,ind]=sort(real(diag(lambda)));
lambda = diag(lambda);
lambda = lambda(ind);
% Sort the right eigenvectors into order
% Largest eigenvalue is lambda(N), smallest is lambda(1)
w = w(:,ind);
% Calculate left eigenvectors
v = inv(w)';
N = length(lambda);


if index==0 || index>N,
    % Add up all components
    list = 1:N;
else
    % Look at only one component
    list = index;
end

% Calculate S
x = 0;
for k=list,
    if k==N,
        x = x + lambda(N)^t * (conj(v(i,k))*w(j,k) * t*w(:,N)/lambda(N) + ...
            sum(w(:,1:N-1).*(ones(N,1)*(conj(v(i,1:N-1))*w(j,N)./(lambda(N)-lambda(1:N-1).'))),2));
    else
        x = x - lambda(k)^t * conj(v(i,k))*w(j,N)' * w(:,k)/(lambda(N)-lambda(k));
    end
end
