#' A closed population function that calculates the length of time of the transient period after MPA implementation for a deterministic model
#'This function calculates the ratio change in a fished poulation after a marine protected area is implemented assuming closed population
#'method is from White et al. 2013 'transient responses of fished populations to marine reserve establishment published in conservation letters
#' assumes a lambda of 1.0 in the MPA
#' @param tf the number of time step to run the population
#' @param maxage max age of the species ie. number of age classes
#' @param Lmat length at maturity
#' @param M the natural mortality rate, if unknown generally use 0.2
#' @param Fi the fishing mortality rate, F, find in stock assessment if don't have more localized estimate
#' @param Linf asymptotic growth rate used in von-Bertallanfy growth equation, can find on fishbase
#' @param k von-bertallanfy growth parameter estimate
#' @param a0 the age at length 0 used in the von-Bertallanfy growth equation
#' @param pW weight length relationship parameter, same as a on fishbase.org but need to divide by 1000 to get in kg not grams
#' @param qW weight length relationship parameter, same as b on fishbase.org
#' @return transient_length is the length of time of the transient duration for the closed population
#' @keywords closed population, population dynamics, Leslie matrix
#' @examples closedpop_time(M=0.2,Fi=0.14,Lfish=25,Linf=37.8,k=0.13,a0=-0.7,maxage=25,pW=9.37e-06,qW=3.172)
#'closed_time()
closedpop_time = function(tf,maxage,Lmat,Lfish,M,Fi, Linf,k,a0,pW,qW) {
a_mat0=(log((Lmat-Linf)/-Linf)/-k)+a0 ##calculate the age at maturity from length
a_mat=round(a_mat0,digits=0)##rounds that age to whole number to input into leslie matrix
a_harv0=(log((Lfish-Linf)/-Linf)/-k)+a0 ##age fished calculated from length fished
a_harv=round(a_harv0,digits=0)##rounds to whole number to input into leslie matrix
amf=a_mat0/a_harv0 ##ratio of age at maturity to age fished
fm_ratio=Fi/M ##ratio fo fishing mortality to natural mortality
# repro and mortality
Ac = c(rep(0,a_mat),rep(1,(maxage-a_mat))) # set age at first reproduction vector
M = c(rep(M,(maxage)),0) # vector of natural mortality at age a
H = c(rep(0,a_harv),rep(Fi,maxage-a_harv)) # vector of fishing mortality at age a
ma = array(0,dim=c(maxage)) # set base array for fecundity at age
La = ma # " length at age
Wa = ma # " weight at age
for (a in 1:(maxage)) {
La[a] = Linf*(1-exp(-k*(a-a0))) # vonBert growth function
Wa[a] = pW*(La[a])^qW # weight-length relationship # Weight-length coefficient kg*cm^q
ma[a] = Wa[a]*Ac[a] # fecundity proportional to weight
}
##Alpha calc function, calculate the alpha (fecundity value)
##needed to keep the unfished population growing at at lambda rate of 1.02
##can adjust that growth rate assumption by changing lambdaGoal
lambdaGoal=1.0
alphafun=function(alpha){ ##function to calculate lambda, which will be used below to calculate correct alpha level
ma_final = alpha*ma ##top row of Leslie matrix
projM_unfished = rbind(ma_final,cbind(diag(c(exp(-M[1:(maxage-1)]))),0)) ##sets up Leslie matrix
lambda = Re(eigen(projM_unfished)$values[1]) ##Lambda calculation is 1st dominant eigenvalue
return(lambda)
}
alphaMin=function(alpha) { ## minimize the alphafun above to find appropriate alpha level
err=(lambdaGoal-alphafun(alpha)) ##error is the difference between the goal lambda and our function of lambda
sum(err^2)/length(err) ##sum of squared error to minimize in optimization function
}
OptimOut=optimize(alphaMin,lower=0,upper=5) ##uses optimization function to determine alpha level
alpha=OptimOut$minimum
##gives correct alpha value
ma_final = alpha*ma # final fecundity values to go across top of Leslie matrix
# construct Leslie matrix by combining ma with diagonal matrix of survival
projM_unfished = rbind(ma_final,cbind(diag(c(exp(-M[1:(maxage-1)]))),0)) # diag creates diagonal matrix with terms exp(-M) for a = 1:maxage-1 and 0 for a = maxage
projM_fished = rbind(ma_final,cbind(diag(c(exp(-(M[1:(maxage-1)]+H[1:(maxage-1)])))),0)) #MPA Leslie matrix
lambda1.unfished = eigen(projM_unfished)$values[1] # lambda1 MPA
lambda2.unfished = eigen(projM_unfished)$values[2] # lambda2 MPA
lambda1.fished = eigen(projM_fished)$values[1] # lambda1 fished
lambda2.fished = eigen(projM_fished)$values[2] # lambda2 fished
# calculate 1st right eigenvector of projM to get SAD, v1
##v1 is the fished stable age distribution
v1 = abs((Re(eigen(projM_fished)$vectors[,1]))) # real portion of eigenvector
N0=rep(1000,maxage)
N0=N0*v1 ##get starting pop vector by multiplying by the fished stable age dist
Nc = matrix(0,tf,maxage) #Initialize vector of population sizes
Nc[1,] = N0 #Put in initial values
for(t in 1:(tf-1)) {
Nc[t+1,1:(maxage)] =projM_unfished%*%Nc[t,]#calculates population changing over time with MPA
}
final.Nc=rowSums(Nc[,a_harv:maxage]) ##include only fished age classes in final abundance ratios
Nratioc=final.Nc/final.Nc[1]
##Calculate weights at length
##Now calculate weights at length
w=pW*(La^qW)
weights=Nc[,a_harv:maxage]%*%w[a_harv:maxage]
Bratio=weights/weights[1]
##Calculate period
P1=2*pi/atan((Im(lambda2.fished)/Re(lambda2.fished)))
##Calculate the rate of convergence, rho
rho=Re(lambda1.fished/abs(lambda2.fished))
##calculate total length of transience
t.t=log(0.01)/-log(rho)
##Calculate theta
w1=abs(Re(eigen(projM_unfished)$vectors[,1]))
theta=acos((N0%*%w1)/(sqrt(sum(N0^2))*sqrt(sum(w1^2))))
ages=seq(1:maxage)
output=list(transient_length=t.t)
return(output)
}
##test function is working with kelp RF parms
#closedpop(maxage=25,Lmat=18,Lfish=25,M=0.2,Fi=0.17, Linf=37.8,k=0.23,a0=-0.7,pW=6.29e-06,qW=3.172)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.