R/09_find_msyr.R

In mcsiple/mmrefpoints: Project Marine Mammal Populations and Calculate Reference Points

#' calculate MSYR
#'
#' @description one of the functions needed for getting MNPL. MSYR is the same as FMSY.
#'
#' @param E.start a starting guess for bycatch mortality rate that will result in FMSY (numeric value)
#' @param lh.params a list of life history parameters (S0, S1plus, nages, AgeMat, lmabdaMax, K1plus, and z)
#' @param fmax Max theoretical fecundity (numeric)
#' 
#' @examples
#' # Set parameters
#' S0.w = 0.5; S1plus.w = 0.944; nages.w = 25; K1plus.w = 9000; AgeMat.w = 18 
#' InitDepl.w = 0.9; z.w = 2.39; lambdaMax.w = 1.04
#' lh.params = list(S0 = S0.w,S1plus = S1plus.w, nages = nages.w,
#' AgeMat = AgeMat.w,K1plus=9000,z=z.w,lambdaMax = lambdaMax.w) 
#' # Get number of individuals per recruit in terms of mature individuals (\eqn{N0.w})
#' NPROut <- npr(S0 = S0.w, S1plus = S1plus.w, nages = nages.w, AgeMat = AgeMat.w, E = 0)
#' N0 <- NPROut$npr # mature numbers per recruit
#' # Get number of individuals per recruit in terms of individuals aged 1+ (\eqn{P0.w})
#' P0 <- NPROut$P1r # 1+ nums per recruit
#'
#' fmax <- getfecmax(lambdaMax = lambdaMax.w, S0 = S0.w, S1plus = S1plus.w, AgeMat = AgeMat.w)
#' Fec0 <- 1.0 / N0
#' A <- (fmax - Fec0) / Fec0
#' fmsy<-find_msyr(E.start=0.01, lh.params=lh.params, fmax=fmax)
#' cat("fmsy =",fmsy,"\n")
#' results <- matrix(0,ncol=2,nrow=101)
#' results[,1] <- c(0:100)*fmsy/50
#' for (II in 1:101)
#'  results[II,2] <- ce(S0 = S0.w, S1plus = S1plus.w, nages = nages.w, AgeMat = AgeMat.w,z = z.w,
#'                      E=results[II,1], A = A, P0 = P0, N0 = N0)
#' plot(results[,1],results[,2],xlab="f",ylab="ce",type="l",yaxs="i",ylim=c(0,max(results[,2])*1.2))
#' abline(v=fmsy)
#' @export
find_msyr <- function(E.start, lh.params, fmax) {
  S0.w <- lh.params$S0
  S1plus.w <- lh.params$S1plus
  nages.w <- lh.params$nages
  AgeMat.w <- lh.params$AgeMat
  lambdaMax.w <- lh.params$lambdaMax
  K1plus.w <- lh.params$K1plus
  z.w <- lh.params$z
  
  if(AgeMat.w > nages.w){warning("Age at maturity cannot be larger than plus group age. Change AgeMat or nages.")}
  if(S0.w < 0 | S0.w >= 1){stop("Calf/pup survival must be between 0 and 1.")}
  if(S1plus.w < 0 | S1plus.w >= 1){stop("Adult survival must be between 0 and 1.")}
  if(K1plus.w < 0){stop("Carrying capacity K1plus must be greater than zero.")}

  unex <- npr(S0 = S0.w, S1plus = S1plus.w, nages = nages.w, AgeMat = AgeMat.w, E = 0)
  N0 <- unex$npr
  P0 <- unex$P1r

  Fec0 <- 1.0 / N0
  A.w <- (fmax - Fec0) / Fec0

  search.limit <- get_f(
    f.start = 0.005,
    S0.w = S0.w,
    S1plus.w = S1plus.w,
    nages.w = nages.w,
    AgeMat.w = AgeMat.w,
    InitDepl.w = 0.000001, # 0 causes strange behavior
    z.w = z.w,
    lambdaMax.w = lambdaMax.w,
    N0.w = N0,
    P0.w = P0
  ) * .95
  # cat("Search limit",search.limit,"\n")

  lims <- logit(c(0.00001, search.limit))
  logit.E <- logit(E.start)

  zero.cross <- stats::uniroot(f = get_diff, interval = lims, tol = 1e-7, 
                               S0 = S0.w, S1plus = S1plus.w, nages = nages.w, 
                               AgeMat = AgeMat.w, A = A.w, z = z.w)
  fmsy <- inv_logit(zero.cross$root)

  MSY <- ce(
    S0 = S0.w,
    S1plus = S1plus.w,
    nages = nages.w,
    AgeMat = AgeMat.w,
    z = z.w,
    E = fmsy,
    A = A.w, P0 = P0, N0 = N0
  )
  if (MSY == 0) cat("Warning FMSY is not correct", "\n")
  # cat("FMSY Final",fmsy,zero.cross$f.root,MSY,"\n")
  return(fmsy)
}