-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcarabids_03_function_GDD.R
128 lines (115 loc) · 6.81 KB
/
carabids_03_function_GDD.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
# Calculate degree days - TrenchR GitHub (citation("TrenchR"))
# (wouldn't install correctly - developer Lauren said this might be an issue)
#
#
# @details Calculate degree days
# @description This function allows you to calculate degree days using single or double sine wave and single or double triangulation approximation. Source: http://ipm.ucanr.edu/WEATHER/ddfigindex.html. Double methods assume symmetry, that is that next day's thermal minima is equal to previous day. Double Sine wave approximation of degree days from Allen JC. 1976. A Modified Sine Wave Method for Calculating Degree Days. Environmental Entomology 5:388–396.
# @param T_min Minimum temperature of the day (°C)
# @param T_max Maximum temperature of the day (°C)
# @param LDT lower developmental threshold (°C). # Nufio et al. used 12*C for grasshoppers
# @param UDT upper developmental threshold (°C). #Nufio et al. used 33*C - was never reached at Niwot
# @param method type of method being used. Current choices: "single.sine","double.sine", "single.triangulation" or "double.triangulation".
# @return degree days (°C)
# @keywords degree days
# @family microclimate functions
# @export
# @examples
# \dontrun{
# degree_days(T_min=7, T_max=14, LDT=12, UDT=33, method="single.sine")
# degree_days(T_min=7, T_max=14, LDT=12, UDT=33, method="single.triangulation")
# }
#
degree_days=function(T_min,T_max,LDT=NA,UDT=NA, method="single.sine"){
stopifnot(T_max>=T_min, method %in% c("single.sine","double.sine", "single.triangulation", "double.triangulation"))
#amplitude
alpha=(T_max-T_min)/2
dd = 0
#Single sine calculation
if (method == "single.sine") {
if (T_min >= UDT && T_max > UDT) { # entirely above both thresholds
dd = (T_max - T_min)
} else if ( T_min > LDT && T_max > UDT) { #Intercepted by upper threshold
theta2=asin((UDT-(T_max+T_min)/2)/alpha)
dd = 1 / pi * (((T_max + T_min) / 2 - LDT) * (theta2 + pi / 2) + (UDT - LDT) *
(pi / 2 - theta2) - alpha * cos(theta2))
} else if (T_min < LDT && T_max > UDT ) { #Intercepted by both thresholds
theta2=asin((UDT-(T_max+T_min)/2)/alpha)
theta1=asin((LDT-(T_max+T_min)/2)/alpha)
dd = 1 / pi * (( ((T_max + T_min) / 2) - LDT) * (theta2 - theta1) + alpha * (cos(theta1) -
cos(theta2)) + (UDT - LDT) * (pi / 2 - theta2))
} else if (T_min > LDT && T_max < UDT ) { #Entirely between both thresholds
dd = ((T_max + T_min) / 2) - LDT
} else if (T_min < LDT && T_max > LDT) { # intercepted by LDT
#theta1=asin((LDT-(T_max+T_min)/2)/alpha)
# credit - http://stackoverflow.com/questions/20998460/unexpected-behavior-in-asin
#It's a floating point issue. The way floating point numbers work is that all
#numbers need to be mapped to the nearest one which can be expressed as a finite
#sum of powers of two and this may lead to small inaccuracies in the expected output
#and can be dependent upon how the numbers are calculated
theta1= asin(pmax(-1,pmin(1,(LDT-(T_max+T_min)/2)/alpha)))
dd = 1 / pi * (( ((T_max + T_min) / 2) - LDT) * ( (pi / 2) - theta1) + alpha * cos(theta1))
} else if (T_min < LDT && T_max <= LDT) { # entirely below both thresholds
dd = 0
}
} #end single sine method
#double sine calculation
if (method == "double.sine") {
if (T_min >= LDT && T_max > UDT) { # entirely above both thresholds
dd = (UDT - LDT) / 2
} else if (T_min > LDT && T_max > UDT) { #Intercepted by upper threshold
theta2=asin((UDT-(T_max+T_min)/2)/alpha)
dd = 1 / (2 * pi) * (((T_max + T_min) / 2 - LDT) * (theta2 + pi / 2) + (UDT -
LDT) * (pi / 2 - theta2) - alpha * cos(theta2))
} else if (T_min < LDT && T_max > UDT) { #Intercepted by both thresholds
theta2=asin((UDT-(T_max+T_min)/2)/alpha)
theta1=asin((LDT-(T_max+T_min)/2)/alpha)
dd = 1 / (2 * pi) * (((T_max + T_min) / 2 - LDT) * (theta2 - theta1) + alpha *
(cos(theta1) - cos(theta2)) + (UDT - LDT) * (pi / 2 - theta2))
} else if (T_min > LDT && T_max < UDT) { #Entirely between both thresholds
dd = 0.5 * ((T_max + T_min) / 2 - LDT)
} else if (T_min < LDT && T_max > LDT) { # intercepted by LDT
theta1= asin(pmax(-1,pmin(1,(LDT-(T_max+T_min)/2)/alpha)))
dd = 1 / (2 * pi) * (((T_max + T_min) / 2 - LDT) * (pi / 2 - theta1) + alpha *
cos(theta1))
} else if (T_min < LDT && T_max <= LDT) { # entirely below both thresholds
dd = 0
}
dd= dd*2
} #end double sine method
#Single triangulation - with simplified formula
if (method == "single.triangulation") {
MT = (T_max+T_min)/2
if (T_min >= UDT && T_max > UDT) { # entirely above both thresholds
dd = (UDT - LDT)
} else if ( T_min > LDT && T_max > UDT) { #Intercepted by upper threshold
dd = (MT-LDT)-((T_max-UDT)^2/((T_max-T_min)*2))
} else if (T_min < LDT && T_max > UDT ) { #Intercepted by both thresholds
dd = ((T_max-LDT)^2-(T_max-UDT)^2)/((T_max-T_min)*2)
} else if (T_min > LDT && T_max < UDT ) { #Entirely between both thresholds
dd = MT-LDT
} else if (T_min < LDT && T_max > LDT) { # intercepted by LDT
dd = (T_max-LDT)^2/((T_max-T_min)*2)
} else if (T_min < LDT && T_max <= LDT) { # entirely below both thresholds
dd = 0
}
} #end single triangulation method
#Double triangulation - with simplified formula
if (method == "double.triangulation") {
MT = (T_max+T_min)/2
if (T_min >= UDT && T_max > UDT) { # entirely above both thresholds
dd = (UDT - LDT)/2
} else if ( T_min > LDT && T_max > UDT) { #Intercepted by upper threshold
dd = (MT-LDT)-((T_max-UDT)^2/((T_max-T_min)*4))
} else if (T_min < LDT && T_max > UDT ) { #Intercepted by both thresholds
dd = ((T_max-LDT)^2-(T_max-UDT)^2)/((T_max-T_min)*4)
} else if (T_min > LDT && T_max < UDT ) { #Entirely between both thresholds
dd = (MT/4)-(LDT/2)
} else if (T_min < LDT && T_max > LDT) { # intercepted by LDT
dd = (T_max-LDT)^2/((T_max-T_min)*4)
} else if (T_min < LDT && T_max <= LDT) { # entirely below both thresholds
dd = 0
}
dd= dd*2
} #end double triangulation method
return(round(dd,2))
}