temperheic:

Documented in thBoundary thUnits

#'Temperheic Aquifer, Boundary, and Units objects
#'
#'@description The \code{temperheic} package relies on three different S3
#'  objects: \code{thAquifer}, \code{thBoundary}, and \code{thUnits}.  The "th"
#'  prefix is short for "temperheic".  A \code{thAquifer} object contains
#'  user-specified and derived aquifer properties.  A \code{thBoundary} object
#'  contains user-specified and derived values that describe a time varying
#'  temperature boudary as a cosine function.  A \code{thUnits} object describes
#'  the units of measure used to specify \code{thAquifer} and \code{thBoundary}
#'  objects.
#'
#'  NOTE THAT THE USER IS RESPONSIBLE FOR ENSURING THAT THE UNITS OF ALL
#'  ARGUMENTS PASSED TO \code{TEMPERHEIC} PACKAGE FUNCTIONS ARE CONSISTENT. See
#'  Details (below) for more information on units.
#'
#'@details \code{thAquifer} objects include values for the arguments provided by
#'  the user (see 'Arguments', above), along with the following derived values:
#'
#'  \itemize{
#'
#'  \item{\code{darcyFlux}}{  Darcy Flux (L t-1) is flux rate of water through
#'  the aquifer along the hydrologic gradient. Calculated as:
#'  \code{hydCond*headGrad}.}
#'
#'  \item{\code{density_bulk}}{  Density of the bulk medium (M L-3) is the
#'  average density of water and sediment, weighted by porosity.  Calculated as:
#'  \code{density_sed*(1-porosity)+density_h2o*porosity}.}
#'
#'  \item{\code{spHeat_bulk}}{  Specific heat of the bulk medium (E M-1 T-1) is
#'  the average specific heat of water and sediments, weighted by porosity.
#'  Calculated as: \code{spHeat_sed*(1-porosity)+spHeat_h2o*porosity}.}
#'
#'  \item{\code{velocity_h2o}}{  Water velocity (L T-1) is the aerial averaged
#'  rate of water movement along the hydrologic gradient in the aquifer.
#'  Calculated as: \code{darcyFlux/porosity}.}
#'
#'  \item{\code{darcyFlux_heat}}{  Heat velocity (L T-1) is the aerially
#'  averaged rate of heat movement along the hydraulic gradient in the aquifer.
#'  Calculated as: \code{darcyFlux_h2o*(volHeatCap_h2o/volHeatCap_bulk)} Note
#'  that the temperheic package represents 1D movement of water and heat, so
#'  that heat gradients and hydrologic gradients are always aligned.}
#'
#'  \item{\code{volHeatCap_h2o, volHeatCap_sed}}{  The volumetric heat capacity
#'  of of water, and of sediment [E L-3 T-1] are calcuated from the
#'  corresponding values for specific heat and density:
#'  \code{spHeat_<x>*density_<x>}.}
#'
#'  \item{\code{volHeatCap_bulk} {  The volumentric heat capacity of the bulk
#'  medium [E L-3 T-1] is calcuated as the porosity-weighted average of
#'  volumetric heat capacity of water and sediment:
#'  volHeatCap_sed*(1-porosity)+volHeatCap_h2o*porosity}.}
#'
#'  }
#'
#'@param hydCond Hydraulic conductivity (L t-1) (either vertical or horizontal)
#'  of the aquifer.
#'@param dispersivity Dispersivity (L) of the aquifer, appropriate for the the
#'  scale of the flow path observed (see Gelher 20XX). ## needs to be updated according to Bons et al., 2015
#'@param headGrad Average head gradient (L L-1) in the same dimension (vertical
#'  or horizonal) for the aquifer.
#'@param porosity Porosity (L3 L-3) of the aquifer; ratio of water volume to
#'  unit aquifer volume under saturated conditions.
#'@param thermCond_sed Thermal conductivity (E t-1 L-1 T-1) of the sediment.
#'@param thermCond_h2o Thermal conductivity (E t-1 L-1 T-1) of water.
#'@param spHeat_sed Specific heat (E M-1 T-1) of the sediment.
#'@param spHeat_h2o Specific heat (E M-1 T-1) of water.
#'@param density_sed Density (M L-3) of sediment.
#'@param density_h2o Density (M L-3) of water.
#'@param specificUnits A \code{thUnits} object generated by calling
#'  \code{thUnits()}. Specifies the units of length (L), mass (M), time (t),
#'  temperature (T), and energy (E) used in the numeric arguments.  A
#'  \code{thUnits} object is used only to label the values provided by the user,
#'  in hopes of helping the user avoid errors.  The \code{temperheic} package
#'  makes no conversions.  The user must ensure that all arguments passed to
#'  \code{thAquifer()} and \code{thBoundary()} have units consistent with
#'  \code{specificUnits}.
#'@return \code{thAquifer()} returns a S3 object of class \code{thAquifer} which
#'  includes all specified and derived parameters for the aquifer
#'@export
thAquifer = function (porosity, thermCond_sed, thermCond_h2o, spHeat_sed, spHeat_h2o, density_sed, density_h2o, specificUnits = thUnits()) {
#  spHeat_bulk = spHeat_sed * (1 - porosity) + spHeat_h2o * porosity # E M-1 T-1
#  density_bulk = density_sed * (1 - porosity) + density_h2o * porosity # E M-2 T-1

  volHeatCap_h2o = spHeat_h2o * density_h2o # E L-3 T-1
  volHeatCap_sed = spHeat_sed * density_sed # E L-3 T-1
  volHeatCap_bulk = volHeatCap_sed * (1 - porosity) + volHeatCap_h2o * porosity # E L-3 T-1

  thermCond_bulk = thermCond_sed * (1 - porosity) + thermCond_h2o * porosity

  newAquifer = .temperheic(
    thEnvir = environment(),
    thClass = "thAquifer",
    generalUnits =  c(
      rep("M L-3", 3),
      "L3 L-3",
      rep("E M-1 T-1", 3),
      rep("E t-1 L-1 T-1", 3),
      rep("E L-3 T-1", 3)
    ),
    specificUnits = specificUnits
  )
  return(newAquifer)
}

#' @rdname thAquifer
#' @details \code{thBoundary} objects include values for the arguments provided
#'   by the user (see 'Arguments', above), along with the following derived
#'   value:
#'
#'   \itemize{
#'
#'   \item{\code{frequency}}{  Frequency (t-1) of the temperature signal at the
#'   boundary.  Calculated as: \code{1/period}.}
#'
#'   }
#'
#' @return \code{thBoundary} returns an S3 object of class \code{thBoundary}, which
#'   includes the specified and derived parameters describing a temperature
#'   signal at the aquifer boundary.
#' @export
thBoundary = function(mean, amplitude, phase, period, specificUnits = thUnits()) {
  frequency = 1/period
  #frequency = (2*pi)/period
  newBoundary = .temperheic(
    thEnvir = environment(),
    thClass = "thBoundary",
    generalUnits = c(
      "T",
      "t-1",
      "T",
      rep("t",2)
    ),
    specificUnits = specificUnits
  )
  class(newBoundary) = c("thSpecifiedBoundary", class(newBoundary))
  return(newBoundary)
}

# CANT FIND ANY REMAINING REFERENCE TO THE NEXT FUNCTION.  fitCosine no longer
# matches signature in this function.  Commenting it out.

# Fit a cos wave to the specified column in the series.  Probably best to use lm
# rather than nls to avoid issues of initial guesses.  Byron has code for this.
# NOTE: an observedBoundary object cannot be calculated if we use our "regress
# one signal against the other" method because this method only calculate
# ampRatio and phaseLag (not amplitude and phase).  So, interestingly, for a
# thObservedSeries, the mean, amplitude, and phase of the "boundary" element
# will be NULL.
#
#.thObservedBoundary = function(observations, period) {
#  return(fitCosine(obserations, period))
# might have to set class or other attributes before returning results of
# fitCosine
#}

#' @rdname thAquifer
#' @details \code{thUnits} objects include unit labels for the parameter values
#'   stored in \code{thAquifer} and \code{thBoundary} objects.  Labels for units
#'   of length (L), mass (M), time (t), temperature (T), and Energy (E) are
#'   specified when creating a \code{thUnits} object.  The \code{thUnits} object
#'   can be stored in a variable and subsequently passed to \code{thAquifer()}
#'   and \code{thBoundary()}.  Note that the \code{thUnits} object is used ONLY
#'   to create lables for \code{thAquifer} and \code{thBoundary} values.  No
#'   attempt is made to convert to common units.  The user is responsible for
#'   ensuring that the units of all values in \code{thAquifer} and
#'   \code{thBoundary} objects agree with the unit labels displayed when the
#'   object is printed.
#' @return \code{thUnits()} return an S3 object of class \code{thUnits}, which
#'   includes specific units of length, mass, time, temperature, and energy. The
#'   values in a \code{thUnits} object are used to create unit labels for
#'   \code{thAquifer} and \code{thBoundary} parameter values.
#' @param L A chararcter string label for units of length.  Default is 'm' for
#'   meters.
#' @param M A chararcter string label for units of mass.  Default is 'kg' for
#'   kilograms.
#' @param t A chararcter string label for units of time.  Default is 's' for
#'   seconds.
#' @param T A chararcter string label for units of temperature.  Default is
#'   "degC" for degrees celsius.
#' @param E A chararcter string label for units of Energy. Default is "kJ" for
#'   kilojoules.
#' @export
thUnits = function(L = "m", M = "kg", t = "s", T = "degC", E = "kJ") {
  #get the list of argument names that can be passed to this function.
  unitNames = names(formals())
  #get the corresponding values, as a list
  unitList = sapply(unitNames, get, envir = environment(), simplify = F
  )

  newUnits = structure(
    unitList,
    class = c("thUnits"),
    longUnitName = c("Length", "Mass", "Time", "Temperature", "Energy")
  )
  return(newUnits)
}

#' @export
thHydro = function(hydCond, dispersivity = 0.001, headGrad, aquifer, specificUnits = thUnits()) {
  if(!is.thUnits(specificUnits)) stop("Units must be specified as a 'thUnits' object; call 'thUnits()' to generate such an object.")
  if(!is.thAquifer(aquifer)) stop("The 'aquifer' argument must be a thAquifer object; call 'thAquifer()' to generate such an object.")
  if(!identical(specificUnits, attr(aquifer, "specificUnits"))) stop("The units of the aquifer are not the same as those specificed in 'specificUnits.'")

  #pull these calculations into functions for reuse in thObservedHydro
  darcyFlux = hydCond * headGrad # Darcy velocity L t-1
  velocity_h2o = (darcyFlux / aquifer$porosity); # Water velocity L t-1
  #advectiveThermVel = darcyFlux * ((aquifer$density_h2o * aquifer$spHeat_h2o) / (aquifer$density_bulk * aquifer$spHeat_bulk))
  advectiveThermVel = darcyFlux * (aquifer$volHeatCap_h2o / aquifer$volHeatCap_bulk)
  #diffusivity_cond = aquifer$thermCond_bulk / (aquifer$density_bulk * aquifer$spHeat_bulk)
  diffusivity_cond = aquifer$thermCond_bulk / aquifer$volHeatCap_bulk
  ### added the Rau exponentto explore Rau's 2012 solution which has vt squared
  ### removed exponent
  diffusivity_disp = dispersivity * advectiveThermVel #* ((aquifer$density_h2o * aquifer$spHeat_h2o) / (aquifer$density_bulk * aquifer$spHeat_bulk))
  diffusivity_effective = diffusivity_disp + diffusivity_cond

  newHydro = .temperheic(
    thEnvir = environment(),
    thClass = "thHydro",
    generalUnits = c(
      rep("L t-1",2),
      rep("L2 t-1",3),
      "L",
      "L L-1",
      rep("L t-1",2),
      "L t-1"
    ),
    specificUnits = specificUnits
  )
  class(newHydro) = c("thSpecifiedHydro", class(newHydro))
  return(newHydro)

}

.thObservedHydro = function(series, headGrad, aquifer) {
  return("We need to implement this.")
}

#' @export
thSignal = function(hydro, boundary) {
  if((!identical(attr(hydro, "specificUnits"), attr(boundary, "specificUnits")))) stop("specificUnits attribute of the thAquifer, thHydro, and thBoundary must be identical.")
  # pecletNumber = (hydro$velocity_h2o * hydro$d50) / hydro$diffusivity_cond
  # diffusivity_disp = hydro$dispersivity * pecletNumber^hydro$exponent
  # diffusivity_effective = hydro$diffusivity_cond + diffusivity_disp

    # after Vogt et al., 2014; i.e. c
  if(hydro$darcyFlux == 0){
    phaseVel = sqrt(2 * hydro$diffusivity_cond * 2 * pi * boundary$frequency)
    thermDecayDist = sqrt(2 * hydro$diffusivity_cond / (2 * pi * boundary$frequency))
  } else{
    # num1 <- 64*(pi^2)*(boundary$frequency^2)*(hydro$diffusivity_effective^2)
    # rad1 <- (1 + (num1/(hydro$advectiveThermVel^4)))^0.5
    # rad2 <- (0.5 + (0.5*rad1))^0.5
    # phaseVel <- hydro$advectiveThermVel * rad2
    # thermDecayDist <- (2*hydro$diffusivity_effective)/(phaseVel-hydro$advectiveThermVel)

    ### updated formulas below based on PDE solutions of Logan, Luce, and Vogt
    rad1A = (hydro$advectiveThermVel^4 + (4 * 2 * pi * boundary$frequency * hydro$diffusivity_effective)^2)^0.5
    rad22A = sqrt(2)/(((rad1A + hydro$advectiveThermVel^2))^0.5 - sqrt(2)*hydro$advectiveThermVel)
    thermDecayDist = 2*hydro$diffusivity_effective * rad22A

    rad1B = (hydro$advectiveThermVel^4 + (4 * 2 * pi * boundary$frequency * hydro$diffusivity_effective)^2)^0.5
    rad22B = sqrt(2)/abs((rad1B - hydro$advectiveThermVel^2))^0.5
    phaseVel = 2*hydro$diffusivity_effective * rad22B  * 2 * pi * boundary$frequency
  }
  pecletNumberDisp = (hydro$advectiveThermVel * thermDecayDist) / hydro$diffusivity_disp
  pecletNumberCond = (hydro$advectiveThermVel * thermDecayDist) / hydro$diffusivity_cond
  pecletNumberEff = (hydro$advectiveThermVel * thermDecayDist) / hydro$diffusivity_effective

  pecletNumberPhaseDisp = (phaseVel * thermDecayDist) / hydro$diffusivity_disp
  pecletNumberPhaseCond = (phaseVel * thermDecayDist) / hydro$diffusivity_cond
  pecletNumberPhaseEff = (phaseVel * thermDecayDist) / hydro$diffusivity_effective

  dispersionDiffusionRatioCond = hydro$diffusivity_disp / hydro$diffusivity_cond
  dispersionDiffusionRatioEff = hydro$diffusivity_disp / hydro$diffusivity_effective
  specificUnits = attr(hydro$aquifer, "specificUnits")

  newSignal = .temperheic(
    thEnvir = environment(),
    thClass = "thSignal",
    generalUnits = c(
      "L t-1 L t2 L-2",
      "L2 t-2 t2 L-2",
      rep("L t-1",3),
      rep("L" ,3)
    ),
    specificUnits = specificUnits
  )
  return(newSignal)

}

.thObservedSignal = function(observedHydro, observedBoundary) {
  return("Need to implement this.")
  # call thSignal and then simply add the class "thObservedSignal" before returning.
}

#' creates a thObject from the variables in the environment passed to the
#' function.
#' @param thEnvir The environment to be searched for variables
#' @param thClass The class of the object to be create; also the name of the
#'   constructor function for the class
#' @param generalUnits A vector of general units (that uses the names from a
#'   specificUnits object) to describe the units of for each numeric value.
#' @param specificUnits A specificUnits object telling what specificUnits the
#'   values are in.
#' @param attr A character vector with the names of the variables in thEnvir
#'   that should be attached as attributes rather than being elements of the
#'   thObject
#' @param attrNames The attribute names to be used for the variables specifiec
#'   by attr.
.temperheic = function(thEnvir, thClass, generalUnits, specificUnits, attrs = NULL, attrNames = attrs) {
  attrNames = c(attrNames, "specificUnits", "thObjectNames")
  attrs = c(attrs, as.character(as.list(match.call())$specificUnits), "thObjects")
  if(!is.thUnits(specificUnits)) stop("Units must be specified as a 'thUnits' object; call 'thUnits()' to generate such an object.")
  elementNames = ls(envir = thEnvir)
  elementNames = elementNames[!(elementNames %in% attrs)]
  newTemperheic = structure(
    sapply(elementNames, function(x) get(x, envir = thEnvir), simplify = F),
    class = c(thClass, "temperheic"),
    units = .thSpecificUnits(generalUnits, specificUnits),
    derivedValueNames = elementNames[!elementNames %in% names(formals(thClass))]
  )
  assign("thObjects", names(newTemperheic)[sapply(newTemperheic, is.temperheic)], envir = thEnvir)
#  if(!is.null(attrs)) {
    attributeList = sapply(attrs, function(x) get(x, envir = thEnvir), simplify = F)
    for (i in 1:length(attributeList)) {
      attr(newTemperheic, attrNames[i]) = attributeList[[i]]
    }
#  }
#  attr(newTemperheic, "specificUnits") = specificUnits
  return(newTemperheic)
}

#' @export
print.temperheic = function(x, ...) {
  # star = c("", "*")
  # starVector = star[as.integer(names(x) %in% attributes(x)[["derivedValues"]])+1]
  isThObject = sapply(x, is.temperheic)
  x[isThObject] = NULL
  derivedNames = attributes(x)[["derivedValueNames"]]
  title = c("User-specified values:\n", "Derived values:\n")
  inDerived = c(FALSE, TRUE)
  for (i in 1:2) {
    printEm = (names(x) %in% derivedNames) == inDerived[i]
    if(any(printEm)){
      cat(title[i])
      cat(paste0("  ", names(x[printEm]), " = ", x[printEm], " (", attributes(x)[["units"]][printEm], ")\n"), sep="")
    }
  }
  if(length(attr(x, "thObjectNames")) > 0) {
    cat("thObjects:", paste0(attr(x, "thObjectNames"), collapse = ", "), "\n")
  }
  cat("Attributes:", paste0(names(attributes(x)), collapse = ", "))
}

#' @export
print.thUnits = function(x, ...) {
  cat(paste0(names(x), " = '", x, "'", collapse = "; "))
}

# CONSIDER ADDING "is" FUNCTIONS FOR "observed" AND "specified" CHILD CLASSES.

#' @export
is.temperheic = function(x) {
  return(is(x, "temperheic"))
}

#' @export
is.thUnits = function(x) {
  return(is(x, "thUnits"))
}

#' @export
is.thAquifer = function(x) {
  return(is(x, "thAquifer"))
}

#' @export
is.thBoundary = function(x) {
  return(is(x, "thBoundary"))
}

#' @export
is.thHydro = function(x) {
  return(is(x, "thHydro"))
}

#' @export
is.thSignal = function(x) {
  return(is(x, "thSignal"))
}