#' Depth of the given soil volume
#'
#' Given a shape and mass for the roots, calculate the depth of a specific non-root volume. A specific solution is implemented for linear root volumes, and a general solution otherwise requires the \code{massLiveRoots.fn} to be specified. The general solution is a binary search algorithm.
#'
#' @param totalRootMassPerArea an integer that is the total mass per area of the roots, generally in g cm-3
#' @param rootDepthMax an integer that is the maximum root depth, generally in cm
#' @param rootDensity a numeric that is the root density in g per cm3
#' @param soilLength a numeric of the unit length of soil volume, generally 1 cm
#' @param soilWidth a numeric of the unit width of the soil column, generally 1cm
#' @param shape a character of the shape of the root distribution, currently only 'linear' is implemented
#' @param ...
#' @param nonRootVolume.arr a numerical array of the cumulative non-rooting volume, assumed to be monotonically increasing
#' @param massLiveRoots.fn a function defining the mass of the live roots, if shape is not linear then this must be defined.
#' @param relTol a numeric the relative tolerance accepted for the general depth solution
#' @param verbose a boolean flag, true prints out a lot of strings to help debug the general solution
#'
#' @return A numeric corresponding to the depth of the specified non root volume
#'
#' @export
calculateDepthOfNonRootVolume <- function(nonRootVolume.arr,
massLiveRoots.fn = NULL,
totalRootMassPerArea,
rootDepthMax,
rootDensity,
shape = 'linear',
soilLength = 1, soilWidth = 1,
relTol = 1e-4,
verbose = FALSE,
...){
####
totalRootMass <- soilLength*soilWidth*totalRootMassPerArea
if(verbose) print(paste('totalRootMass = ', totalRootMass))
totalRootVolume <- totalRootMass/rootDensity
if(verbose) print(paste('totalRootVolume = ', totalRootVolume))
if(totalRootVolume > soilLength*soilWidth*rootDepthMax){
stop('Bad root volume')
}
if(totalRootMassPerArea == 0){
return(nonRootVolume.arr*soilLength*soilWidth)
}
##If the shape is linear then solve the non-root volume depth explicitly
if(shape == 'linear'){
if(verbose) print('We are linear.')
rootWidth <- totalRootVolume*2/(rootDepthMax*soilLength)
#nonRootVolume = ((rootWidth/rootDepthMax*depth^2)/2 + depth*(soilWidth-rootWidth))*soilLength
# 0 = rootWidth / (2*rootDepthMax) * depth ^2 +
# (soilWidth-rootWidth) * depth - nonRootVolume/soilLength ##solve for depth
coef1 <- rootWidth / (2*rootDepthMax)
coef2 <- soilWidth-rootWidth
coef3 <- -nonRootVolume.arr/soilLength
ansDepth <- (-coef2 + sqrt(coef2^2-4*coef1*coef3))/(2*coef1)
#correct for beyond root zone
behondRootZone <- nonRootVolume.arr > soilLength*soilWidth*rootDepthMax - totalRootVolume
if(verbose) print(paste0('nonRootVolume.arr: [', paste0(nonRootVolume.arr, collapse = ', '), ']'))
if(verbose) print(paste0('nrv comparison: [', paste0(soilLength*soilWidth*rootDepthMax - totalRootVolume, collapse = ', '), ']'))
if(verbose) print(paste0('beyondRootZon: [', paste0(behondRootZone, collapse = ', '), ']'))
ansDepth[behondRootZone] <- (rootDepthMax +
(nonRootVolume.arr - soilLength*soilWidth*rootDepthMax + totalRootVolume ) /
(soilLength*soilWidth)) [behondRootZone] #treat as a square
return(ansDepth)
}else{ ##Otherwise apply a general binary search algorithm
#What is the non rooting volumne down to rooting max
nonRootVolumeToRootMax <- (soilLength*soilWidth*rootDepthMax) - totalRootVolume
previousDepth <- 0
possibleDepth.arr <- rep(NA, length=length(nonRootVolume.arr))
for(nonRootVolume.index in (1:length(nonRootVolume.arr))){
if(verbose)cat(paste('layerIndex: ', nonRootVolume.index, '----\n'))
nonRootVolume <- nonRootVolume.arr[nonRootVolume.index]
if(verbose)cat(paste('\ttarget nonRootVolume=', nonRootVolume, '\n'))
if(nonRootVolume < relTol){
if(verbose)cat('no volume\n')
possibleDepth.arr[nonRootVolume] <- previousDepth
next
}
#If we overfill that root zone volume
if(nonRootVolumeToRootMax <= nonRootVolume){
if(verbose)cat('outside root zone\n')
#Subtract the root zone volumne, find the depth beyond and add to the max rooting depth
possibleDepth.arr[nonRootVolume.index] <- (nonRootVolume-nonRootVolumeToRootMax)/(soilLength*soilWidth) + rootDepthMax
previousDepth <- possibleDepth.arr[nonRootVolume.index]
next
}
#incrament <- (min(rootDepthMax, nonRootVolume/(soilLength*soilWidth)) - previousDepth)
#possibleDepth <- min(rootDepthMax, nonRootVolume/(soilLength*soilWidth)) - incrament / 2
incrament <- rootDepthMax - previousDepth
possibleDepth <- previousDepth + incrament/2
maxSearchDepth <- ceiling(log(relTol, base=0.5)) #search down to relTol*depth
for(ii in 1:maxSearchDepth){
if(verbose)cat(paste('ii=', ii, '\n\tincrament=', incrament, '\n\tpossibleDepth=', possibleDepth))
rootVolume <- massLiveRoots.fn(layerTop = 0, layerBottom = possibleDepth,
totalRootMassPerArea = totalRootMassPerArea,
rootDepthMax=rootDepthMax,
soilLength=soilLength, soilWidth=soilWidth, shape=shape
) / rootDensity
if(verbose)cat(paste('\n\trootVolume_test=', rootVolume))
nonRootVolume_test <- (possibleDepth * soilLength * soilWidth) - rootVolume
if(verbose)cat(paste('\n\tnonRootVolume_test=', nonRootVolume_test))
if( (abs(nonRootVolume-nonRootVolume_test) / nonRootVolume) < relTol){
if(verbose)cat('\n\tdone\n')
break
}else{
incrament <- incrament / 2
if(verbose)cat(paste('\n\trelTol=', abs(nonRootVolume-nonRootVolume_test) / nonRootVolume))
#Should we go up?
if(nonRootVolume > nonRootVolume_test){
if(verbose)cat('\n\tup\n')
possibleDepth <- possibleDepth + incrament / 2
}else{
if(verbose)cat('\n\tdown\n')
possibleDepth <- possibleDepth - incrament / 2
}#if-else up/down
}#if-else relTol
}#for-loop refining search
if(ii == maxSearchDepth){
if(verbose) warning(paste('Target volumne precision [', relTol, '] is not reached'))
}
previousDepth <- possibleDepth
possibleDepth.arr[nonRootVolume.index] <- possibleDepth
}#for-loop across array
return(possibleDepth.arr)
}#if checking for unknown shape
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.