%first run the code below to define the function York1994mod1()
%then you can run the code below to get scaled Leslie matrices as were used in the Holmes et al 2007 paper

%In the paper I rounded the scaling factors to 2 digits
pert.juv = .94;
pert.adult = 1.07;
pert.fec = .64;
pert.a = 1; pert.t = 1;  %nothing, leftover stuff in the York1994mod1() function
A = York1994mod1(pert);
max(eig(A))

%I would have gotten a slightly smaller max(eig) if I had used more digits
pert.juv = .9350;
pert.adult = 1.0676;
pert.fec = .6412;
pert.a = 1;
pert.t = 1;
A = York1994mod1(pert);
max(eig(A))


%to show the original matrix
pert.juv = 1;
pert.adult = 1;
pert.fec = 1;
pert.a = 1;
pert.t = 1;
A = York1994mod1(pert);  %now since all scaling factors are 1, this is the 1970s matrix
max(eig(A))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = York1994mod1(pert);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%York1994mod1() returns the corrected York 1994 lifehisotry matrix from Weibull
%Note that fecundity is already corrected for survivorship during pup period (but not surv from age 0 to 1)
%the survivorships are 0 to 1, 1 to 2, 2 to 3 etc
%the population vector starts with pups
%This is a modified Leslie matrix such that n_i including pups is number at the start of the breeding season
%First line is s_i f_i+1 as described on pg 40 in York 1994

%Estimated late-term pregnancy rates; this is like C&P table 26 perc mature * birth rate
%note that animal may be age 3 at implantation, but age 4 when in late-term pregnancy
%age
%group                  fit.mod1             fit.mod2
%(age at implantation)
%------------------------------------------------------
%3 (age 4 at late-term)   0.096                0.168
%4                        0.339                0.493
%5                        0.443                0.497
%6                        0.559                0.607
%7-9                      0.657                0.752
%10-16                    0.770                0.535
%16-20                    0.514                0.316
%21-30                    0.000                0.000
%------------------------------------------------------

if(nargin == 0) pert.fec = 1; pert.juv = 1; pert.a = 1; pert.t = 1; pert.adult = 1; end

%female birth rate at age i+1 is the late-term preg rate of animals impregnanted at age i
% get pregnant at age i and give birth (late-term preg) at age i+1
sexratio = 0.5;
sn = 0.949; %neonatal survival
f=sn*sexratio*[ 0 0  0  0.096 0.339 0.443 0.559 0.657*ones(1,3) 0.770*ones(1,6) 0.514*ones(1,5) zeros(1,11)];
f = pert.fec*f;
%note that this is f_i+1; i've removed the f_0 term so this is f_1 f_2 f_3 ....
%note that birth starts at age 4 (implantation at age 3)

%pups year i+1 is n_i+1 * f_i+1 = n_i * s_i * f_i+1
%note s_i is surv age i to age i+1
%x*s accomplishes this below

x=[ f
   
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
];

xx = 1:28;
a = 1.367*pert.a;
t = 8.416*pert.t;

%from York 1994
adultsurv = (gammainc(((xx + 1)./t).^a, 1/a) - gammainc( (xx./t).^a, 1/a ) ) ./ ...
   (gammainc( (xx./t).^a, 1/a ) - gammainc( ((xx-1)./t).^a, 1/a ) );

adultsurv = [adultsurv*pert.adult 0]; %die at age 32

juv1surv = 0.806* pert.juv;
juvsurv = [juv1surv juv1surv+(adultsurv(1)-juv1surv)/3 juv1surv+2*(adultsurv(1)-juv1surv)/3 ]; %constant set so that max(eig) = 1 per York 1994

s= diag([juvsurv adultsurv]);

x=x*s;  %modified Leslie so that pup in n_0 are pups at start of breeding season t
