function sense=near_term(A,i,j,t,n0,output)
% sense = near_term(A,i,j,t,n0,output)
%
% A function that calculates the near-term sensitivity of a
% matrix model a time t after a perturbation, given the initial conditions.
% Uses the approach described in Yearsley (2004) 
% "Transient population dynamics and  short-term sensitivity analysis 
% of matrix population models", Ecological Modelling, Vol 177, p245-258

% There is no numerical difference between this and near_term_FG.m
% near_term_FG.m uses the results of Fox and Gurevitch "Population numbers count:
% Tools for near-term demographic analysis.", Am. Nat., 156, p242--256
%
% Input arguments:
%   A is the Leslie matrix
%   i,j are the indices describing the element of the Leslie matrix which is
%     being perturbed
%   t is the the time after the perturbation for which the sensitivity is
%     being calculated
%   n0 is the initial population structure
%   output determines the kind of sesnitivity being calculated
%     output = 0  sensitivity of average growth rate for each stage
%     output = 1  sensitivity of average growth rate
%     output = 2  sensitivity of the population size [default]
%
% Output:
%   The sensitivity of the population size, or the growth rate of the 

% Jon Yearsley 01/6/05
% j.yearsley@macaulay.ac.uk

if nargin<5,
    error('!! Incorrect number of input arguments');
end
if nargin==5,
    output = 2;
end
if output>2,
    disp('!! Assuming required output is sensitivity of population size');
    disp('   If this is not your intention, reset the argument "output"');
    output = 2;
end

% resize n0
n0 = reshape(n0,length(n0),1);

% Calculate sensitivity of population size of each class
dnda = delta_S(A,i,j,t) * n0;

if output==0,
    % Calculate sensitivity of population growth rate for each stage
    dnda_t2 = delta_S(A,i,j,t+1) * n0;
    nt = A^t * n0;
    nt2 = A*nt;
    dlambdada_t = (dnda_t2 - dnda .* nt2./nt)./ nt;
    sense = dlambdada_t;
elseif output==1,
    % Calculate sensitivity of population growth rate (sum over all stages)
    dnda_t2 = sum(delta_S(A,i,j,t+1) * n0);
    nt = A^t * n0;
    nt2 = A*nt;
    dlambdada_t = (dnda_t2 - sum(dnda) .* sum(nt2)./sum(nt))./ sum(nt);
    sense = dlambdada_t;
else
    sense = dnda;
end
if any(imag(sense)>1e-10),
    warning('Sensitivities are complex! Possible error.')
end

% Sensitivities should not be complex (although numerical errors may
% introdcue a small complex component)
sense = real(sense);
