-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathget_lj.R
144 lines (125 loc) · 4.5 KB
/
get_lj.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
## get_lj
# Gets scaled length at metamorphosis
##
#function [lj, lp, lb, info] =
#' Gets scaled length at metamorphosis
#'
#' Type M-acceleration: Isomorph, but V1-morph between vHb and vHj
#' This routine obtaines scaled length at metamorphosis lj given scaled muturity at metamorphosis vHj.
#' The theory behind get_lj, is discussed in the comments to DEB3.
#' If scaled length at birth (third input) is not specified, it is computed (using automatic initial estimate);
#' if it is specified. however, is it just copied to the (third) output.
#' The code assumes vHb < vHj < vHp (see first input).
#'
#' @param p 6-vector with parameters: g, k, l_T, v_H^b, v_H^j, v_H^p; if p is a 5-vector, output lp is empty
#' @param f optional scalar with scaled functional responses (default 1)
#' @param lb0 optional scalar with scaled length at birth
#'
#' @return [lj, lp, lb, info] lj: scalar with scaled length at metamorphosis; * lp: scalar with scaled length at puberty; lb: scalar with scaled length at birth; info: indicator equals 1 if successful
#'
#' @export
#'
#' @importFrom deSolve ode
#'
get_lj=function(p, f, lb0) {
# created at 2010/02/10 by Bas Kooijman,
# modified 2014/03/03 Starrlight Augustine, 2015/01/18 Bas Kooijman
## Syntax
# [lj, lp, lb, info] = <../get_lj.m *get_lj*>(p, f, lb0)
## Description
# Type M-acceleration: Isomorph, but V1-morph between vHb and vHj
# This routine obtaines scaled length at metamorphosis lj given scaled muturity at metamorphosis vHj.
# The theory behind get_lj, is discussed in the comments to DEB3.
# If scaled length at birth (third input) is not specified, it is computed (using automatic initial estimate);
# if it is specified. however, is it just copied to the (third) output.
# The code assumes vHb < vHj < vHp (see first input).
#
# Input
#
# * p: 6-vector with parameters: g, k, l_T, v_H^b, v_H^j, v_H^p
#
# if p is a 5-vector, output lp is empty
#
# * f: optional scalar with scaled functional responses (default 1)
# * lb0: optional scalar with scaled length at birth
#
# Output
#
# * lj: scalar with scaled length at metamorphosis
# * lp: scalar with scaled length at puberty
# * lb: scalar with scaled length at birth
# * info: indicator equals 1 if successful
## Remarks
# Similar to <get_lj1.html get_lj1>, which uses root finding, rather than integration
# Scaled length l = L/ L_m with L_m = kap {p_Am}/ [p_M] where {p_Am} is of embryo,
# because the amount of acceleration is food-dependent so the value after metamorphosis is not a parameter.
# {p_Am} and v increase between maturities v_H^b and v_H^j till
# {p_Am} l_j/ l_b and v l_j/ l_b at metamorphosis.
# After metamorphosis l increases from l_j till l_i = f l_j/ l_b - l_T.
# Scaled length l can thus be larger than 1.
# See <get_lj_foetus.html *get_lj_foetus*> for foetal development.
## Example of use
# get_lj([.5, .1, .1, .01, .2, .3])
n_p = length(p)
# unpack pars
g = p[1] # energy investment ratio
k = p[2] # k_J/ k_M, ratio of maturity and somatic maintenance rate coeff
lT = p[3] # scaled heating length {p_T}/[p_M]
vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_H^b = M_H^b/ {J_EAm} = E_H^b/ {p_Am}
vHj = p[5] # v_H^j = U_H^j g^2 kM^3/ (1 - kap) v^2; U_H^j = M_H^j/ {J_EAm} = E_H^j/ {p_Am}
if (n_p > 5){
vHp = p[6] # v_H^p = U_H^j g^2 kM^3/ (1 - kap) v^2; U_B^j = M_H^j/ {J_EAm}
}
if (!exists('f')) {
f = 1
} else if (is.na(f)){
f = 1
}
# get lb
if (!exists('lb0')){
lb0 = NA
}
if (is.na(lb0)){
lbinfo = get_lb(c(g, k, vHb), f, lb0)
lb=lbinfo[1]
info=lbinfo[2]
} else {
info = 1
lb = lb0
}
if (lT > f - lb) { # no growth is possible at birth
lj = NA
lp = NA
info = 0
}
# no acceleration
if (vHb == vHj){
if (n_p == 5){ # no puberty specified
lbinfo = get_lb(p[c(1, 2, 4)], f, lb0)
lb=lbinfo[1]
info=lbinfo[2]
lj = lb
lp = NA
} # else { # puberty specified
# [lp, lb, info] = get_tp(p([1 2 3 4 6]), f, lb0);
# lj = lb;
# }
}
# get lj
#library(deSolve)
tspan=seq(vHb, vHj, length=100)
lini=c(l=lb)
out = deSolve::ode(lini, tspan, dget_l_V1, parms=c(k, lT, g, f, lb), method="ode45")
vH=out[,1]
lj=out[length(out[,2]),2]
sM = lj/ lb
# get lp
if (n_p > 5){
tspan=seq(vHj, vHp, length=100)
lini=lj
out = ode(lini, tspan, dget_l_ISO, parms=c(k, lT, g, f, sM), method="ode45")
vH=out[,1]
lp=out[length(out[,2]),2]
} else lp = NA
return(c(lj, lp, lb, info))
}