-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathget_tj.R
173 lines (150 loc) · 6.08 KB
/
get_tj.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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
## get_tj
# Gets scaled age at metamorphosis
##
# function [tj, tp, tb, lj, lp, lb, li, rj, rB, info] =
get_tj = function(p, f, lb0=NA){
# created at 2011/04/25 by Bas Kooijman,
# modified 2014/03/03 Starrlight Augustine, 2015/01/18 Bas Kooijman
## Syntax
# [tj, tp, tb, lj, lp, lb, li, rj, rB, info] = <../get_tj.m *get_tj*> (p, f, lb0)
## Description
# Obtains scaled ages at metamorphosis, puberty, birth and the scaled lengths at these ages;
# Food density is assumed to be constant.
# Multiply the result with the somatic maintenance rate coefficient to arrive at unscaled ages.
# Metabolic acceleration occurs between birth and metamorphosis, see also get_ts.
# Notice j-p-b sequence in output, due to the name of the routine
#
# Input
#
# * p: 5 or 6-vector with parameters: g, k, l_T, v_H^b, v_H^j v_H^p
#
# Last value is optional. If ommitted: outputs tp and lp are empty
#
# * f: optional scalar with functional response (default f = 1)
# * lb0: optional scalar with scaled length at birth
# or optional 2-vector with scaled length, l, and scaled maturity, vH
# for a juvenile that is now exposed to f, but previously at another f
#
# Output
#
# * tj: scaled with age at metamorphosis \tau_j = a_j k_M
#
# if length(lb0)==2, tj is the scaled time till metamorphosis
#
# * tp: scaled with age at puberty \tau_p = a_p k_M
#
# if length(lb0)==2, tp is the scaled time till puberty
#
# * tb: scaled with age at birth \tau_b = a_b k_M
# * lj: scaled length at end of V1-stage
# * lp: scaled length at puberty
# * lb: scaled length at birth
# * li: ultimate scaled length
# * rj: scaled exponential growth rate between s and j
# * rB: scaled von Bertalanffy growth rate between b and s and between j and i
# * info: indicator equals 1 if successful, 0 otherwise
## Remarks
# See <get_tj_foetus.html 8get_tj_foetus*> in case of foetal development
## Example of use
# get_tj([.5, .1, 0, .01, .05, .2])
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]Lm
vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_B^b = M_H^b/ {J_EAm}
vHj = p[5] # v_H^j = U_H^j g^2 kM^3/ (1 - kap) v^2; U_B^j = M_H^j/ {J_EAm}
if (n_p > 5){
vHp = p[6]} else vHp = 0 # v_H^p = U_H^p g^2 kM^3/ (1 - kap) v^2; U_B^p = M_H^p/ {J_EAm}
if (!exists('f')) {
f = 1
}
if (is.na(f)){
f = 1
}
if (!exists('lb0')){
lb0 = NA
}
# no acceleration
if (vHb == vHj){
if (vHp == 0){ # no puberty specified
tblbinfo=get_tb(p[c(1:3)], f, lb0)
tj = tblbinfo[1]
tp = NA
lp = NA
lj = tblbinfo[2]
li = f - lT
rj = 0
rB = 1/ 3/ (1 + f/ g)
} # else { # puberty specified
# [tp, tb, lp, lb, info] = get_tp(p([1 2 3 4 6]), f, lb0);
# tj = tb; lj = lb; li = f - lT; rj = 0; rB = 1/ 3/ (1 + f/ g);
# }
}
# maintenance ratio k = 1: maturity thresholds coincide with length thresholds
if (k == 1 & f * (f - lT)^2 > vHp * k){ # constant maturity density, reprod possible
lb = vHb^(1/3) # scaled length at birth
tblbinfo = get_tb(p[c(1,2,4)], f, lb) # scaled age at birth
tb = tblbinfo[1]
lj = vHj^(1/3) # scaled length at metamorphosis
sM = lj/ lb # acceleration factor
rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j
tj = tb + (log(sM)) * 3/ rj # scaled age at metamorphosis
lp = vHp^(1/3) # scaled length at puberty
li = f * sM - lT # scaled ultimate length
rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i
tp = tj + (log ((li - lj)/ (li - lp)))/ rB # scaled age at puberty
info = 1
} else if (k == 1 & f * (f - lT)^2 > vHj * k) { # constant maturity density, metam possible
lb = vHb^(1/3) # scaled length at birth
tblbinfo = get_tb(p[c(1,2,4)], f, lb) # scaled age at birth
tb = tblbinfo[1]
lj = vHj^(1/3) # scaled length at metamorphosis
sM = lj/ lb # acceleration factor
rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j
tj = tb + (log(sM)) * 3/ rj # scaled age at metamorphosis
lp = vHp^(1/3) # scaled length at puberty
li = f * sM - lT # scaled ultimate length
rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i
tp = 1e20 # scaled age at puberty
info = 1
}
if (is.na(lb0)){
tblbinfo = get_tb (p[c(1, 2, 4)], f, lb0)
} else {tblbinfo = get_tb (p[C(1, 2, 4)], f, lb0[1]) }
tb=tblbinfo[1]
lb=tblbinfo[2]
info_tb=tblbinfo[3]
ljlplbinfo = get_lj(p, f, lb)
lj=ljlplbinfo[1]
lp=ljlplbinfo[2]
lb=ljlplbinfo[3]
info_tj=ljlplbinfo[3]
sM = lj/lb # acceleration factor
rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j
tj = tb + log(sM) * 3/ rj # scaled age at metamorphosis
rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i
li = f * sM - lT # scaled ultimate length
if (is.na(lp)) { # length(p) < 6
tp = NA
} else if (li <= lp) { # reproduction is not possible
tp = 1e20; # tau_p is never reached
lp = 1; # lp is nerver reached
} else { # reproduction is possible
if (length(lb0) != 2) { # lb0 is absent, empty or a scalar
tp = tj + log((li - lj)/ (li - lp))/ rB
} else{ # lb0 = l and t for a juvenile
tb = NaN
l = lb0[1]
tp = log((li - l)/ (li - lp))/ rB;
}
}
if (is.na(tp)) info = info_tb
info = min(info_tb, info_tj)
if (!is.finite(tp) || !is.finite(tj)){ # tj and tp must be real and positive
info = 0
} else if (tp < 0 || tj < 0){
info = 0
}
return(c(tj, tp, tb, lj, lp, lb, li, rj, rB, info))
}