R/Barton_Choubey_plot.R

In maxgav13/GMisc: Miscellaneous Functions for Geology (and its areas)

#' @title Barton-Choubey criterion plot
#' @description Draws the Barton-Choubey criterion for the given data and for the given level of stress (unit weight and depth).
#' @param BC The result from running \code{Barton_Choubey()}
#' @param units The units to use in the plot, either "kPa" or "MPa". Default is "kPa". If "MPa" is desired you can use units = "" or units = "MPa"
#' @export
#' @return A ggplot object for the Barton-Choubey criterion
#' @import ggplot2
#' @references Barton, N. & Choubey, V. (1977). The shear strength of rock joints in theory and practice. Rock Mechanichs, 10: 1-54.
#' @examples
#' JRC = 10
#' JCS = 30
#' phi.r = 26
#' unit.weight = 18.6
#' depth = 40
#' BC = Barton_Choubey(JRC, JCS, phi.r, unit.weight, depth)
#' Barton_Choubey_plot(BC)
#'
Barton_Choubey_plot = function(BC, units = "kPa") {

  unit.factor = ifelse(units == "kPa", 1000, 1)

  c = BC$parameters$c
  phi = BC$parameters$phi

  if (units == "kPa") {
    shear = function(x) {c*unit.factor + tan(rads(phi)) * x}
  } else {
    shear = function(x) {c + tan(rads(phi)) * x}
  }

  labelphi = sprintf("phi*minute == %0.1f*degree", phi)
  if (units == "kPa") {
    labelc = sprintf("c*minute == %0.1f~kPa", c*unit.factor)
  } else {
    labelc = sprintf("c*minute == %0.4f~MPa", c)
  }

  if (units == "kPa") {
    q = ggplot(BC$dat, aes(.data$sig_n*unit.factor, .data$tau*unit.factor)) +
      geom_line(col = "red") +
      geom_point(data = BC$stress.level, aes(.data$x*unit.factor,
                                             .data$y*unit.factor),
                 col = "blue", size = 2) +
      stat_function(fun = shear, col = "blue") +
      labs(x = expression(paste("Normal stress ", sigma, " (kPa)")),
           y = expression(paste("Shear stress ", tau, " (kPa)"))) +
      annotate("text", x = max(BC$dat$sig_n*unit.factor)/4, y = max(BC$dat$tau*unit.factor)/1.25,
               label = labelphi, parse = T)+
      annotate("text", x = max(BC$dat$sig_n*unit.factor)/4, y = 0.925*max(BC$dat$tau*unit.factor)/1.25,
               label = labelc, parse = T)+
      theme_bw()+
      theme(axis.title = element_text(size = 16),
            axis.text = element_text(size = 14))
  } else {
    q = ggplot(BC$dat, aes(.data$sig_n, .data$tau)) +
      geom_line(col = "red") +
      geom_point(data = BC$stress.level, aes(.data$x, .data$y), col = "blue", size = 2) +
      stat_function(fun = shear, col = "blue") +
      labs(x = expression(paste("Normal stress ", sigma, " (MPa)")),
           y = expression(paste("Shear stress ", tau, " (MPa)"))) +
      annotate("text", x = max(BC$dat$sig_n)/4, y = max(BC$dat$tau)/1.25,
               label = labelphi, parse = T)+
      annotate("text", x = max(BC$dat$sig_n)/4, y = 0.925*max(BC$dat$tau)/1.25,
               label = labelc, parse = T)+
      theme_bw()+
      theme(axis.title = element_text(size = 16),
            axis.text = element_text(size = 14))
  }
  return(q)
}