R/ode_functs.R

In Echinophoria/DEB_output: Dynamic Energy Budget exploration tools

beta0 <- function (x0,x1){
  x03 = x0^(1/3); x13 = x1^(1/3); a3 = sqrt(3)
  f1 = -3 * x13 + a3 * atan((1+2 * x13)/a3) - log (as.complex(x13 - 1))+
    log(1 + x13 + x13^2)/2
  f0 = -3 * x03 + a3 * atan((1+2 * x03)/a3) - log (as.complex(x03 - 1))+
    log(1 + x03 + x03^2)/2
  f = f1 - f0
  return(f)
} 

#'@export
dget_l_ISO <- function (vH, l, parms) {  # growth for isometric adult)
  if (length(parms)< 3){
    warning ("parms vector must include at least k,g,f (same order)")
    stop()
  }
  if (length(parms)==3){
    k = parms[1]
    g = parms[2]
    f = parms[3]
    sM = 1
    }else{
      k = parms[1]
      g = parms[2]
      f = parms[3]
      sM = parms[4]
      }

  r = g * (f * sM - l)/ l / (f + g)
  dl = l * r / 3
  dvH = f * l^2 * (sM - l * r / g) - k * vH
  dl = dl / dvH
  return(list(dl))
}  # get length asuming isomorphic growth (VB)

#'@export
dget_l_V1 <- function (vH, l, parms){   # growth for v1 morph (for length at metamorphosis)
  k = parms[1]
  g= parms[2]
  f=parms[3]
  lref=parms[4]
  r=((g*(f - lref))/lref)/(g+f) # specific growth rate
  dl = r * l/3 # d/dt l 
  dvH = f * l^3 * (1/lref - r/g) - k * vH; # d/dt vH
  dl = dl/dvH;
  return(list(dl))
}   # get length at metamorphosis asumes V1-morph growth

#'@export
dget_LUH <- function (a, LUH, p) {
  kap = p$kap
  v = p$v
  kJ= p$k_J
  g = p$g
  Lm = p$L_m
  L = LUH[1]  # cm, structural length
  U = LUH[2]  # d cm^2, scaled reserves M_E/{J_EAm}
  H = LUH[3]  # d cm^2, scaled maturity M_H/{J_EAm}
  eL3 = U * v
  gL3 = g * L^3
  SC = L^2 * (1 + L/(g*Lm))*g*eL3/(gL3 + eL3)
  dL = v * (eL3 - L^4 / Lm)/(3 * (eL3+gL3))
  dU = -SC
  dH = (1-kap)* SC - kJ * H
  return(list(c(dL,dU,dH)))
}  # get structural length and reserves at age during embryonic development

#'@export
get_uE0 <- function (p, lb, f){ # initial scaled reserves
  f=f        # -, scaled functional response for the mother
  lb=lb      # cm, scaled length at birth
  g = p$g    # -, energy investment ratio
  v = p$v    # cm2/d, energy conductance
  k_M=p$k_M  # 1/d, somatic maintenance rate coefficient
  
  
  xb = g/ (f + g)
  u0 = Re( (3 * g/ (3 * g * xb^(1/ 3)/ lb - beta0(0, xb)))^3)
  uE0 = u0 * v^2/ g^2/ k_M^3 # d.cm^2, initial scaled reserve
  return(uE0)
}

#'@export
random.sample<-function(n,mean,sd,lowerBound,upperBound){
  range <- upperBound - lowerBound
  m <- (mean-lowerBound) / range #mapping mean to 0-1 range
  s <- sd / range #mapping sd to 0-1 range
  a <- (m^2 - m^3 - m*s^2)/s^2 #calculating alpha for rbeta 
  b <- (m-2*m^2+m^3-s^2+m*s^2)/s^2 #calculating beta for rbeta
  data <- rbeta(n,a,b)  #generating data
  data <- lowerBound + data * range #remaping to given bounds
  return(data)
}