R/martincurveb.R

In kevans27/aqeco: Basic Functions for Aquatic Ecologists

martincurveb = function(des, flu, rds = NULL, tds = NULL){

  #Function to solve for exponent of logarithmic fit
  solveit <- function(fz, fz0, z, z0){
    b <- log10(fz0/fz)/log10(z/z0)
    return(b)
  }

  bs <- c()

  #Check if the user defined test and reference depths (tds & rds);
  #if not, use all depths for all variables
  if(is.null(tds)){
    todes <- des
    matchflu <- flu
    refflu <- flu
    reldes <- des
  }else{todes <- des[which(des %in% tds)]
        matchflu <- flu[which(des %in% tds)]
        refflu <- flu[which(des %in% rds)]
        reldes <- des[which(des %in% rds)]
  }


  #Now match everything on the list of target depths to their fluxes
  for(i in 1:length(todes)){
    fz = matchflu[i]
    z = todes[i]
    #Find the indexes of the reference depths which are deeper than the target depth
    inds <- which(reldes > z)
    #If no deeper depths than target depth, move on
    if(length(inds) == 0){
      break
    }
    #Calculate b for each reference depth for your target depth
    for(j in inds){
      z0 = reldes[j]
      fz0 = refflu[j]
      b = solveit(fz, fz0, z, z0)
      bs <- c(bs, b)
    }
  }

  #Find out some stuff about your b values and return those as a list
  lowest <- min(bs)
  highest <- max(bs)
  avest <- mean(bs)
  middl <- median(bs)
  returnvals <- list("Minimum" = lowest, "Maximum" = highest, "Mean/Average" = avest, "Median" = middl, "AllVals" = bs)
  return(returnvals)
}