################################
# 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)