clear all
t=[0:65]';
N=[18 22 26 16 19 21 18 22 25 31 30 34 31 25 21 24 21 28 24 26 32 33 36 39 32 33 42 44 43 48 50 56 57 59 51 49 49 57 69 72 75 76 78 73 73 75 86 97 110 134 138 146 146 132 136 143 133 158 160 180 183 188 180 176 185 194]'
size(N)
X=[ones(size(t)) t];
p=X\log(N)
q=exp(X*p);
plot(t,N,'r*',t,q,'-og')
xlabel('years after 1938')
ylabel('whooping crane population')
legend('data','malthus fitted model')
N1938=exp(p(1))
N2004=exp(p(1)+66*p(2))
N2005=exp(p(1)+67*p(2))
N2006=exp(p(1)+68*p(2))

K=1000;
c=log(N./(K-N))
d=X\c
R=exp(d(2))
N0=(K*exp(d(1)))/(1+exp(d(1)))
Nf=(K*N0*exp(d(2)*t))./(K-N0+N0*exp(d(2)*t));

plot(t,N,'r*',t,Nf,'-ob')
xlabel('years after 1938')
ylabel('whooping crane population')
legend('data','Beverton-Holt fitted model')

%prediction for 2004, 2005, 2006

Ne2004=(K*N0*exp(d(2)*66))./(K-N0+N0*exp(d(2)*66))
Ne2005=(K*N0*exp(d(2)*67))./(K-N0+N0*exp(d(2)*67))
Ne2006=(K*N0*exp(d(2)*68))./(K-N0+N0*exp(d(2)*68))