################################
# Discrete-time Logistic model #
################################
# Start with classic model: N(t+1)=N(t) + R*N(t)
# and ... N(t+1)-N(t)=R*N(t)
# Modify this to allow R to vary with N(t)
# R=R0*(1-N(t)/K) ... where R0=intrinsic growth rate
# when N(t) <<<< K, R ~= R0
# when N(t) = K, R=0
R0= 0.3; N0=2; K=150; t=45
N=matrix(0,t,1)
N[1]=N0
R=matrix(0,(t-1),1)
deltaN=matrix(0,(t-1),1)
for (i in 1:(t-1)){
R[i]=R0*(1-N[i]/K)
N[i+1]=N[i]+R[i]*N[i]
deltaN[i]=N[i+1]-N[i]}
plot(N[1:(t-1)],deltaN,xlab="N",ylab="N[t+1]-N[t]",
ylim=c(min(deltaN),max(R0*K/4)),
typ='l')
abline(h=0,col='gray')
abline(h=R0*K/4,col='gray')
abline(v=K/2,col='gray')
abline(h=8,col='gray')
# Now find the 2 population sizes for which a harvest of 8
# intersects the density-dependent annual recruitment curve
locator(2) # Click on the 2 points where the horizontal line
# at y=8 intersects the recruitment curve
# after you click twice, you'll see the x-y coordinates
# for the 2 points you clicked on. We're interested in
# the x-coordinate values
##########################################################
# QUOTA HARVEST with different starting population sizes #
##########################################################
############################
# Scenario A = START AT K #
############################
t=100
H=8
N=matrix(0,t,1)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# fill in appropriate value for starting value of N on next line #
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
N[1]=
R=matrix(0,(t-1),1)
deltaN=matrix(0,(t-1),1)
for (i in 1:(t-1)){
R[i]=R0*(1-N[i]/K)
N[i+1]=N[i]+R[i]*N[i]-H
deltaN[i]=N[i+1]-N[i]}
windows()
par(mfrow=c(3,1))
plot(1:t,N,xlab="time",type='o',ylim=c(0,(K+10)))
abline(h=c(K,K/2,0),col='gray')
abline(h=115,col='gray',lty=2)
####################################
# Scenario B = START JUST ABOVE X1 #
####################################
N=matrix(0,t,1)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# fill in appropriate value for starting value of N on next line #
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
N[1]=
R=matrix(0,(t-1),1)
deltaN=matrix(0,(t-1),1)
for (i in 1:(t-1)){
R[i]=R0*(1-N[i]/K)
N[i+1]=N[i]+R[i]*N[i]-H
deltaN[i]=N[i+1]-N[i]}
plot(1:t,N,xlab="time",type='o',ylim=c(0,(K+10)))
abline(h=c(K,K/2,0),col='gray')
abline(h=115,col='gray',lty=2)
####################################
# Scenario B = START JUST BELOW X1 #
####################################
N=matrix(0,t,1)
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# fill in appropriate value for starting value of N on next line #
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
N[1]=
R=matrix(0,(t-1),1)
deltaN=matrix(0,(t-1),1)
for (i in 1:(t-1)){
R[i]=R0*(1-N[i]/K)
N[i+1]=N[i]+R[i]*N[i]-H
deltaN[i]=N[i+1]-N[i]}
plot(1:t,N,xlab="time",type='o',ylim=c(0,(K+10)))
abline(h=c(K,K/2,0),col='gray')
abline(h=115,col='gray',lty=2)