Mercurial > repos > ecology > srs_preprocess_s2

comparison prosail-master/R/Lib_PROSAIL.R @ 0:a990ea26442f draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/Marie59/Sentinel_2A/srs_tools commit b32737c1642aa02cc672534e42c5cb4abe0cd3e7
author ecology
date Mon, 09 Jan 2023 13:54:29 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:a990ea26442f
1 # ============================================================================= =
2 # prosail
3 # Lib_PROSAIL.R
4 # ============================================================================= =
5 # PROGRAMMERS:
6 # Jean-Baptiste FERET <jb.feret@teledetection.fr >
7 # Florian de BOISSIEU <fdeboiss@gmail.com >
8 # copyright 2019 / 11 Jean-Baptiste FERET
9 # ============================================================================= =
10 # This Library includes functions dedicated to PROSAIL simulation
11 # SAIL versions available are 4SAIL and 4SAIL2
12 # ============================================================================= =
13
14 #" computes bidirectional reflectance factor based on outputs from PROSAIL and sun position
15 #"
16 #" The direct and diffuse light are taken into account as proposed by:
17 #" Francois et al. (2002) conversion of 400-1100 nm vegetation albedo
18 #" measurements into total shortwave broadband albedo using a canopy
19 #" radiative transfer model, Agronomie
20 #" Es = direct
21 #" Ed = diffuse
22 #"
23 #" @param rdot numeric. Hemispherical-directional reflectance factor in viewing direction
24 #" @param rsot numeric. Bi-directional reflectance factor
25 #" @param tts numeric. Solar zenith angle
26 #" @param specatm_sensor list. direct and diffuse radiation for clear conditions
27 #" @return brf numeric. Bidirectional reflectance factor
28 #" @export
29 compute_brf <- function(rdot, rsot, tts, specatm_sensor) {
30
31 ############################## #
32 ## direct / diffuse light ##
33 ############################## #
34 es <- specatm_sensor$Direct_Light
35 ed <- specatm_sensor$Diffuse_Light
36 rd <- pi / 180
37 skyl <- 0.847 - 1.61 * sin((90 - tts) * rd) + 1.04 * sin((90 - tts) * rd) * sin((90 - tts) * rd) # diffuse radiation (Francois et al., 2002)
38 pardiro <- (1 - skyl) * es
39 pardifo <- skyl * ed
40 brf <- (rdot * pardifo + rsot * pardiro) / (pardiro + pardifo)
41 return(brf)
42 }
43
44 #" Performs PROSAIL simulation based on a set of combinations of input parameters
45 #" @param spec_sensor list. Includes optical constants required for PROSPECT
46 #" refractive index, specific absorption coefficients and corresponding spectral bands
47 #" @param input_prospect list. PROSPECT input variables
48 #" @param n numeric. Leaf structure parameter
49 #" @param chl numeric. chlorophyll content (microg.cm-2)
50 #" @param car numeric. carotenoid content (microg.cm-2)
51 #" @param ant numeric. anthocyain content (microg.cm-2)
52 #" @param brown numeric. brown pigment content (Arbitrary units)
53 #" @param ewt numeric. Equivalent Water Thickness (g.cm-2)
54 #" @param lma numeric. Leaf Mass per Area (g.cm-2)
55 #" @param prot numeric. protein content (g.cm-2)
56 #" @param cbc numeric. nonprotcarbon-based constituent content (g.cm-2)
57 #" @param alpha numeric. Solid angle for incident light at surface of leaf (simulation of roughness)
58 #" @param typelidf numeric. Type of leaf inclination distribution function
59 #" @param lidfa numeric.
60 #" if typelidf == 1, controls the average leaf slope
61 #" if typelidf == 2, corresponds to average leaf angle
62 #" @param lidfb numeric.
63 #" if typelidf == 1, unused
64 #" if typelidf == 2, controls the distribution"s bimodality
65 #" @param lai numeric. Leaf Area Index
66 #" @param q numeric. Hot Spot parameter
67 #" @param tts numeric. Sun zeith angle
68 #" @param tto numeric. Observer zeith angle
69 #" @param psi numeric. Azimuth Sun / Observer
70 #" @param rsoil numeric. Soil reflectance
71 #" @param fraction_brown numeric. Fraction of brown leaf area
72 #" @param diss numeric. Layer dissociation factor
73 #" @param cv numeric. vertical crown cover percentage
74 #" = % ground area covered with crowns as seen from nadir direction
75 #" @param zeta numeric. Tree shape factor
76 #" = ratio of crown diameter to crown height
77 #" @param sailversion character. choose between 4SAIL and 4SAIL2
78 #" @param brownvegetation list. Defines optical properties for brown vegetation, if not nULL
79 #" - WVL, reflectance, Transmittance
80 #" - Set to nULL if use PROSPECT to generate it
81 #"
82 #" @return list. rdot, rsot, rddt, rsdt
83 #" rdot: hemispherical-directional reflectance factor in viewing direction
84 #" rsot: bi-directional reflectance factor
85 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
86 #" rddt: bi-hemispherical reflectance factor
87 #" @import prospect
88 #" @export
89 pro4sail <- function(spec_sensor, input_prospect = nULL, n = 1.5, chl = 40.0,
90 car = 8.0, ant = 0.0, brown = 0.0, ewt = 0.01,
91 lma = 0.008, prot = 0.0, cbc = 0.0, alpha = 40.0,
92 typelidf = 2, lidfa = nULL, lidfb = nULL, lai = nULL,
93 q = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
94 fraction_brown = 0.0, diss = 0.0, cv = 1, zeta = 1,
95 sailversion = "4SAIL", brownvegetation = nULL) {
96
97 ############################ #
98 # LEAF OPTICAL PROPERTIES ##
99 ############################ #
100 if (is.null(input_prospect)) {
101 input_prospect <- data.frame("chl" = chl, "car" = car, "ant" = ant, "brown" = brown, "ewt" = ewt,
102 "lma" = lma, "prot" = prot, "cbc" = cbc, "n" = n, "alpha" = alpha)
103 }
104 greenvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
105 n = input_prospect$n[1],
106 chl = input_prospect$chl[1],
107 car = input_prospect$car[1],
108 ant = input_prospect$ant[1],
109 brown = input_prospect$brown[1],
110 ewt = input_prospect$ewt[1],
111 lma = input_prospect$lma[1],
112 prot = input_prospect$prot[1],
113 cbc = input_prospect$cbc[1],
114 alpha = input_prospect$alpha[1])
115
116 if (sailversion == "4SAIL2") {
117 # 4SAIL2 requires one of the following combination of input parameters
118 # Case #1: valid optical properties for brown vegetation
119 if (!is.null(brownvegetation)) {
120 # need to define reflectance and Transmittance for brownvegetation
121 if (length(grep("reflectance", names(brownvegetation))) == 0 || length(grep("Transmittance", names(brownvegetation))) == 0) {
122 message("Please define brownvegetation as a list including reflectance and Transmittance")
123 stop()
124 }
125 # check if spectral domain for optical properties of brown vegetation match
126 # with spectral domain for optical properties of green vegetation
127 if (length(setdiff(spec_sensor$lambda, brownvegetation$wvl)) > 0) {
128 message("Please define same spectral domain for brownvegetation and SpecPROSPECT")
129 stop()
130 }
131 if (length(unique(lengths(input_prospect))) == 1) {
132 if (!unique(lengths(input_prospect)) == 1) {
133 message("brownvegetation defined along with multiple leaf chemical properties")
134 message("Only first set of leaf chemical properties will be used to simulate green vegetation")
135 }
136 }
137 # if no leaf optical properties brown vegetation defined
138 } else if (is.null(brownvegetation)) {
139 # if all PROSPECT input parameters have the same length
140 if (length(unique(lengths(input_prospect))) == 1) {
141 # if all PROSPECT input parameters are unique (no possibility to simulate 2 types of leaf optics)
142 if (unique(lengths(input_prospect)) == 1) {
143 # if fraction_brown set to 0, then assign green vegetation optics to brown vegetation optics
144 if (fraction_brown == 0) {
145 brownvegetation <- greenvegetation
146 # else run 4SAIL
147 } else {
148 message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
149 message("Currently one set is defined. will run 4SAIL instead of 4SAIL2")
150 sailversion <- "4SAIL"
151 }
152 # if all PROSPECT parameters have at least 2 elements
153 } else if (unique(lengths(input_prospect)) >= 2) {
154 # compute leaf optical properties
155 brownvegetation <- prospect::PROSPECT(SpecPROSPECT = spec_sensor,
156 n = input_prospect$n[2],
157 chl = input_prospect$chl[2],
158 car = input_prospect$car[2],
159 ant = input_prospect$ant[2],
160 brown = input_prospect$brown[2],
161 ewt = input_prospect$ewt[2],
162 lma = input_prospect$lma[2],
163 prot = input_prospect$prot[2],
164 cbc = input_prospect$cbc[2],
165 alpha = input_prospect$alpha[2])
166 if (unique(lengths(input_prospect)) > 2) {
167 message("4SAIL2 needs two sets of optical properties for green and brown vegetation")
168 message("Currently more than 2 sets are defined. will only use the first 2")
169 }
170 }
171 }
172 }
173 }
174 if (sailversion == "4SAIL") {
175 if (length(unique(lengths(input_prospect))) == 1) {
176 if (unique(lengths(input_prospect)) > 1) {
177 message("4SAIL needs only one set of optical properties")
178 message("Currently more than one set of leaf chemical constituents is defined.")
179 message("Will run 4SAIL with the first set of leaf chemical constituents")
180 }
181 }
182 }
183
184 if (sailversion == "4SAIL") {
185 # run 4SAIL
186 ref <- foursail(leafoptics = greenvegetation,
187 typelidf, lidfa, lidfb, lai, q, tts, tto, psi, rsoil)
188 } else if (sailversion == "4SAIL2") {
189 # run 4SAIL2
190 ref <- foursail2(leafgreen = greenvegetation, leafbrown = brownvegetation,
191 typelidf, lidfa, lidfb, lai, q, tts, tto, psi, rsoil,
192 fraction_brown, diss, cv, zeta)
193 }
194 return(ref)
195 }
196
197 #" Performs PROSAIL simulation based on a set of combinations of input parameters
198 #" @param leafoptics list. Includes leaf optical properties (reflectance and transmittance)
199 #" and corresponding spectral bands
200 #" @param typelidf numeric. Type of leaf inclination distribution function
201 #" @param lidfa numeric.
202 #" if typelidf == 1, controls the average leaf slope
203 #" if typelidf == 2, corresponds to average leaf angle
204 #" @param lidfb numeric.
205 #" if typelidf == 1, unused
206 #" if typelidf == 2, controls the distribution"s bimodality
207 #" @param lai numeric. Leaf Area Index
208 #" @param q numeric. Hot Spot parameter
209 #" @param tts numeric. Sun zeith angle
210 #" @param tto numeric. Observer zeith angle
211 #" @param psi numeric. Azimuth Sun / Observer
212 #" @param rsoil numeric. Soil reflectance
213 #"
214 #" @return list. rdot, rsot, rddt, rsdt
215 #" rdot: hemispherical-directional reflectance factor in viewing direction
216 #" rsot: bi-directional reflectance factor
217 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
218 #" rddt: bi-hemispherical reflectance factor
219 #" @export
220
221 foursail <- function(leafoptics, typelidf = 2, lidfa = nULL, lidfb = nULL, lai = nULL,
222 q = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL) {
223
224 ############################## #
225 # LEAF OPTICAL PROPERTIES ##
226 ############################## #
227 rho <- leafoptics$Reflectance
228 tau <- leafoptics$Transmittance
229
230 # Geometric quantities
231 rd <- pi / 180
232 cts <- cos(rd * tts)
233 cto <- cos(rd * tto)
234 ctscto <- cts * cto
235 tants <- tan(rd * tts)
236 tanto <- tan(rd * tto)
237 cospsi <- cos(rd * psi)
238 dso <- sqrt(tants * tants + tanto * tanto - 2. * tants * tanto * cospsi)
239
240 # Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
241 if (typelidf == 1) {
242 foliar_distrib <- dladgen(lidfa, lidfb)
243 lidf <- foliar_distrib$lidf
244 litab <- foliar_distrib$litab
245
246 } else if (typelidf == 2) {
247 foliar_distrib <- campbell(lidfa)
248 lidf <- foliar_distrib$lidf
249 litab <- foliar_distrib$litab
250 }
251
252 # angular distance, compensation of shadow length
253 # Calculate geometric factors associated with extinction and scattering
254 # Initialise sums
255 ks <- 0
256 ko <- 0
257 bf <- 0
258 sob <- 0
259 sof <- 0
260
261 # Weighted sums over LIDF
262 na <- length(litab)
263 for (i in 1:na) {
264 ttl <- litab[i] # leaf inclination discrete values
265 ctl <- cos(rd * ttl)
266 # SAIL volume scattering phase function gives interception and portions to be
267 # multiplied by rho and tau
268 resvolscatt <- volscatt(tts, tto, psi, ttl)
269 chi_s <- resvolscatt$chi_s
270 chi_o <- resvolscatt$chi_o
271 frho <- resvolscatt$frho
272 ftau <- resvolscatt$ftau
273
274 #********************************************************************************
275 #* SUITS SYSTEM coEFFICIEnTS
276 #*
277 #* ks : Extinction coefficient for direct solar flux
278 #* ko : Extinction coefficient for direct observed flux
279 #* att : Attenuation coefficient for diffuse flux
280 #* sigb : Backscattering coefficient of the diffuse downward flux
281 #* sigf : Forwardscattering coefficient of the diffuse upward flux
282 #* sf : Scattering coefficient of the direct solar flux for downward diffuse flux
283 #* sb : Scattering coefficient of the direct solar flux for upward diffuse flux
284 #* vf : Scattering coefficient of upward diffuse flux in the observed direction
285 #* vb : Scattering coefficient of downward diffuse flux in the observed direction
286 #* w : Bidirectional scattering coefficient
287 #********************************************************************************
288
289 # Extinction coefficients
290 ksli <- chi_s / cts
291 koli <- chi_o / cto
292
293 # Area scattering coefficient fractions
294 sobli <- frho * pi / ctscto
295 sofli <- ftau * pi / ctscto
296 bfli <- ctl * ctl
297 ks <- ks + ksli * lidf[i]
298 ko <- ko + koli * lidf[i]
299 bf <- bf + bfli * lidf[i]
300 sob <- sob + sobli * lidf[i]
301 sof <- sof + sofli * lidf[i]
302 }
303
304 # Geometric factors to be used later with rho and tau
305 sdb <- 0.5 * (ks + bf)
306 sdf <- 0.5 * (ks - bf)
307 dob <- 0.5 * (ko + bf)
308 dof <- 0.5 * (ko - bf)
309 ddb <- 0.5 * (1. + bf)
310 ddf <- 0.5 * (1. - bf)
311
312 # Here rho and tau come in
313 sigb <- ddb * rho + ddf * tau
314 sigf <- ddf * rho + ddb * tau
315 att <- 1 - sigf
316 m2 <- (att + sigb) * (att - sigb)
317 m2[which(m2 <= 0)] <- 0
318 m <- sqrt(m2)
319
320 sb <- sdb * rho + sdf * tau
321 sf <- sdf * rho + sdb * tau
322 vb <- dob * rho + dof * tau
323 vf <- dof * rho + dob * tau
324 w <- sob * rho + sof * tau
325
326 # Here the LAI comes in
327 # Outputs for the case LAI = 0
328 if (lai < 0) {
329 tss <- 1
330 too <- 1
331 tsstoo <- 1
332 rdd <- 0
333 tdd <- 1
334 rsd <- 0
335 tsd <- 0
336 rdo <- 0
337 tdo <- 0
338 rso <- 0
339 rsos <- 0
340 rsod <- 0
341
342 rddt <- rsoil
343 rsdt <- rsoil
344 rdot <- rsoil
345 rsodt <- 0 * rsoil
346 rsost <- rsoil
347 rsot <- rsoil
348 } else {
349 # Other cases (LAI > 0)
350 e1 <- exp(-m * lai)
351 e2 <- e1 * e1
352 rinf <- (att - m) / sigb
353 rinf2 <- rinf * rinf
354 re <- rinf * e1
355 denom <- 1. - rinf2 * e2
356
357 j1ks <- jfunc1(ks, m, lai)
358 j2ks <- jfunc2(ks, m, lai)
359 j1ko <- jfunc1(ko, m, lai)
360 j2ko <- jfunc2(ko, m, lai)
361
362 ps <- (sf + sb * rinf) * j1ks
363 qs <- (sf * rinf + sb) * j2ks
364 pv <- (vf + vb * rinf) * j1ko
365 qv <- (vf * rinf + vb) * j2ko
366
367 rdd <- rinf * (1. - e2) / denom
368 tdd <- (1. - rinf2) * e1 / denom
369 tsd <- (ps - re * qs) / denom
370 rsd <- (qs - re * ps) / denom
371 tdo <- (pv - re * qv) / denom
372 rdo <- (qv - re * pv) / denom
373
374 tss <- exp(-ks * lai)
375 too <- exp(-ko * lai)
376 z <- jfunc3(ks, ko, lai)
377 g1 <- (z - j1ks * too) / (ko + m)
378 g2 <- (z - j1ko * tss) / (ks + m)
379
380 tv1 <- (vf * rinf + vb) * g1
381 tv2 <- (vf + vb * rinf) * g2
382 t1 <- tv1 * (sf + sb * rinf)
383 t2 <- tv2 * (sf * rinf + sb)
384 t3 <- (rdo * qs + tdo * ps) * rinf
385
386 # Multiple scattering contribution to bidirectional canopy reflectance
387 rsod <- (t1 + t2 - t3) / (1. - rinf2)
388
389 # Treatment of the hotspot-effect
390 alf <- 1e6
391 # Apply correction 2 / (K + k) suggested by F.-M. Breon
392 if (q > 0) {
393 alf <- (dso / q) * 2. / (ks + ko)
394 }
395 if (alf > 200) {
396 # inserted H. Bach 1 / 3 / 04
397 alf <- 200
398 }
399 if (alf == 0) {
400 # The pure hotspot - no shadow
401 tsstoo <- tss
402 sumint <- (1 - tss) / (ks * lai)
403 } else {
404 # Outside the hotspot
405 fhot <- lai * sqrt(ko * ks)
406 # Integrate by exponential Simpson method in 20 steps
407 # the steps are arranged according to equal partitioning
408 # of the slope of the joint probability function
409 x1 <- 0
410 y1 <- 0
411 f1 <- 1
412 fint <- (1. - exp(-alf)) * 0.05
413 sumint <- 0
414 for (i in 1:20) {
415 if (i < 20) {
416 x2 <- -log(1. - i * fint) / alf
417 } else {
418 x2 <- 1
419 }
420 y2 <- -(ko + ks) * lai * x2 + fhot * (1. - exp(-alf * x2)) / alf
421 f2 <- exp(y2)
422 sumint <- sumint + (f2 - f1) * (x2 - x1) / (y2 - y1)
423 x1 <- x2
424 y1 <- y2
425 f1 <- f2
426 }
427 tsstoo <- f1
428 }
429 # Bidirectional reflectance
430 # Single scattering contribution
431 rsos <- w * lai * sumint
432 # Total canopy contribution
433 rso <- rsos + rsod
434 # Interaction with the soil
435 dn <- 1. - rsoil * rdd
436 # rddt: bi-hemispherical reflectance factor
437 rddt <- rdd + tdd * rsoil * tdd / dn
438 # rsdt: directional-hemispherical reflectance factor for solar incident flux
439 rsdt <- rsd + (tsd + tss) * rsoil * tdd / dn
440 # rdot: hemispherical-directional reflectance factor in viewing direction
441 rdot <- rdo + tdd * rsoil * (tdo + too) / dn
442 # rsot: bi-directional reflectance factor
443 rsodt <- rsod + ((tss + tsd) * tdo + (tsd + tss * rsoil * rdd) * too) * rsoil / dn
444 rsost <- rsos + tsstoo * rsoil
445 rsot <- rsost + rsodt
446 }
447 my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt)
448 return(my_list)
449 }
450
451 #" Performs pro4sail2 simulation based on a set of combinations of input parameters
452 #" @param leafgreen list. includes relfectance and transmittance for vegetation #1 (e.g. green vegetation)
453 #" @param leafbrown list. includes relfectance and transmittance for vegetation #2 (e.g. brown vegetation)
454 #" @param typelidf numeric. Type of leaf inclination distribution function
455 #" @param lidfa numeric.
456 #" if typelidf == 1, controls the average leaf slope
457 #" if typelidf == 2, corresponds to average leaf angle
458 #" @param lidfb numeric.
459 #" if typelidf == 1, unused
460 #" if typelidf == 2, controls the distribution"s bimodality
461 #" @param lai numeric. Leaf Area Index
462 #" @param hot numeric. Hot Spot parameter = ratio of the correlation length of leaf projections in the horizontal plane and the canopy height (doi:10.1016 / j.rse.2006.12.013)
463 #" @param tts numeric. Sun zeith angle
464 #" @param tto numeric. Observer zeith angle
465 #" @param psi numeric. Azimuth Sun / Observer
466 #" @param rsoil numeric. Soil reflectance
467 #" @param fraction_brown numeric. Fraction of brown leaf area
468 #" @param diss numeric. Layer dissociation factor
469 #" @param cv numeric. vertical crown cover percentage
470 #" = % ground area covered with crowns as seen from nadir direction
471 #" @param zeta numeric. Tree shape factor
472 #" = ratio of crown diameter to crown height
473 #"
474 #" @return list. rdot, rsot, rddt, rsdt
475 #" rdot: hemispherical-directional reflectance factor in viewing direction
476 #" rsot: bi-directional reflectance factor
477 #" rsdt: directional-hemispherical reflectance factor for solar incident flux
478 #" rddt: bi-hemispherical reflectance factor
479 #" alfast: canopy absorptance for direct solar incident flux
480 #" alfadt: canopy absorptance for hemispherical diffuse incident flux
481 #" @export
482
483 foursail2 <- function(leafgreen, leafbrown,
484 typelidf = 2, lidfa = nULL, lidfb = nULL,
485 lai = nULL, hot = nULL, tts = nULL, tto = nULL, psi = nULL, rsoil = nULL,
486 fraction_brown = 0.5, diss = 0.5, cv = 1, zeta = 1) {
487
488 # This version does not include non-Lambertian soil properties.
489 # original codes do, and only need to add the following variables as input
490 rddsoil <- rdosoil <- rsdsoil <- rsosoil <- rsoil
491
492 # Geometric quantities
493 rd <- pi / 180
494
495 # Generate leaf angle distribution from average leaf angle (ellipsoidal) or (a, b) parameters
496 if (typelidf == 1) {
497 foliar_distrib <- dladgen(lidfa, lidfb)
498 lidf <- foliar_distrib$lidf
499 litab <- foliar_distrib$litab
500
501 } else if (typelidf == 2) {
502 foliar_distrib <- campbell(lidfa)
503 lidf <- foliar_distrib$lidf
504 litab <- foliar_distrib$litab
505 }
506
507 if (lai < 0) {
508 message("Please define positive LAI value")
509 rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
510 alfast <- alfadt <- 0 * rsoil
511 } else if (lai == 0) {
512 tss <- too <- tsstoo <- tdd <- 1.0
513 rdd <- rsd <- tsd <- rdo <- tdo <- 0.0
514 rso <- rsos <- rsod <- rsodt <- 0.0
515 rddt <- rsdt <- rdot <- rsost <- rsot <- rsoil
516 alfast <- alfadt <- 0 * rsoil
517 } else if (lai > 0) {
518 cts <- cos(rd * tts)
519 cto <- cos(rd * tto)
520 ctscto <- cts * cto
521 tants <- tan(rd * tts)
522 tanto <- tan(rd * tto)
523 cospsi <- cos(rd * psi)
524 dso <- sqrt(tants * tants + tanto * tanto - 2.0 * tants * tanto * cospsi)
525
526 # Clumping effects
527 cs <- co <- 1.0
528 if (cv <= 1.0) {
529 cs <- 1.0 - (1.0 - cv)^(1.0 / cts)
530 co <- 1.0 - (1.0 - cv)^(1.0 / cto)
531 }
532 overlap <- 0.0
533 if (zeta > 0.0) {
534 overlap <- min(cs * (1.0 - co), co * (1.0 - cs)) * exp(-dso / zeta)
535 }
536 fcd <- cs * co + overlap
537 fcs <- (1.0 - cs) * co - overlap
538 fod <- cs * (1.0 - co) - overlap
539 fos <- (1.0 - cs) * (1.0 - co) + overlap
540 fcdc <- 1.0 - (1.0 - fcd)^(0.5 / cts + 0.5 / cto)
541
542 # Part depending on diss, fraction_brown, and leaf optical properties
543 # First save the input fraction_brown as the old fraction_brown, as the following change is only artificial
544 # Better define an fraction_brown that is actually used: fb, so that the input is not modified!
545
546 fb <- fraction_brown
547 # if only green leaves
548 if (fraction_brown == 0.0) {
549 fb <- 0.5
550 leafbrown$Reflectance <- leafgreen$Reflectance
551 leafbrown$Transmittance <- leafgreen$Transmittance
552 }
553 if (fraction_brown == 1.0) {
554 fb <- 0.5
555 leafgreen$Reflectance <- leafbrown$Reflectance
556 leafgreen$Transmittance <- leafbrown$Transmittance
557 }
558 s <- (1.0 - diss) * fb * (1.0 - fb)
559 # rho1 && tau1 : green foliage
560 # rho2 && tau2 : brown foliage (bottom layer)
561 rho1 <- ((1 - fb - s) * leafgreen$Reflectance + s * leafbrown$Reflectance) / (1 - fb)
562 tau1 <- ((1 - fb - s) * leafgreen$Transmittance + s * leafbrown$Transmittance) / (1 - fb)
563 rho2 <- (s * leafgreen$Reflectance + (fb - s) * leafbrown$Reflectance) / fb
564 tau2 <- (s * leafgreen$Transmittance + (fb - s) * leafbrown$Transmittance) / fb
565
566 # angular distance, compensation of shadow length
567 # Calculate geometric factors associated with extinction and scattering
568 # Initialise sums
569 ks <- ko <- bf <- sob <- sof <- 0
570
571 # Weighted sums over LIDF
572
573 for (i in 1:seq_along(litab)) {
574 ttl <- litab[i]
575 ctl <- cos(rd * ttl)
576 # SAIL volscatt function gives interception coefficients
577 # and two portions of the volume scattering phase function to be
578 # multiplied by rho and tau, respectively
579 resvolscatt <- volscatt(tts, tto, psi, ttl)
580 chi_s <- resvolscatt$chi_s
581 chi_o <- resvolscatt$chi_o
582 frho <- resvolscatt$frho
583 ftau <- resvolscatt$ftau
584 # Extinction coefficients
585 ksli <- chi_s / cts
586 koli <- chi_o / cto
587 # Area scattering coefficient fractions
588 sobli <- frho * pi / ctscto
589 sofli <- ftau * pi / ctscto
590 bfli <- ctl * ctl
591 ks <- ks + ksli * lidf[i]
592 ko <- ko + koli * lidf[i]
593 bf <- bf + bfli * lidf[i]
594 sob <- sob + sobli * lidf[i]
595 sof <- sof + sofli * lidf[i]
596 }
597 # Geometric factors to be used later in combination with rho and tau
598 sdb <- 0.5 * (ks + bf)
599 sdf <- 0.5 * (ks - bf)
600 dob <- 0.5 * (ko + bf)
601 dof <- 0.5 * (ko - bf)
602 ddb <- 0.5 * (1. + bf)
603 ddf <- 0.5 * (1. - bf)
604
605 # LAIs in two layers
606 lai1 <- (1 - fb) * lai
607 lai2 <- fb * lai
608
609 tss <- exp(-ks * lai)
610 ck <- exp(-ks * lai1)
611 alf <- 1e6
612 if (hot > 0.0) {
613 alf <- (dso / hot) * 2.0 / (ks + ko)
614 }
615 if (alf > 200.0) {
616 alf <- 200.0 # inserted H. Bach 1 / 3 / 04
617 }
618 if (alf == 0.0) {
619 # The pure hotspot
620 tsstoo <- tss
621 s1 <- (1 - ck) / (ks * lai)
622 s2 <- (ck - tss) / (ks * lai)
623 } else {
624 # Outside the hotspot
625 fhot <- lai * sqrt(ko * ks)
626 # Integrate 2 layers by exponential simpson method in 20 steps
627 # the steps are arranged according to equal partitioning
628 # of the derivative of the joint probability function
629 x1 <- y1 <- 0.0
630 f1 <- 1.0
631 ca <- exp(alf * (fb - 1.0))
632 fint <- (1.0 - ca) * .05
633 s1 <- 0.0
634 for (istep in 1:20) {
635 if (istep < 20) {
636 x2 <- -log(1. - istep * fint) / alf
637 } else {
638 x2 <- 1. - fb
639 }
640 y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
641 f2 <- exp(y2)
642 s1 <- s1 + (f2 - f1) * (x2 - x1) / (y2 - y1)
643 x1 <- x2
644 y1 <- y2
645 f1 <- f2
646 }
647 fint <- (ca - exp(-alf)) * .05
648 s2 <- 0.0
649 for (istep in 1:20) {
650 if (istep < 20) {
651 x2 <- -log(ca - istep * fint) / alf
652 } else {
653 x2 <- 1.0
654 }
655 y2 <- -(ko + ks) * lai * x2 + fhot * (1.0 - exp(-alf * x2)) / alf
656 f2 <- exp(y2)
657 s2 <- s2 + (f2 - f1) * (x2 - x1) / (y2 - y1)
658 x1 <- x2
659 y1 <- y2
660 f1 <- f2
661 }
662 tsstoo <- f1
663 }
664
665 # Calculate reflectances and transmittances
666 # Bottom layer
667 tss <- exp(-ks * lai2)
668 too <- exp(-ko * lai2)
669 sb <- sdb * rho2 + sdf * tau2
670 sf <- sdf * rho2 + sdb * tau2
671
672 vb <- dob * rho2 + dof * tau2
673 vf <- dof * rho2 + dob * tau2
674
675 w2 <- sob * rho2 + sof * tau2
676
677 sigb <- ddb * rho2 + ddf * tau2
678 sigf <- ddf * rho2 + ddb * tau2
679 att <- 1.0 - sigf
680 m2 <- (att + sigb) * (att - sigb)
681 m2[m2 < 0] <- 0
682 m <- sqrt(m2)
683 which_ncs <- which(m > 0.01)
684 which_cs <- which(m <= 0.01)
685
686 tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
687 rsod <- 0 * m
688 if (length(which_ncs) > 0) {
689 resncs <- nonconservativescattering(m[which_ncs], lai2, att[which_ncs], sigb[which_ncs],
690 ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
691 tdd[which_ncs] <- resncs$tdd
692 rdd[which_ncs] <- resncs$rdd
693 tsd[which_ncs] <- resncs$tsd
694 rsd[which_ncs] <- resncs$rsd
695 tdo[which_ncs] <- resncs$tdo
696 rdo[which_ncs] <- resncs$rdo
697 rsod[which_ncs] <- resncs$rsod
698 }
699 if (length(which_cs) > 0) {
700 rescs <- conservativescattering(m[which_cs], lai2, att[which_cs], sigb[which_cs],
701 ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
702 tdd[which_cs] <- rescs$tdd
703 rdd[which_cs] <- rescs$rdd
704 tsd[which_cs] <- rescs$tsd
705 rsd[which_cs] <- rescs$rsd
706 tdo[which_cs] <- rescs$tdo
707 rdo[which_cs] <- rescs$rdo
708 rsod[which_cs] <- rescs$rsod
709 }
710
711 # Set background properties equal to those of the bottom layer on a black soil
712 rddb <- rdd
713 rsdb <- rsd
714 rdob <- rdo
715 rsodb <- rsod
716 tddb <- tdd
717 tsdb <- tsd
718 tdob <- tdo
719 toob <- too
720 tssb <- tss
721 # Top layer
722 tss <- exp(-ks * lai1)
723 too <- exp(-ko * lai1)
724
725 sb <- sdb * rho1 + sdf * tau1
726 sf <- sdf * rho1 + sdb * tau1
727
728 vb <- dob * rho1 + dof * tau1
729 vf <- dof * rho1 + dob * tau1
730
731 w1 <- sob * rho1 + sof * tau1
732
733 sigb <- ddb * rho1 + ddf * tau1
734 sigf <- ddf * rho1 + ddb * tau1
735 att <- 1.0 - sigf
736
737 m2 <- (att + sigb) * (att - sigb)
738 m2[m2 < 0] <- 0
739 m <- sqrt(m2)
740 which_ncs <- which(m > 0.01)
741 which_cs <- which(m <= 0.01)
742
743 tdd <- rdd <- tsd <- rsd <- tdo <- rdo <- 0 * m
744 rsod <- 0 * m
745 if (length(which_ncs) > 0) {
746 resncs <- nonconservativescattering(m[which_ncs], lai1, att[which_ncs], sigb[which_ncs],
747 ks, ko, sf[which_ncs], sb[which_ncs], vf[which_ncs], vb[which_ncs], tss, too)
748 tdd[which_ncs] <- resncs$tdd
749 rdd[which_ncs] <- resncs$rdd
750 tsd[which_ncs] <- resncs$tsd
751 rsd[which_ncs] <- resncs$rsd
752 tdo[which_ncs] <- resncs$tdo
753 rdo[which_ncs] <- resncs$rdo
754 rsod[which_ncs] <- resncs$rsod
755 }
756 if (length(which_cs) > 0) {
757 rescs <- conservativescattering(m[which_cs], lai1, att[which_cs], sigb[which_cs],
758 ks, ko, sf[which_cs], sb[which_cs], vf[which_cs], vb[which_cs], tss, too)
759 tdd[which_cs] <- rescs$tdd
760 rdd[which_cs] <- rescs$rdd
761 tsd[which_cs] <- rescs$tsd
762 rsd[which_cs] <- rescs$rsd
763 tdo[which_cs] <- rescs$tdo
764 rdo[which_cs] <- rescs$rdo
765 rsod[which_cs] <- rescs$rsod
766 }
767
768 # combine with bottom layer reflectances and transmittances (adding method)
769 rn <- 1.0 - rdd * rddb
770 tup <- (tss * rsdb + tsd * rddb) / rn
771 tdn <- (tsd + tss * rsdb * rdd) / rn
772 rsdt <- rsd + tup * tdd
773 rdot <- rdo + tdd * (rddb * tdo + rdob * too) / rn
774 rsodt <- rsod + (tss * rsodb + tdn * rdob) * too + tup * tdo
775
776 rsost <- (w1 * s1 + w2 * s2) * lai
777
778 rsot <- rsost + rsodt
779
780 # Diffuse reflectances at the top and the bottom are now different
781 rddt_t <- rdd + tdd * rddb * tdd / rn
782 rddt_b <- rddb + tddb * rdd * tddb / rn
783
784 # Transmittances of the combined canopy layers
785 tsst <- tss * tssb
786 toot <- too * toob
787 tsdt <- tss * tsdb + tdn * tddb
788 tdot <- tdob * too + tddb * (tdo + rdd * rdob * too) / rn
789 tddt <- tdd * tddb / rn
790
791 # Apply clumping effects to vegetation layer
792 rddcb <- cv * rddt_b
793 rddct <- cv * rddt_t
794 tddc <- 1 - cv + cv * tddt
795 rsdc <- cs * rsdt
796 tsdc <- cs * tsdt
797 rdoc <- co * rdot
798 tdoc <- co * tdot
799 tssc <- 1 - cs + cs * tsst
800 tooc <- 1 - co + co * toot
801
802 # new weight function fcdc for crown contribution (W. Verhoef, 22-05-08)
803 rsoc <- fcdc * rsot
804 tssooc <- fcd * tsstoo + fcs * toot + fod * tsst + fos
805 # Canopy absorptance for black background (W. Verhoef, 02-03-04)
806 alfas <- 1. - tssc - tsdc - rsdc
807 alfad <- 1. - tddc - rddct
808 # Add the soil background
809 rn <- 1 - rddcb * rddsoil
810 tup <- (tssc * rsdsoil + tsdc * rddsoil) / rn
811 tdn <- (tsdc + tssc * rsdsoil * rddcb) / rn
812
813 rddt <- rddct + tddc * rddsoil * tddc / rn
814 rsdt <- rsdc + tup * tddc
815 rdot <- rdoc + tddc * (rddsoil * tdoc + rdosoil * tooc) / rn
816 rsot <- rsoc + tssooc * rsosoil + tdn * rdosoil * tooc + tup * tdoc
817
818 # Effect of soil background on canopy absorptances (W. Verhoef, 02-03-04)
819 alfast <- alfas + tup * alfad
820 alfadt <- alfad * (1. + tddc * rddsoil / rn)
821 }
822 my_list <- list("rdot" = rdot, "rsot" = rsot, "rddt" = rddt, "rsdt" = rsdt,
823 "alfast" = alfast, "alfadt" = alfadt)
824 return(my_list)
825 }
826
827
828
829 #" computes non conservative scattering conditions
830 #" @param m numeric.
831 #" @param lai numeric. Leaf Area Index
832 #" @param att numeric.
833 #" @param sigb numeric.
834 #" @param ks numeric.
835 #" @param ko numeric.
836 #" @param sf numeric.
837 #" @param sb numeric.
838 #" @param vf numeric.
839 #" @param vb numeric.
840 #" @param tss numeric.
841 #" @param too numeric.
842 #"
843 #" @return list. tdd, rdd, tsd, rsd, tdo, rdo, rsod
844 #"
845 #" @export
846 nonconservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
847
848 e1 <- exp(-m * lai)
849 e2 <- e1 * e1
850 rinf <- (att - m) / sigb
851 rinf2 <- rinf * rinf
852 re <- rinf * e1
853 denom <- 1. - rinf2 * e2
854
855 j1ks <- jfunc1(ks, m, lai)
856 j2ks <- jfunc2(ks, m, lai)
857 j1ko <- jfunc1(ko, m, lai)
858 j2ko <- jfunc2(ko, m, lai)
859
860 ps <- (sf + sb * rinf) * j1ks
861 qs <- (sf * rinf + sb) * j2ks
862 pv <- (vf + vb * rinf) * j1ko
863 qv <- (vf * rinf + vb) * j2ko
864
865 tdd <- (1. - rinf2) * e1 / denom
866 rdd <- rinf * (1. - e2) / denom
867 tsd <- (ps - re * qs) / denom
868 rsd <- (qs - re * ps) / denom
869 tdo <- (pv - re * qv) / denom
870 rdo <- (qv - re * pv) / denom
871
872 z <- jfunc2(ks, ko, lai)
873 g1 <- (z - j1ks * too) / (ko + m)
874 g2 <- (z - j1ko * tss) / (ks + m)
875
876 tv1 <- (vf * rinf + vb) * g1
877 tv2 <- (vf + vb * rinf) * g2
878
879 t1 <- tv1 * (sf + sb * rinf)
880 t2 <- tv2 * (sf * rinf + sb)
881 t3 <- (rdo * qs + tdo * ps) * rinf
882
883 # Multiple scattering contribution to bidirectional canopy reflectance
884 rsod <- (t1 + t2 - t3) / (1. - rinf2)
885 my_list <- list("tdd" = tdd, "rdd" = rdd, "tsd" = tsd,
886 "rsd" = rsd, "tdo" = tdo, "rdo" = rdo,
887 "rsod" = rsod)
888 return(my_list)
889 }
890
891 #" computes conservative scattering conditions
892 #" @param m numeric.
893 #" @param lai numeric. Leaf Area Index
894 #" @param att numeric.
895 #" @param sigb numeric.
896 #" @param ks numeric.
897 #" @param ko numeric.
898 #" @param sf numeric.
899 #" @param sb numeric.
900 #" @param vf numeric.
901 #" @param vb numeric.
902 #" @param tss numeric.
903 #" @param too numeric.
904 #"
905 #" @return list. tdd, rdd, tsd, rsd, tdo, rdo, rsod
906 #"
907 #" @export
908 conservativescattering <- function(m, lai, att, sigb, ks, ko, sf, sb, vf, vb, tss, too) {
909
910 # near or complete conservative scattering
911 j4 <- jfunc4(m, lai)
912 amsig <- att - sigb
913 apsig <- att + sigb
914 rtp <- (1 - amsig * j4) / (1 + amsig * j4)
915 rtm <- (-1 + apsig * j4) / (1 + apsig * j4)
916 rdd <- 0.5 * (rtp + rtm)
917 tdd <- 0.5 * (rtp - rtm)
918
919 dns <- ks * ks - m * m
920 dno <- ko * ko - m * m
921 cks <- (sb * (ks - att) - sf * sigb) / dns
922 cko <- (vb * (ko - att) - vf * sigb) / dno
923 dks <- (-sf * (ks + att) - sb * sigb) / dns
924 dko <- (-vf * (ko + att) - vb * sigb) / dno
925 ho <- (sf * cko + sb * dko) / (ko + ks)
926
927 rsd <- cks * (1 - tss * tdd) - dks * rdd
928 rdo <- cko * (1 - too * tdd) - dko * rdd
929 tsd <- dks * (tss - tdd) - cks * tss * rdd
930 tdo <- dko * (too - tdd) - cko * too * rdd
931 # Multiple scattering contribution to bidirectional canopy reflectance
932 rsod <- ho * (1 - tss * too) - cko * tsd * too - dko * rsd
933
934 my_list <- list("tdd" = tdd, "rdd" = rdd, "tsd" = tsd,
935 "rsd" = rsd, "tdo" = tdo, "rdo" = rdo,
936 "rsod" = rsod)
937 return(my_list)
938 }
939
940
941
942
943
944
945 #" computes the leaf angle distribution function value (freq)
946 #"
947 #" Ellipsoidal distribution function characterised by the average leaf
948 #" inclination angle in degree (ala)
949 #" Campbell 1986
950 #" @param ala average leaf angle
951 #" @return foliar_distrib list. lidf and litab
952 #" @export
953 campbell <- function(ala) {
954
955 tx1 <- c(10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88., 90.)
956 tx2 <- c(0., 10., 20., 30., 40., 50., 60., 70., 80., 82., 84., 86., 88.)
957 litab <- (tx2 + tx1) / 2
958 n <- length(litab)
959 tl1 <- tx1 * (pi / 180)
960 tl2 <- tx2 * (pi / 180)
961 excent <- exp(-1.6184e-5 * ala**3 + 2.1145e-3 * ala**2 - 1.2390e-1 * ala + 3.2491)
962 sum0 <- 0
963
964 freq <- c()
965 for (i in 1:n) {
966 x1 <- excent / (sqrt(1. + excent**2. * tan(tl1[i])**2))
967 x2 <- excent / (sqrt(1. + excent**2. * tan(tl2[i])**2))
968 if (excent == 1) {
969 freq[i] <- abs(cos(tl1[i]) - cos(tl2[i]))
970 } else {
971 alpha <- excent / sqrt(abs(1 - excent**2))
972 alpha2 <- alpha**2
973 x12 <- x1**2
974 x22 <- x2**2
975 alpx1 <- 0 * alpha2
976 alpx2 <- 0 * alpha2
977 almx1 <- 0 * alpha2
978 almx2 <- 0 * alpha2
979 if (excent > 1) {
980 alpx1 <- sqrt(alpha2[excent > 1] + x12[excent > 1])
981 alpx2[excent > 1] <- sqrt(alpha2[excent > 1] + x22[excent > 1])
982 dum <- x1 * alpx1 + alpha2 * log(x1 + alpx1)
983 freq[i] <- abs(dum - (x2 * alpx2 + alpha2 * log(x2 + alpx2)))
984 } else {
985 almx1 <- sqrt(alpha2 - x12)
986 almx2 <- sqrt(alpha2 - x22)
987 dum <- x1 * almx1 + alpha2 * asin(x1 / alpha)
988 freq[i] <- abs(dum - (x2 * almx2 + alpha2 * asin(x2 / alpha)))
989 }
990 }
991 }
992 sum0 <- sum(freq)
993 freq0 <- freq / sum0
994 foliar_distrib <- list("lidf" = freq0, "litab" = litab)
995 return(foliar_distrib)
996 }
997
998 #" computes the leaf angle distribution function value (freq)
999 #"
1000 #" Using the original bimodal distribution function initially proposed in SAIL
1001 #" references
1002 #" ----------
1003 #" (Verhoef1998) Verhoef, Wout. Theory of radiative transfer models applied
1004 #" in optical remote sensing of vegetation canopies.
1005 #" nationaal Lucht en Ruimtevaartlaboratorium, 1998.
1006 #" http: / / library.wur.nl / WebQuery / clc / 945481.
1007 #" @param a controls the average leaf slope
1008 #" @param b controls the distribution"s bimodality
1009 #" LIDF type a b
1010 #" Planophile 1 0
1011 #" Erectophile -1 0
1012 #" Plagiophile 0 -1
1013 #" Extremophile 0 1
1014 #" Spherical -0.35 -0.15
1015 #" Uniform 0 0
1016 #" requirement: ||lidfa|| + ||lidfb|| < 1
1017 #"
1018 #" @return foliar_distrib list. lidf and litab
1019 #" @export
1020 dladgen <- function(a, b) {
1021 litab <- c(5., 15., 25., 35., 45., 55., 65., 75., 81., 83., 85., 87., 89.)
1022 freq <- c()
1023 for (i1 in 1:8) {
1024 t <- i1 * 10
1025 freq[i1] <- dcum(a, b, t)
1026 }
1027 for (i2 in 9:12) {
1028 t <- 80. + (i2 - 8) * 2.
1029 freq[i2] <- dcum(a, b, t)
1030 }
1031 freq[13] <- 1
1032 for (i in 13:2) {
1033 freq[i] <- freq[i] - freq[i - 1]
1034 }
1035 foliar_distrib <- list("lidf" = freq, "litab" = litab)
1036 return(foliar_distrib)
1037 }
1038
1039 #" dcum function
1040 #" @param a numeric. controls the average leaf slope
1041 #" @param b numeric. controls the distribution"s bimodality
1042 #" @param t numeric. angle
1043 #" @return f
1044 #" @export
1045 dcum <- function(a, b, t) {
1046 rd <- pi / 180
1047 if (a >= 1) {
1048 f <- 1 - cos(rd * t)
1049 } else {
1050 eps <- 1e-8
1051 delx <- 1
1052 x <- 2 * rd * t
1053 p <- x
1054 while (delx >= eps) {
1055 y <- a * sin(x) + .5 * b * sin(2. * x)
1056 dx <- .5 * (y - x + p)
1057 x <- x + dx
1058 delx <- abs(dx)
1059 }
1060 f <- (2. * y + p) / pi
1061 }
1062 return(f)
1063 }
1064
1065 #" J1 function with avoidance of singularity problem
1066 #"
1067 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
1068 #" @param l numeric.
1069 #" @param t numeric. Leaf Area Index
1070 #" @return jout numeric.
1071 #" @export
1072 jfunc1 <- function(k, l, t) {
1073 # J1 function with avoidance of singularity problem
1074 del <- (k - l) * t
1075 jout <- 0 * l
1076 jout[which(abs(del) > 1e-3)] <- (exp(-l[which(abs(del) > 1e-3)] * t) - exp(-k * t)) / (k - l[which(abs(del) > 1e-3)])
1077 jout[which(abs(del) <= 1e-3)] <- 0.5 * t * (exp(-k * t) + exp(-l[which(abs(del) <= 1e-3)] * t)) * (1 - del[which(abs(del) <= 1e-3)] * del[which(abs(del) <= 1e-3)] / 12)
1078 return(jout)
1079 }
1080
1081 #" J2 function with avoidance of singularity problem
1082 #"
1083 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
1084 #" @param l numeric.
1085 #" @param t numeric. Leaf Area Index
1086 #" @return jout numeric.
1087 #" @export
1088 jfunc2 <- function(k, l, t) {
1089 # J2 function
1090 jout <- (1. - exp(-(k + l) * t)) / (k + l)
1091 return(jout)
1092 }
1093
1094 #" J3 function with avoidance of singularity problem
1095 #"
1096 #" @param k numeric. Extinction coefficient for direct (solar or observer) flux
1097 #" @param l numeric.
1098 #" @param t numeric. Leaf Area Index
1099 #" @return jout numeric.
1100 #" @export
1101 jfunc3 <- function(k, l, t) {
1102 out <- (1. - exp(-(k + l) * t)) / (k + l)
1103 return(out)
1104 }
1105
1106
1107 #" j4 function for treating (near) conservative scattering
1108 #"
1109 #" @param m numeric. Extinction coefficient for direct (solar or observer) flux
1110 #" @param t numeric. Leaf Area Index
1111 #" @return jout numeric.
1112 #" @export
1113 jfunc4 <- function(m, t) {
1114
1115 del <- m * t
1116 out <- 0 * del
1117 out[del > 1e-3] <- (1 - exp(-del)) / (m * (1 + exp(-del)))
1118 out[del <= 1e-3] <- 0.5 * t * (1. - del * del / 12.)
1119 return(out)
1120 }
1121
1122
1123 #" compute volume scattering functions and interception coefficients
1124 #" for given solar zenith, viewing zenith, azimuth and leaf inclination angle.
1125 #"
1126 #" @param tts numeric. solar zenith
1127 #" @param tto numeric. viewing zenith
1128 #" @param psi numeric. azimuth
1129 #" @param ttl numeric. leaf inclination angle
1130 #" @return res list. includes chi_s, chi_o, frho, ftau
1131 #" @export
1132 volscatt <- function(tts, tto, psi, ttl) {
1133 #********************************************************************************
1134 #* chi_s = interception functions
1135 #* chi_o = interception functions
1136 #* frho = function to be multiplied by leaf reflectance rho
1137 #* ftau = functions to be multiplied by leaf transmittance tau
1138 #********************************************************************************
1139 # Wout Verhoef, april 2001, for CROMA
1140
1141 rd <- pi / 180
1142 costs <- cos(rd * tts)
1143 costo <- cos(rd * tto)
1144 sints <- sin(rd * tts)
1145 sinto <- sin(rd * tto)
1146 cospsi <- cos(rd * psi)
1147 psir <- rd * psi
1148 costl <- cos(rd * ttl)
1149 sintl <- sin(rd * ttl)
1150 cs <- costl * costs
1151 co <- costl * costo
1152 ss <- sintl * sints
1153 so <- sintl * sinto
1154
1155 #c ..............................................................................
1156 #c betas -bts- and betao -bto- computation
1157 #c Transition angles (beta) for solar (betas) and view (betao) directions
1158 #c if thetav + thetal > pi / 2, bottom side of the leaves is observed for leaf azimut
1159 #c interval betao + phi<leaf azimut<2pi-betao + phi.
1160 #c if thetav + thetal<pi / 2, top side of the leaves is always observed, betao=pi
1161 #c same consideration for solar direction to compute betas
1162 #c ..............................................................................
1163
1164 cosbts <- 5
1165 if (abs(ss) > 1e-6) {
1166 cosbts <- -cs / ss
1167 }
1168 cosbto <- 5
1169 if (abs(so) > 1e-6) {
1170 cosbto <- -co / so
1171 }
1172
1173 if (abs(cosbts) < 1) {
1174 bts <- acos(cosbts)
1175 ds <- ss
1176 } else {
1177 bts <- pi
1178 ds <- cs
1179 }
1180 chi_s <- 2. / pi * ((bts - pi * .5) * cs + sin(bts) * ss)
1181 if (abs(cosbto) < 1) {
1182 bto <- acos(cosbto)
1183 doo <- so
1184 } else if (tto < 90) {
1185 bto <- pi
1186 doo <- co
1187 } else {
1188 bto <- 0
1189 doo <- -co
1190 }
1191 chi_o <- 2. / pi * ((bto - pi * .5) * co + sin(bto) * so)
1192
1193 #c ..............................................................................
1194 #c computation of auxiliary azimut angles bt1, bt2, bt3 used
1195 #c for the computation of the bidirectional scattering coefficient w
1196 #c .............................................................................
1197
1198 btran1 <- abs(bts - bto)
1199 btran2 <- pi - abs(bts + bto - pi)
1200
1201 if (psir <= btran1) {
1202 bt1 <- psir
1203 bt2 <- btran1
1204 bt3 <- btran2
1205 } else {
1206 bt1 <- btran1
1207 if (psir <= btran2) {
1208 bt2 <- psir
1209 bt3 <- btran2
1210 } else {
1211 bt2 <- btran2
1212 bt3 <- psir
1213 }
1214 }
1215 t1 <- 2. * cs * co + ss * so * cospsi
1216 t2 <- 0
1217 if (bt2 > 0) {
1218 t2 <- sin(bt2) * (2. * ds * doo + ss * so * cos(bt1) * cos(bt3))
1219 }
1220
1221 denom <- 2. * pi * pi
1222 frho <- ((pi - bt2) * t1 + t2) / denom
1223 ftau <- (-bt2 * t1 + t2) / denom
1224
1225 if (frho < 0) {
1226 frho <- 0
1227 }
1228 if (ftau < 0) {
1229 ftau <- 0
1230 }
1231 res <- list("chi_s" = chi_s, "chi_o" = chi_o, "frho" = frho, "ftau" = ftau)
1232 return(res)
1233 }