#'@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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.