#############################################################################################
##Script to modelize the population dynamic of Tibetan antelope considering an anthropogenic
##Allee effect
#############################################################################################
##
##Population dynamics modelizes with a fixed growth rate (0.08) and a variable hunting effort (40-
##80).
##

##List of variable
#Ninit : Initial population size
#K : Carrying capacity
#r : Growth rate
#Einit : Hunting effort
#Q : Catchability coefficient
#Ccinit : Cost per unit catch
#Ce : Cost per unit effort
#Nr : Threshold rarity of Tibetan antelope
#Pinit : Price per unit catch
#p1 : Price
#w : Constant
#b : Constant
#alphainit : Measure of how rapidly poachers respond to changes in profit
#atc : Number of antelopes captures
#temps : Time


E_estim <- seq(40,80,1) #generate 37 values between 40 and 76

dynpop_NAAE<-matrix(nrow=101,ncol=41,dimnames=list(1:101,1:41))

for (j in 1:41){
calcule.par<-function(nbmax=101, Ninit=1000000, K=1000000, r=0.08, Einit=E_estim[j], Q=0.002,
                      Ccinit=213, Ce=426000, Nr=100000, Pinit=818.86161, p1=250, w=568.862, b=1.2,
                      alphainit=2.449e-11, atc=385, temps=7){

N=numeric(nbmax)
E=numeric(nbmax)
Cost_catch=numeric(nbmax)
P=numeric(nbmax)
alpha=numeric(nbmax)
Year=numeric(nbmax)

N[1]=Ninit
E[1]=Einit
Cost_catch[1]=Ccinit
P[1]=Pinit
alpha[1]=alphainit
Year=1950:2050

for(i in 2:nbmax){
# Inclusion of Poacher pressure
N[i]<-N[i-1]+(r*N[i-1]*(1-N[i-1]/K))-Q*E[i-1]*N[i-1]
E[i]<-E[i-1]+alpha[i-1]*(P[i-1]*Q*E[i-1]*N[i-1]-Ce*E[i-1])

#Cost component
Cost_catch[i]<- Ce/(N[i]*Q)

#Price component
P[i]<-p1 + (w*(K/N[i])^b)

#Hunting response component
alpha[i]<-E[i]/(atc*temps*(P[i]*N[i]-Cost_catch[i]*N[i]))

}

return(list(N=N, E=E ,Cost_catch=Cost_catch, P=P, alpha=alpha, Year=Year))
}

dynpop <- calcule.par()
dynpop_NAAE[,j]<- dynpop$N
}


library("matrixStats")
mean_dynpop_NAAE<-rowMeans(dynpop_NAAE)
sd_dynpop_NAAE<-rowSds(dynpop_NAAE)
Year<-1950:2050
dynpop_final<-cbind(Year,mean_dynpop_NAAE,sd_dynpop_NAAE)