R/get_tb.R

In Echinophoria/DEB_output: Dynamic Energy Budget exploration tools

#'@export
get_tb <- function (p, f) {
  
  ##### Obtaining lb and tb. Scaled length and age at birth
  # input --> p: a list containing all the DEB parameters
  # f --> scaled functional response for which these parameters should be calculated
  
  
  if (f>1 || f<=0){
    warning("f must be a value between 0 and 1")
    stop()
  }
  
  if (length(p)< 3){
    warning ("p vector must include at least g, k, v_Hb (same order)")
    stop()
  }else{
    g = p$g
    k = p$k
    v_Hb = p$v_Hb
  }
  
  
  ##### calculating lb --> scaled length at birth ####
  
  
  n = 1000 + round(1000 * max(0, k - 1))
  xb = g/ (g + f)
  xb3 = xb ^ (1/3)
  x = seq(1e-6, xb, length=n)
  dx = xb/ n
  x3 = x ^ (1/3)
  lb = v_Hb^(1/3) # exact solution for k==1, to be used as starting value.
  b = beta0(x, xb)/ (3 * g)
  t0 = xb * g * v_Hb
  i = 0
  norm = 1 # make sure that we start Newton Raphson procedure
  ni = 100 # max number of iterations
  
  while ((i < ni) & (norm > 1e-8)) {
    l = x3 / (xb3/ lb - b);
    s = (k - x) / (1 - x) * l/ g / x;
    y = exp( - dx * cumsum(s))
    vb = y[n];
    r = (g + l)
    rv = r / y;
    t = t0/ lb^3/ vb - dx * sum(rv);
    dl = xb3/ lb^2 * l^2 / x3
    dlnl = dl / l;
    dv = y * exp( - dx * cumsum(s * dlnl));
    dvb = dv[n]
    dlnv = dv / y
    dlnvb = dlnv[n];
    dr = dl 
    dlnr = dr / r;
    dt = - t0/ lb^3/ vb * (3/ lb + dlnvb) - dx * sum((dlnr - dlnv) * rv);
    lb = lb - t/ dt # Newton Raphson step
    lb = Re(lb)
    norm = Re(t^2)
    i = i + 1   
  }
  
  ##### obtaining tb --> scaled age at birth ####
  
  dget_tb <- function(x, ab, xb){
    f = x^(-2/3) / (1 - x) / (ab - beta0(x, xb))
    return(f)
  }  # get scaled age at birth
  
  xb = g/ (f + g);
  ab = 3 * g * xb^(1/ 3)/ lb
  tb = 3 * quad(dget_tb, 1e-15, xb,  tol=1e-15, trace = FALSE, ab, xb);
  tb = Re(tb)
  
  
  ## create output
  
  ltr = list()  #named list containign the sought values (lb,tb)
  ltr["lb"]=lb   # scaled length at birth
  ltr["tb"]=tb   # scaled age at birth

  
  return(ltr)
}