% PROGRAM tbar_ceiling
%    Plots mean time to extinction vs. K and starting population 
%    size for the ceiling model, using expressions for the mean
%    time to extinction from Lande, Am. Nat. 142:911-927 (1993),
%    Foley, Conservation Biology 1:124-137 (1994), and 
%    Middleton et al., Theor. Pop. Biol. 48:277-305 (1995).
% 
%   Note: Eq. 4.3 and 4.4 and Box 4.1 in Morris & Doak 2002 
%   assume that mu is not zero. This program has been modified 
%   to account for the possibility that mu might be zero, using 
%   formulas from Foley 1994
%
%   Updated 1/24/07

mu=0.1;				% mean log population growth rate
sigma2=0.1;			% variance of log population growth rate	
c=mu/sigma2;

% Plot mean time to extinction (N=1) vs. K for populations 
% starting at K, according to Equation 4.4
K=1:50;		% vector of carrying capacities
if mu~=0
    Tbar=(K.^(2*c)-1-2*c*log(K))/(2*mu*c); % Eq. 4.4
else
    Tbar=log(K).^2/sigma2;  % See eq. 7 in Foley 1994
end
subplot(2,1,1);     % as first of 2 stacked figures on page,
plot(K,Tbar)        % plot Tbar vs K
xlabel('Carrying capacity, K','FontSize',12)
ylabel('Mean time to extinction, Tbar','FontSize',12)
title('Populations starting at K','FontSize',14)

% Plot mean time to extinction (N=1) vs. Nc for populations
% starting at various Nc<=K, according to Equation 4.3
K=50;				% carrying capacity
k=log(K);           % log of carrying capacity
Nc=1:K;             % vector of initial population sizes
d=log(Nc);          % log distance to extinction threshold
if mu~=0
    Tbar=(exp(2*c*k)*(1-exp(-2*c*d))-2*c*d)/(2*mu*c); % Eq. 4.3
else
    Tbar=(2*d/sigma2).*(k-d/2); % see eq. 6 in Foley 1994
end
subplot(2,1,2);     % as second of 2 stacked figures on page,
plot(Nc,Tbar)       % plot Tbar vs Nc
xlabel('Initial population size, Nc','FontSize',12)
ylabel('Mean time to extinction, Tbar','FontSize',12)
title('Populations starting below K','FontSize',14)