R/plot_eha_log2.R

#' @title Plots an log2 spaced EHA spectra
#'
#' @description Plot wavelet scalogram using the outcome of the \code{\link{eha_log2}} function.
#'
#'@param eha_log2 eha_log2 object created using the \code{\link{eha_log2}} function.
#'@param plot_opt plot options are  "Power", "Amplitude", "Probability" or "F_test"
#'@param lowerPeriod Lowest period value which will be plotted
#'@param upperPeriod Highest period value which will be plotted
#'@param n.levels  Number of color levels \code{Default=100}.
#'@param palette_name Name of the color palette which is used for plotting.
#'The color palettes than can be chosen depends on which the R package is specified in
#'the color_brewer parameter. The included R packages from which palettes can be chosen
#'from are; the 'RColorBrewer', 'grDevices', 'ColorRamps' and 'Viridis' R packages.
#'There are many options to choose from so please
#'read the documentation of these packages \code{Default=rainbow}.
#'The R package 'viridis' has the color palette options: “magma”, “plasma”,
#'“inferno”, “viridis”, “mako”, and “rocket”  and “turbo”
#'To see the color palette options of the The R pacakge 'RColorBrewer' run
#'the RColorBrewer::brewer.pal.info() function
#'The R package 'colorRamps' has the color palette options:"blue2green",
#'"blue2green2red", "blue2red",	"blue2yellow", "colorRamps",	"cyan2yellow",
#'"green2red", "magenta2green", "matlab.like", "matlab.like2" and	"ygobb"
#'The R package 'grDevices' has the built in  palette options:"rainbow",
#'"heat.colors", "terrain.colors","topo.colors" and "cm.colors"
#'To see even more color palette options of the The R pacakge 'grDevices' run
#'the grDevices::hcl.pals() function
#'The R package 'scico' has the color palette options: “acton”, “bam”,“bamako”,
#' “bamO”, “batlow”, “batlowK”,“batlowW”,“berlin”,“bilbao”,”broc”,”brocO”,
#' ”buda”,”bukavu”,”cork”,”CorkO”,”davos”,”devon”,”fes”,”Glasgow”,”grayC”,
#' “hawaii”,”imola”,”lajolla”,”lapaz”,”lipari”,”lisbon”,”manague”,”navia”,
#' ”nuuk”,”oleron”,”oslo”,”roma”,”romaO”,”Tofino”,”Tokyo”,”turku”,”Vanimo”,
#' ”vik”,”vikO”
#'The R package 'Viridis' has the color palette options: “magma”, “plasma”,
#'“inferno”, “viridis”, “mako”, and “rocket”  and “turbo”
#'@param color_brewer Name of the R package from which the color palette is chosen from.
#'The included R packages from which palettes can be chosen
#'are; the RColorBrewer, grDevices, ColorRamps,scico and Viridis R packages.
#'There are many options to choose from so please
#'read the documentation of these packages. "\code{Default=grDevices}
#'@param useRaster Plot as a raster or vector image \code{Default=TRUE}.
#'WARNING plotting as a vector image is computationally intensive.
#'@param periodlab Label for the y-axis \code{Default="Period (metres)"}.
#'@param x_lab Label for the x-axis \code{Default="depth (metres)"}.
#'@param keep_editable Keep option to add extra features after plotting  \code{Default=FALSE}
#'@param dev_new Opens a new plotting window to plot the plot, this guarantees a "nice" looking plot however when plotting in an R markdown
#'document the plot might not plot  \code{Default=TRUE}
#'@param plot_dir The direction of the proxy record which is assumed for tuning if time increases with increasing depth/time values
#'(e.g. bore hole data which gets older with increasing depth ) then plot_dir should be set to TRUE
#'if time decreases with depth/time values (eg stratospheric logs where 0m is the bottom of the section)
#'then plot_dir should be set to FALSE \code{plot_dir=TRUE}
#'@param add_lines Add  lines to the wavelet plot input should be matrix with first axis being depth/time the columns after that
#'should be period values  \code{Default=NULL}
#'@param add_points Add points to the wavelet plot input should be matrix with first axis being depth/time and columns after that
#'should be period values \code{Default=NULL}
#'@param add_abline_h Add horizontal lines to the plot. Specify the lines as a vector e.g. c(2,3,5,6)  \code{Default=NULL}
#'@param add_abline_v Add vertical lines to the plot. Specify the lines as a vector e.g. c(2,3,5,6)  \code{Default=NULL}
#'@param add_MTM_peaks Add the MTM peak periods as horizontal lines \code{Default=FALSE}
#'@param add_data Plot the data on top of the wavelet \code{Default=TRUE}
#'@param add_avg Plot the average wavelet spectral power to the side of the wavelet \code{Default=FALSE}
#'@param add_MTM Add the MTM  plot next to the wavelet plot \code{Default=FALSE}
#'@param mtm_siglvl select the significance level  (0-1) for the MTM spectrum \code{Default=0.95}
#'@param demean_mtm Remove mean from data before conducting the MTM analysis \code{Default=TRUE}
#'@param detrend_mtm Remove mean from data before conducting the MTM analysis \code{Default=TRUE}
#'@param padfac_mtm Pad factor for the MTM analysis \code{Default=5}
#'@param tbw_mtm time bandwidth product of the MTM analysis  \code{Default=3}
#'@param plot_horizontal plot the wavelet horizontal or vertical eg y axis is depth or y axis power \code{Default=TRUE}
#'@param pval_abline p value of the abline
#'@param pval_cutoff p value cutoff below which no data is displayed
#'
#' @return
#' The output is a plot of a EHA spectra.
#' if add_MTM_peaks = TRUE then the output of the MTM analysis will given as matrix
#'
#' @author
#' Code based on the wt.image" functions of the 'WaveletComp' R package
#' The EHA and MTM analysis parts are from the astrochron R package of Meyers et al., (2012)
#'
#' @references
#'Angi Roesch and Harald Schmidbauer (2018). WaveletComp: Computational
#'Wavelet Analysis. R package version 1.1.
#'\url{https://CRAN.R-project.org/package=WaveletComp}
#'
#'S.R. Meyers, 2012, Seeing Red in Cyclic Stratigraphy: Spectral Noise Estimation for
#'Astrochronology: Paleoceanography, 27, PA3228, <doi:10.1029/2012PA002307>
#'
#'S.R. Meyers, 2019 Cyclostratigraphy and the problem of astrochronologic testing, Earth-Science Reviews,
#'Volume 190, <doi.org/10.1016/j.earscirev.2018.11.015.>
#'
#'@examples
#' \donttest{
#'#Example 1. A plot of a wavelet spectra using the Total Solar Irradiance
#'# data set of Steinhilber et al., (2012)
#'
#'TSI_eha_log2 <- eha_log2(data = TSI,
#'win =  8192,
#'tbw = 4,
#'demean = TRUE,
#'detrend = TRUE,
#'upperPeriod = 8192,
#'lowerPeriod = 50,
#'pad = NULL,
#'padding = "noise")
#'
#'plot_eha_log2(
#'eha_log2 = TSI_eha_log2,
#'plot_opt = "Amplitude",
#'lowerPeriod = 50,
#'upperPeriod = 8192,
#'n.levels = 100,
#'palette_name = "rainbow",
#'color_brewer = "grDevices",
#'useRaster = TRUE,
#'periodlab = "Period (metres)",
#'x_lab = "depth (metres)",
#'keep_editable = FALSE,
#'dev_new = TRUE,
#'plot_dir = TRUE,
#'add_lines = NULL,
#'add_points = NULL,
#'add_abline_h = NULL,
#'add_abline_v = NULL,
#'add_MTM_peaks = FALSE,
#'add_data = TRUE,
#'add_avg = FALSE,
#'add_MTM = FALSE,
#'mtm_siglvl = 0.95,
#'demean_mtm = TRUE,
#'detrend_mtm = TRUE,
#'padfac_mtm = 5,
#'tbw_mtm = 3,
#'plot_horizontal = TRUE)
#'
#'#Example 2. A plot of a wavelet spectra using the magnetic susceptibility
#'#data set of Pas et al., (2018)
#'mag_eha_log2 <- eha_log2(data = mag,
#'win =  50,
#'tbw = 4,
#'demean = TRUE,
#'detrend = TRUE,
#'upperPeriod = 50,
#'lowerPeriod = 1,
#'pad = NULL,
#'padding = "noise")
#'
#'plot_eha_log2(
#'eha_log2 = mag_eha_log2,
#'plot_opt = "Amplitude",
#'lowerPeriod = 1,
#'upperPeriod = 50,
#'n.levels = 100,
#'palette_name = "rainbow",
#'color_brewer = "grDevices",
#'useRaster = TRUE,
#'periodlab = "Period (metres)",
#'x_lab = "depth (metres)",
#'keep_editable = FALSE,
#'dev_new = TRUE,
#'plot_dir = TRUE,
#'add_lines = NULL,
#'add_points = NULL,
#'add_abline_h = NULL,
#'add_abline_v = NULL,
#'add_MTM_peaks = FALSE,
#'add_data = TRUE,
#'add_avg = FALSE,
#'add_MTM = FALSE,
#'mtm_siglvl = 0.95,
#'demean_mtm = TRUE,
#'detrend_mtm = TRUE,
#'padfac_mtm = 5,
#'tbw_mtm = 3,
#'plot_horizontal = TRUE)
#'
#'
#'#Example 3. A plot of a wavelet spectra using the greyscale
#'# data set of Zeeden et al., (2013)
#'grey_eha_log2 <- eha_log2(data = grey,
#'win =  20,
#'tbw = 4,
#'demean = TRUE,
#'detrend = TRUE,
#'upperPeriod = 20,
#'lowerPeriod = 1,
#'pad = NULL,
#'padding = "noise")
#'
#'
#'plot_eha_log2(
#'eha_log2 = grey_eha_log2,
#'plot_opt = "Amplitude",
#'lowerPeriod = 1,
#'upperPeriod = 20,
#'n.levels = 100,
#'palette_name = "rainbow",
#'color_brewer = "grDevices",
#'useRaster = TRUE,
#'periodlab = "Period (metres)",
#'x_lab = "depth (metres)",
#'keep_editable = FALSE,
#'dev_new = TRUE,
#'plot_dir = TRUE,
#'add_lines = NULL,
#'add_points = NULL,
#'add_abline_h = NULL,
#'add_abline_v = NULL,
#'add_MTM_peaks = FALSE,
#'add_data = TRUE,
#'add_avg = FALSE,
#'add_MTM = FALSE,
#'mtm_siglvl = 0.95,
#'demean_mtm = TRUE,
#'detrend_mtm = TRUE,
#'padfac_mtm = 5,
#'tbw_mtm = 3,
#'plot_horizontal = TRUE)
#'}
#'
#' @export
#' @importFrom stats quantile
#' @importFrom graphics par
#' @importFrom graphics image
#' @importFrom graphics axis
#' @importFrom graphics mtext
#' @importFrom graphics text
#' @importFrom graphics box
#' @importFrom graphics polygon
#' @importFrom graphics layout
#' @importFrom graphics title
#' @importFrom grDevices rgb
#' @importFrom astrochron mtm
#' @importFrom DescTools Closest
#' @importFrom graphics abline
#' @importFrom RColorBrewer brewer.pal.info
#' @importFrom RColorBrewer brewer.pal
#' @importFrom grDevices colorRampPalette
#' @importFrom colorRamps blue2green
#' @importFrom colorRamps blue2green2red
#' @importFrom colorRamps blue2red
#' @importFrom colorRamps blue2yellow
#' @importFrom colorRamps cyan2yellow
#' @importFrom colorRamps green2red
#' @importFrom colorRamps magenta2green
#' @importFrom colorRamps matlab.like
#' @importFrom colorRamps matlab.like2
#' @importFrom colorRamps ygobb
#' @importFrom viridis viridis
#' @importFrom viridis magma
#' @importFrom viridis plasma
#' @importFrom viridis inferno
#' @importFrom viridis cividis
#' @importFrom viridis mako
#' @importFrom viridis rocket
#' @importFrom viridis turbo
#' @importFrom grDevices rainbow
#' @importFrom grDevices heat.colors
#' @importFrom grDevices terrain.colors
#' @importFrom grDevices topo.colors
#' @importFrom grDevices cm.colors
#' @importFrom grDevices hcl.colors
#' @importFrom RColorBrewer brewer.pal.info
#' @importFrom RColorBrewer brewer.pal
#' @importFrom grDevices colorRampPalette
#' @importFrom colorRamps blue2green
#' @importFrom colorRamps blue2green2red
#' @importFrom colorRamps blue2red
#' @importFrom colorRamps blue2yellow
#' @importFrom colorRamps cyan2yellow
#' @importFrom colorRamps green2red
#' @importFrom colorRamps magenta2green
#' @importFrom colorRamps matlab.like
#' @importFrom colorRamps matlab.like2
#' @importFrom colorRamps ygobb
#' @importFrom viridis viridis
#' @importFrom viridis magma
#' @importFrom viridis plasma
#' @importFrom viridis inferno
#' @importFrom viridis cividis
#' @importFrom viridis mako
#' @importFrom viridis rocket
#' @importFrom viridis turbo
#' @importFrom grDevices rainbow
#' @importFrom grDevices heat.colors
#' @importFrom grDevices terrain.colors
#' @importFrom grDevices topo.colors
#' @importFrom grDevices cm.colors
#' @importFrom grDevices hcl.colors
#' @importFrom scico scico

plot_eha_log2 <- function (eha_log2 = NULL, plot_opt = "Amplitude", lowerPeriod = NULL,
                           upperPeriod = NULL, n.levels = 100, palette_name = "rainbow",
                           color_brewer = "grDevices", useRaster = TRUE, periodlab = "Period (metres)",
                           x_lab = "depth (metres)", keep_editable = FALSE, dev_new = TRUE,
                           plot_dir = TRUE, add_lines = NULL, add_points = NULL, add_abline_h = NULL,
                           add_abline_v = NULL, add_MTM_peaks = FALSE, add_data = TRUE,
                           add_avg = FALSE, pval_abline = c(0.1, 0.05), pval_cutoff = c(0.1),
                           add_MTM = FALSE, mtm_siglvl = 0.95, demean_mtm = TRUE, detrend_mtm = TRUE,
                           padfac_mtm = 5, tbw_mtm = 3, plot_horizontal = TRUE)
{

  if (plot_opt == "Amplitude") {
    mat <- eha_log2$amp
    average <- eha_log2$amp.avg
  }
  else if (plot_opt == "Power") {
    mat <- eha_log2$pwr
    average <- eha_log2$pwr.avg
  }
  else if (plot_opt == "Probability") {
    mat <- eha_log2$prob
    average <- eha_log2$prob.avg
  }
  else if (plot_opt == "F_test") {
    mat <- eha_log2$f_test
    average <- eha_log2$f_test.avg
  }
  else {
    mat <- eha_log2$amp
    average <- eha_log2$amp.avg
  }
  if (keep_editable == FALSE) {
    oldpar <- par(no.readonly = TRUE)
    on.exit(par(oldpar))
  }
  maximum.level = max(mat, na.rm = TRUE)
  power_max_mat.levels = quantile(mat, na.rm = TRUE, probs = seq(from = 0,
                                                                 to = 1, length.out = n.levels + 1))
  if (color_brewer == "RColorBrewer") {
    key.cols <- rev(colorRampPalette(brewer.pal(brewer.pal.info[palette_name,
                                                                1], palette_name))(n.levels))
  }
  if (color_brewer == "colorRamps") {
    color_brewer_Sel <- paste("colorRamps::", palette_name,
                              "(n=n.levels)")
    key.cols = eval(parse(text = color_brewer_Sel))
  }
  if (color_brewer == "grDevices") {
    if (palette_name == "rainbow") {
      color_brewer_Sel <- "grDevices::rainbow(n=n.levels, start = 0, end = 0.7)"
      key.cols <- rev(eval(parse(text = color_brewer_Sel)))
    }
    else if (palette_name == "heat.colors" | palette_name ==
             "terrain.colors" | palette_name == "topo.colors" |
             palette_name == "cm.colors") {
      color_brewer_Sel <- paste("grDevices::", palette_name,
                                "(n=n.levels, start = 0, end = 1)")
      key.cols <- rev(eval(parse(text = color_brewer_Sel)))
    }
    else {
      key.cols <- hcl.colors(n = n.levels, palette = palette_name,
                             alpha = NULL, rev = FALSE, fixup = TRUE)
    }
  }

  if (color_brewer== "scico"){
    color_brewer_Sel <- paste(
      "scico::scico(n = n.levels, palette = '",
      palette_name,
      "', direction = -1)",
      sep = ""
    )
    key.cols = rev(eval(parse(text = color_brewer_Sel)))
  }



  if (color_brewer == "viridis") {
    color_brewer_Sel <- paste("viridis::", palette_name,
                              "(n=n.levels,direction = -1)")
    key.cols = rev(eval(parse(text = color_brewer_Sel)))
  }
  periodtck = 0.02
  periodtcl = 0.5
  main = NULL
  lwd = 2
  lwd.axis = 1
  legend.params = list(width = 1.2, shrink = 0.9, mar = 5.1,
                       n.ticks = 6, label.digits = 3, label.format = "f", lab = NULL,
                       lab.line = 2.5)
  key.marks = round(seq(from = 0, to = 1, length.out = legend.params$n.ticks) *
                      n.levels)
  key.labels = formatC(as.numeric(power_max_mat.levels), digits = legend.params$label.digits,
                       format = legend.params$label.format)[key.marks + 1]
  if (dev_new == TRUE & plot_horizontal == TRUE) {
    dev.new(width = 15, height = 7, noRStudioGD = TRUE)
  }
  if (dev_new == TRUE & plot_horizontal == FALSE) {
    dev.new(width = 7, height = 10, noRStudioGD = TRUE)
  }
  depth <- eha_log2$depth
  y_axis <- eha_log2$log2_period
  depth <- as.numeric(depth)
  y_axis <- as.numeric(y_axis)
  MTM_res_2 <- matrix(data = NA, ncol = 2, nrow = 1)
  if (add_MTM == TRUE | add_MTM_peaks == TRUE) {
    MTM_res_1 <- mtm(cbind(eha_log2$x, eha_log2$y), tbw = tbw_mtm,
                     padfac = padfac_mtm, output = 1, siglevel = mtm_siglvl,
                     genplot = FALSE, verbose = FALSE, demean = demean_mtm,
                     detrend = detrend_mtm)
    MTM_res_2 <- mtm(cbind(eha_log2$x, eha_log2$y), tbw = tbw_mtm,
                     padfac = padfac_mtm, output = 3, siglevel = mtm_siglvl,
                     genplot = FALSE, verbose = FALSE, demean = demean_mtm,
                     detrend = detrend_mtm)
    MTM_res_1$period <- 1/MTM_res_1[, 1]
    MTM_res_1$Peak <- 0
    if (is.na(MTM_res_2[1, 2]) == FALSE) {
      for (i in 1:nrow(MTM_res_2)) {
        row_nr <- DescTools::Closest(MTM_res_1[, 1],
                                     MTM_res_2[i, 1], which = TRUE)
        row_nr <- row_nr[1]
        MTM_res_1[row_nr, 10] <- 1
      }
    }
    mtm_res <- MTM_res_1
  }
  if (plot_dir != TRUE) {
    xlim_vals = rev(c(min(eha_log2$x), max(eha_log2$x)))
  }
  else {
    xlim_vals = c(min(eha_log2$x), max(eha_log2$x))
  }
  if (is.null(lowerPeriod) == TRUE) {
    lowerPeriod <- min(eha_log2$Period)
  }
  if (is.null(upperPeriod) == TRUE) {
    upperPeriod <- max(eha_log2$Period)
  }
  ylim_vals = c(lowerPeriod, upperPeriod)
  if (add_data == TRUE & add_avg == FALSE & add_MTM == FALSE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 3, 0, 2), nrow = 2, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      0.25), widths = c(1, 4))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(4, 4, 2, 0))
    plot(x = eha_log2$y, y = eha_log2$x, type = "l", yaxs = "i",
         xlab = "proxy value", ylab = x_lab, ylim = xlim_vals)
    if (is.null(add_abline_h) != TRUE) {
      abline(h = add_abline_h)
    }
    par(new = FALSE, mar = c(3, 0, 2, 2))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "")
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 2, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 0, 2, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = "", yaxt = "n",
          xaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    period.tick = unique(trunc(  as.numeric(eha_log2$log2_period)))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == FALSE & add_avg == FALSE & add_MTM == FALSE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 2), nrow = 2, ncol = 1,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      4), widths = c(1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(2, 4, 2, 3))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 2, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 4, 2, 3))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = x_lab,
          xaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
  }
  if (add_data == FALSE & add_avg == FALSE & add_MTM == TRUE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(2, 1, 3, 0, 4, 0, 5, 0), nrow = 4,
                            ncol = 2, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      1, 1, 4), widths = c(4, 1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(2, 2, 2, 3))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 2, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(mar = c(0, 4, 2, 2))
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 2], type = "l",
         xaxs = "i", xaxt = "n", ylab = "MTM power", xlab = "",
         log = "x", xlim = ylim_vals)
    for (i in 5:8) {
      lines(x = 1/MTM_res_1[, 1], y = MTM_res_1[, i],
            xlim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(mar = c(0, 4, 0, 2))
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 4], type = "l",
         xaxs = "i", xaxt = "n", ylab = "Ar. conf lvl", xlab = "",
         log = "x", ylim = c(80, 101), xaxs = "i", xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 3], type = "l",
         xaxs = "i", xaxt = "n", ylab = "Har. conf lvl",
         xlab = "", log = "xy", ylim = c(80, 101), yaxs = "i",
         xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(new = FALSE, mar = c(4, 4, 0, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = x_lab,
          xaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(as.numeric(eha_log2$log2_period)))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == FALSE & add_avg == TRUE & add_MTM == TRUE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(2, 1, 3, 0, 4, 0, 5, 0, 6,
                              0), nrow = 5, ncol = 2, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      1, 1, 1, 4), widths = c(4, 1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(2, 2, 2, 3))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 2, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(mar = c(0, 4, 2, 2))
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 2], type = "l",
         xaxs = "i", xaxt = "n", ylab = "MTM power", xlab = "",
         log = "x", xlim = ylim_vals)
    for (i in 5:8) {
      lines(x = 1/MTM_res_1[, 1], y = MTM_res_1[, i],
            xlim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(mar = c(0, 4, 0, 2))
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 4], type = "l",
         xaxs = "i", xaxt = "n", ylab = "Ar. conf lvl", xlab = "",
         log = "x", ylim = c(80, 101), xaxs = "i", xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 3], type = "l",
         xaxs = "i", xaxt = "n", ylab = "Har. conf lvl",
         xlab = "", log = "xy", ylim = c(80, 101), yaxs = "i",
         xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = average, x = 2^eha_log2$log2_period, log = "x",
         type = "l", xaxs = "i", xaxt = "n", ylab = "Wt. power",
         xaxs = "i", xlim = ylim_vals)
    par(new = FALSE, mar = c(4, 4, 0, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = x_lab,
          xaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == TRUE & add_avg == TRUE & add_MTM == TRUE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 2, 0, 3, 0, 4, 0, 5, 6,
                              7), nrow = 5, ncol = 2, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      1, 1, 1, 4), widths = c(1, 4))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 1, 2, 6), xpd = FALSE)
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 1, at = mean(key.marks), line = 4,
          font = par()$font.axis, cex = par()$cex.axis, las = 1)
    box(lwd = lwd.axis)
    par(mar = c(0, 0, 2, 2), xpd = FALSE)
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 2], type = "l",
         xaxs = "i", xaxt = "n", ylab = "MTM power", xlab = "",
         log = "x", xlim = ylim_vals)
    for (i in 5:8) {
      lines(x = 1/MTM_res_1[, 1], y = MTM_res_1[, i],
            xlim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(mar = c(0, 0, 0, 2), xpd = FALSE)
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 4], type = "l",
         xaxs = "i", xaxt = "n", ylab = "", xlab = "", log = "x",
         ylim = c(80, 101), xaxs = "i", xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    title(ylab = "Ar. conf lvl", xpd = NA)
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 3], type = "l",
         xaxs = "i", xaxt = "n", ylab = "", xlab = "", log = "xy",
         ylim = c(80, 101), yaxs = "i", xlim = ylim_vals)
    title(ylab = "Har. conf lvl", xpd = NA)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = average, x = 2^eha_log2$log2_period, log = "x",
         type = "l", xaxs = "i", xaxt = "n", ylab = "Wt. power",
         xaxs = "i", xlim = ylim_vals)
    title(ylab = "Wt. power", xpd = NA)
    par(mar = c(4, 4, 0, 0), xpd = FALSE)
    plot(x = eha_log2$y, y = eha_log2$x, type = "l", yaxs = "i",
         xlab = "proxy value", ylab = x_lab, ylim = xlim_vals)
    par(new = FALSE, mar = c(4, 0, 0, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, xaxt = "n",
          yaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(as.numeric(eha_log2$log2_period)))
    period.tick[period.tick < log2(as.numeric(eha_log2$log2_period[1]))] = NA
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == TRUE & add_avg == TRUE & add_MTM == FALSE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 2, 3, 4), nrow = 2, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      4), widths = c(1, 4))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(2, 0.5, 2, 6), xpd = NA)
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(2, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 2, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 1, at = 1, line = 0.5, font = par()$font.axis,
          cex = par()$cex.axis, las = 1, xpd = NA)
    box(lwd = lwd.axis)
    par(mar = c(0, 0, 2, 2))
    plot(y = average, x = 2^eha_log2$log2_period, log = "x",
         type = "l", xaxs = "i", xaxt = "n", ylab = "Wt. power",
         xlab = "", xaxs = "i", xlim = ylim_vals)
    title(ylab = "Wt. power", xpd = NA)
    par(mar = c(4, 4, 0, 0), xpd = TRUE)
    plot(x = eha_log2$y, y = eha_log2$x, type = "l", yaxs = "i",
         xlab = "proxy value", ylab = x_lab, ylim = xlim_vals)
    par(new = FALSE, mar = c(4, 0, 0, 2), xpd = FALSE)
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = "", xaxt = "n",
          yaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == TRUE & add_avg == FALSE & add_MTM == TRUE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 2, 0, 3, 0, 4, 5, 6), nrow = 4,
                            ncol = 2, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      1, 1, 4), widths = c(1, 4))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 1, 2, 6), xpd = NA)
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 1, at = mean(key.marks), line = 4,
          font = par()$font.axis, cex = par()$cex.axis, las = 1)
    box(lwd = lwd.axis)
    par(mar = c(0, 0, 2, 2))
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 2], type = "l",
         xaxs = "i", xaxt = "n", ylab = "MTM power", xlab = "",
         log = "x", xlim = ylim_vals)
    for (i in 5:8) {
      lines(x = 1/MTM_res_1[, 1], y = MTM_res_1[, i],
            xlim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(mar = c(0, 0, 0, 2), xpd = FALSE)
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 4], type = "l",
         xaxs = "i", xaxt = "n", ylab = "", xlab = "", log = "x",
         ylim = c(80, 101), xaxs = "i", xlim = ylim_vals)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    title(ylab = "Ar. conf lvl", xpd = NA)
    plot(x = 1/MTM_res_1[, 1], y = MTM_res_1[, 3], type = "l",
         xaxs = "i", xaxt = "n", ylab = "", xlab = "", log = "xy",
         ylim = c(80, 101), yaxs = "i", xlim = ylim_vals)
    title(ylab = "Har. conf lvl", xpd = NA)
    abline(h = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(mar = c(4, 4, 0, 0))
    plot(x = eha_log2$y, y = eha_log2$x, type = "l", yaxs = "i",
         xlab = "proxy value", ylab = x_lab, ylim = xlim_vals)
    par(new = FALSE, mar = c(4, 0, 0, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, xaxt = "n",
          yaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == FALSE & add_avg == TRUE & add_MTM == FALSE &
      plot_horizontal == FALSE) {
    layout.matrix <- matrix(c(1, 2, 0, 3), nrow = 2, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1,
                                                      4), widths = c(1, 4))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(2, 2, 2, 3))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 1, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 1, at = median(key.marks),
          line = 3, las = 2, font = par()$font.axis, cex = par()$cex.axis,
          xpd = NA)
    box(lwd = lwd.axis)
    par(mar = c(0, 4, 2, 2))
    plot(y = average, x = 2^eha_log2$log2_period, log = "x",
         type = "l", xaxs = "i", xaxt = "n", ylab = "Wt. power",
         xaxs = "i", xlim = ylim_vals)
    if (is.null(add_abline_v) != TRUE) {
      abline(v = (add_abline_v))
    }
    par(new = FALSE, mar = c(4, 4, 0, 2))
    image(y = eha_log2$depth, x = eha_log2$log2_period,
          z = (mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, xlab = periodlab, ylab = x_lab,
          xaxt = "n", main = main, ylim = xlim_vals, xlim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(1, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 1, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(y = add_lines[,
                                                       1], x = log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(y = add_points[,
                                                          1], x = log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(v = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = (add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = log2(add_abline_v))
    }
  }
  if (add_data == TRUE & add_avg == FALSE & add_MTM == TRUE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(1, 1, 1, 1, 1, 1, 2, 2, 2,
                              3, 3, 3, 3, 3, 3, 4, 5, 6), nrow = 2, ncol = 9,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(0.25,
                                                      1), widths = c(1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 4, 2, 0))
    plot(y = eha_log2$y, x = eha_log2$x, type = "l", xaxs = "i",
         xlab = "", ylab = "proxy value", xaxt = "n", xlim = xlim_vals)
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(3, 2, 2, 2), mgp = c(2, 1,
                                                  0))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "Power",
          ylab = "")
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.2)
    mtext(key.labels, side = 1, at = key.marks, line = 0.1,
          las = 2, cex = 0.75)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 4, 0, 0), mgp = c(2, 1,
                                                  0))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(4, 0, 0, 0.5))
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 2], type = "l",
         yaxs = "i", yaxt = "n", xlab = "MTM power", ylab = "",
         log = "y", ylim = ylim_vals)
    for (i in 5:8) {
      lines(y = 1/MTM_res_1[, 1], x = MTM_res_1[, i],
            ylim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 4], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Ar. conf lvl", ylab = "",
         log = "y", xlim = c(80, 101), xaxs = "i", ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 3], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Har. conf lvl",
         ylab = "", log = "xy", xlim = c(80, 101), xaxs = "i",
         ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
  }
  if (add_data == TRUE & add_avg == TRUE & add_MTM == FALSE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(1, 2, 4, 3), nrow = 2, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(0.25,
                                                      1), widths = c(8, 2))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 4, 2, 0))
    plot(y = eha_log2$y, x = eha_log2$x, type = "l", xaxs = "i",
         xlab = "", ylab = "proxy value", xaxt = "n", xlim = xlim_vals)
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(3, 2, 2, 2), mgp = c(2, 1,
                                                  0))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "Power",
          ylab = "")
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.2)
    mtext(key.labels, side = 1, at = key.marks, line = 0.1,
          las = 2, cex = 0.75)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 0, 0, 0.5))
    plot(x = average, y = 2^eha_log2$log2_period, log = "y",
         type = "l", yaxs = "i", yaxt = "n", xlab = "Wt. power",
         xaxs = "i", ylim = ylim_vals)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "black", lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = add_abline_h)
    }
    par(new = FALSE, mar = c(4, 4, 0, 0), xpd = FALSE)
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_data == TRUE & add_avg == FALSE & add_MTM == FALSE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(1, 0, 3, 2), nrow = 2, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(0.25,
                                                      1), widths = c(8, 2))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 4, 2, 2))
    plot(y = eha_log2$y, x = eha_log2$x, type = "l", xaxs = "i",
         xlab = "", ylab = "proxy value", xaxt = "n", xlim = xlim_vals)
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(4, 0, 2, 5))
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 4, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 4, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 4, 0, 2))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_data == TRUE & add_avg == TRUE & add_MTM == TRUE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(1, 1, 1, 1, 1, 1, 2, 2, 2,
                              2, 3, 3, 3, 3, 3, 3, 4, 5, 6, 7), nrow = 2, ncol = 10,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(0.25,
                                                      1), widths = c(1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(0, 4, 2, 0))
    plot(y = eha_log2$y, x = eha_log2$x, type = "l", xaxs = "i",
         xlab = "", ylab = "proxy value", xaxt = "n", xlim = xlim_vals)
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(3, 2, 2, 2), mgp = c(2, 1,
                                                  0))
    image(x = seq(from = 0, to = n.levels), y = 1, z = t(matrix(power_max_mat.levels,
                                                                nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "Power",
          ylab = "")
    axis(1, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.2)
    mtext(key.labels, side = 1, at = key.marks, line = 0.1,
          las = 2, cex = 0.75)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 4, 0, 0), mgp = c(2, 1,
                                                  0))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
    par(new = FALSE, mar = c(4, 0, 0, 0.5))
    plot(x = average, y = 2^eha_log2$log2_period, log = "y",
         type = "l", yaxs = "i", yaxt = "n", xlab = "Wt. power",
         xaxs = "i", ylim = ylim_vals)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "black", lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = add_abline_h)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 2], type = "l",
         yaxs = "i", yaxt = "n", xlab = "MTM power", ylab = "",
         log = "y", ylim = ylim_vals)
    for (i in 5:8) {
      lines(y = 1/MTM_res_1[, 1], x = MTM_res_1[, i],
            ylim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 4], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Ar. conf lvl", ylab = "",
         log = "y", xlim = c(80, 101), xaxs = "i", ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 3], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Har. conf lvl",
         ylab = "", log = "xy", xlim = c(80, 101), xaxs = "i",
         ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
  }
  if (add_data == FALSE & add_avg == TRUE & add_MTM == FALSE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(3, 2, 1), nrow = 1, ncol = 3,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1),
                     widths = c(6, 1, 1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(4, 0, 2, 5))
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 4, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 4, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 0, 2, 0.5))
    plot(x = average, y = 2^eha_log2$log2_period, log = "y",
         type = "l", yaxs = "i", yaxt = "n", xlab = "Wt. power",
         xaxs = "i", ylim = ylim_vals)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "black", lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = add_abline_h)
    }
    par(new = FALSE, mar = c(4, 4, 2, 0), mgp = c(2, 1,
                                                  0))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_data == FALSE & add_avg == TRUE & add_MTM == TRUE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(6, 5, 4, 3, 2, 1), nrow = 1,
                            ncol = 6, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1),
                     widths = c(6, 1, 1, 1, 1, 1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(4, 0, 2, 5))
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 4, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 4, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 0, 2, 0.5))
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 3], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Har. conf lvl",
         ylab = "", log = "xy", xlim = c(80, 101), xaxs = "i",
         ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 4], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Ar. conf lvl", ylab = "",
         log = "y", xlim = c(80, 101), xaxs = "i", ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 2], type = "l",
         yaxs = "i", yaxt = "n", xlab = "MTM power", ylab = "",
         log = "y", ylim = ylim_vals)
    for (i in 5:8) {
      lines(y = 1/MTM_res_1[, 1], x = MTM_res_1[, i],
            ylim = c(min(y_axis), max(y_axis)), lty = 3,
            col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(x = average, y = 2^eha_log2$log2_period, log = "y",
         type = "l", yaxs = "i", yaxt = "n", xlab = "Wt. power",
         xaxs = "i", ylim = ylim_vals)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "black", lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = add_abline_h)
    }
    par(new = FALSE, mar = c(4, 4, 2, 0), mgp = c(2, 1,
                                                  0))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_data == FALSE & add_avg == FALSE & add_MTM == TRUE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(5, 4, 3, 2, 1), nrow = 1,
                            ncol = 5, byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1),
                     widths = c(7, 1, 1, 1, 1))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(4, 0, 2, 5))
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 4, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 4, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 0, 2, 0.5))
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 3], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Har. conf lvl",
         ylab = "", log = "xy", xlim = c(80, 101), xaxs = "i",
         ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 4], type = "l",
         yaxs = "i", yaxt = "n", xlab = "Ar. conf lvl", ylab = "",
         log = "y", xlim = c(80, 101), xaxs = "i", ylim = ylim_vals)
    abline(v = c(90, 95, 99), lty = 3, col = "grey", lwd = 2)
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    plot(y = 1/MTM_res_1[, 1], x = MTM_res_1[, 2], type = "l",
         yaxs = "i", yaxt = "n", xlab = "MTM power", ylab = "",
         log = "y", ylim = ylim_vals)
    for (i in 5:8) {
      lines(y = 1/MTM_res_1[, 1], x = MTM_res_1[, i],
            ylim = ylim_vals, lty = 3, col = "grey", lwd = 2)
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = 1/MTM_res_2[, 1], col = "red", lty = 3)
    }
    par(new = FALSE, mar = c(4, 4, 2, 0), mgp = c(2, 1,
                                                  0))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_data == FALSE & add_avg == FALSE & add_MTM == FALSE &
      plot_horizontal == TRUE) {
    layout.matrix <- matrix(c(2, 1), nrow = 1, ncol = 2,
                            byrow = TRUE)
    graphics::layout(mat = layout.matrix, heights = c(1),
                     widths = c(10, 2.25))
    power_max_mat.levels = quantile(mat, probs = seq(from = 0,
                                                     to = 1, length.out = n.levels + 1))
    par(mar = c(4, 0, 2, 5))
    image(y = seq(from = 0, to = n.levels), x = 1, z = (matrix(power_max_mat.levels,
                                                               nrow = 1)), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, yaxt = "n", xaxt = "n", xlab = "",
          ylab = "", )
    axis(4, lwd = lwd.axis, at = key.marks, labels = NA,
         tck = 0.02, tcl = 1.24)
    mtext(key.labels, side = 4, at = key.marks, line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    mtext(c("Power"), side = 4, at = mean(key.marks), line = 0.5,
          las = 2, font = par()$font.axis, cex = par()$cex.axis)
    box(lwd = lwd.axis)
    par(new = FALSE, mar = c(4, 4, 2, 0.5))
    image(x = eha_log2$depth, y = eha_log2$log2_period,
          z = t(mat), col = key.cols, breaks = power_max_mat.levels,
          useRaster = TRUE, ylab = periodlab, xlab = x_lab,
          axes = TRUE, yaxt = "n", main = main, xlim = xlim_vals,
          ylim = log2(ylim_vals))
    box(lwd = lwd.axis)
    period.tick = unique(trunc(eha_log2$log2_period))
    period.tick = na.omit(period.tick)
    period.tick.label = 2^(period.tick)
    axis(2, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    axis(4, lwd = lwd.axis, at = period.tick, labels = NA,
         tck = periodtck, tcl = periodtcl)
    mtext(period.tick.label, side = 2, at = period.tick,
          las = 2, line = par()$mgp[2] - 0.5, font = par()$font.axis,
          cex = par()$cex.axis)
    if (is.null(add_lines) != TRUE) {
      for (i in 2:ncol(add_lines)) lines(add_lines[, 1],
                                         log2(add_lines[, i]))
    }
    if (is.null(add_points) != TRUE) {
      for (i in 2:ncol(add_points)) points(add_points[,
                                                      1], log2(add_points[, i]))
    }
    if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
        FALSE) {
      abline(h = log2(1/MTM_res_2[, 1]), col = "black",
             lty = 3)
    }
    if (is.null(add_abline_h) != TRUE) {
      abline(h = log2(add_abline_h))
    }
    if (is.null(add_abline_v) != TRUE) {
      abline(v = add_abline_v)
    }
  }
  if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) == FALSE) {
    return(invisible(mtm_res))
  }
}


### old function keep for now #

#
# plot_eha_log2 <- function (eha_log2 = NULL,
#                            plot_opt = "Amplitude",
#                            lowerPeriod = NULL,
#                            upperPeriod = NULL,
#                            n.levels = 100,
#                            palette_name = "rainbow",
#                            color_brewer = "grDevices",
#                            useRaster = TRUE,
#                            periodlab = "Period (metres)",
#                            x_lab = "depth (metres)",
#                            keep_editable = FALSE,
#                            dev_new = TRUE,
#                            plot_dir = TRUE,
#                            add_lines = NULL,
#                            add_points = NULL,
#                            add_abline_h = NULL,
#                            add_abline_v = NULL,
#                            add_MTM_peaks = FALSE,
#                            add_data = TRUE,
#                            add_avg = FALSE,
#                            pval_abline = c(0.1, 0.05),
#                            pval_cutoff = c(0.1),
#                            add_MTM = FALSE,
#                            mtm_siglvl = 0.95,
#                            demean_mtm = TRUE,
#                            detrend_mtm = TRUE,
#                            padfac_mtm = 5,
#                            tbw_mtm = 3,
#                            plot_horizontal = TRUE)
# {
#   if (plot_opt == "Amplitude") {
#     mat <- eha_log2$amp
#     average <- eha_log2$amp.avg} else if (plot_opt == "Power") {
#       mat <- eha_log2$pwr
#       average <- eha_log2$pwr.avg} else if (plot_opt == "Probability") {
#         mat <- eha_log2$prob
#         average <- eha_log2$prob.avg} else if (plot_opt == "F_test") {
#           mat <- eha_log2$f_test
#           average <- eha_log2$f_test.avg} else{
#             mat <- eha_log2$amp
#             average <- eha_log2$amp.avg
#           }
#
#
#   if (keep_editable == FALSE) {
#     oldpar <- par(no.readonly = TRUE)
#     on.exit(par(oldpar))
#   }
#   maximum.level = max(mat,na.rm=TRUE)
#   power_max_mat.levels = quantile(mat,na.rm=TRUE, probs = seq(
#     from = 0,
#     to = 1,
#     length.out = n.levels + 1
#   ))
#   if (color_brewer == "RColorBrewer") {
#     key.cols <- rev(colorRampPalette(brewer.pal(brewer.pal.info[palette_name, 1], palette_name))(n.levels))
#   }
#   if (color_brewer == "colorRamps") {
#     color_brewer_Sel <- paste("colorRamps::", palette_name, "(n=n.levels)")
#     key.cols = eval(parse(text = color_brewer_Sel))
#   }
#   if (color_brewer == "grDevices") {
#     if (palette_name == "rainbow") {
#       color_brewer_Sel <- "grDevices::rainbow(n=n.levels, start = 0, end = 0.7)"
#       key.cols <- rev(eval(parse(text = color_brewer_Sel)))
#     }
#     else if (palette_name == "heat.colors" | palette_name ==
#              "terrain.colors" | palette_name == "topo.colors" |
#              palette_name == "cm.colors") {
#       color_brewer_Sel <- paste("grDevices::",
#                                 palette_name,
#                                 "(n=n.levels, start = 0, end = 1)")
#       key.cols <- rev(eval(parse(text = color_brewer_Sel)))
#     }
#     else {
#       key.cols <- hcl.colors(
#         n = n.levels,
#         palette = palette_name,
#         alpha = NULL,
#         rev = FALSE,
#         fixup = TRUE
#       )
#     }
#   }
#   if (color_brewer == "viridis") {
#     color_brewer_Sel <- paste("viridis::", palette_name, "(n=n.levels,direction = -1)")
#     key.cols = rev(eval(parse(text = color_brewer_Sel)))
#   }
#   periodtck = 0.02
#   periodtcl = 0.5
#   main = NULL
#   lwd = 2
#   lwd.axis = 1
#   legend.params = list(
#     width = 1.2,
#     shrink = 0.9,
#     mar = 5.1,
#     n.ticks = 6,
#     label.digits = 3,
#     label.format = "f",
#     lab = NULL,
#     lab.line = 2.5
#   )
#   key.marks = round(seq(
#     from = 0,
#     to = 1,
#     length.out = legend.params$n.ticks
#   ) *
#     n.levels)
#   key.labels = formatC(
#     as.numeric(power_max_mat.levels),
#     digits = legend.params$label.digits,
#     format = legend.params$label.format
#   )[key.marks + 1]
#   if (dev_new == TRUE & plot_horizontal == TRUE) {
#     dev.new(width = 15,
#             height = 7,
#             noRStudioGD = TRUE)
#   }
#   if (dev_new == TRUE & plot_horizontal == FALSE) {
#     dev.new(width = 7,
#             height = 10,
#             noRStudioGD = TRUE)
#   }
#
#   depth <-  eha_log2$depth
#   y_axis <- eha_log2$log2_period
#   depth <- as.numeric(depth)
#   y_axis <- as.numeric(y_axis)
#   MTM_res_2 <- matrix(data = NA,
#                       ncol = 2,
#                       nrow = 1)
#   if (add_MTM == TRUE | add_MTM_peaks == TRUE) {
#     MTM_res_1 <- mtm(
#       cbind(eha_log2$x, eha_log2$y),
#       tbw = tbw_mtm,
#       padfac = padfac_mtm,
#       output = 1,
#       siglevel = mtm_siglvl,
#       genplot = FALSE,
#       verbose = FALSE,
#       demean = demean_mtm,
#       detrend = detrend_mtm
#     )
#     MTM_res_2 <- mtm(
#       cbind(eha_log2$x, eha_log2$y),
#       tbw = tbw_mtm,
#       padfac = padfac_mtm,
#       output = 3,
#       siglevel = mtm_siglvl,
#       genplot = FALSE,
#       verbose = FALSE,
#       demean = demean_mtm,
#       detrend = detrend_mtm
#     )
#     MTM_res_1$period <- 1 / MTM_res_1[, 1]
#     MTM_res_1$Peak <- 0
#     if (is.na(MTM_res_2[1, 2]) == FALSE) {
#       for (i in 1:nrow(MTM_res_2)) {
#         row_nr <- DescTools::Closest(MTM_res_1[, 1], MTM_res_2[i, 1], which = TRUE)
#         row_nr <- row_nr[1]
#         MTM_res_1[row_nr, 10] <- 1
#       }
#     }
#     mtm_res <- MTM_res_1
#   }
#   if (plot_dir != TRUE) {
#     xlim_vals = rev(c(min(eha_log2$x), max(eha_log2$x)))
#   }else {
#     xlim_vals = c(min(eha_log2$x), max(eha_log2$x))
#   }
#   if (is.null(lowerPeriod) == TRUE) {
#     lowerPeriod <- min(eha_log2$Period)
#   }
#   if (is.null(upperPeriod) == TRUE) {
#     upperPeriod <- max(eha_log2$Period)
#   }
#
#
#   ylim_vals = c(lowerPeriod, upperPeriod)
#   if (add_data == TRUE & add_avg == FALSE & add_MTM == FALSE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(c(1, 3, 0, 2),
#                             nrow = 2,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 0.25),
#                      widths = c(1, 4))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(4, 4, 2, 0))
#     plot(
#       x = eha_log2$y,
#       y = eha_log2$x,
#       type = "l",
#       yaxs = "i",
#       xlab = "proxy value",
#       ylab = x_lab,
#       ylim = xlim_vals
#     )
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = add_abline_h)
#     }
#     par(new = FALSE, mar = c(3, 0, 2, 2))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = ""
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 2,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 0, 2, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = "",
#       yaxt = "n",
#       xaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#     period.tick = unique(trunc(eha_log2$log2_period))
#     #period.tick[period.tick < log2(eha_log2$Period[1])] = NA
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#
#   if (add_data == FALSE & add_avg == FALSE & add_MTM == FALSE &
#       plot_horizontal == FALSE) {
#     #graphics.off()
#     layout.matrix <- matrix(c(1, 2),
#                             nrow = 2,
#                             ncol = 1,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 4),
#                      widths = c(1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(2, 4, 2, 3))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 2,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 4, 2, 3))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = x_lab,
#       xaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#   }
#   if (add_data == FALSE & add_avg == FALSE & add_MTM == TRUE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(
#       c(2, 1, 3, 0, 4, 0, 5, 0),
#       nrow = 4,
#       ncol = 2,
#       byrow = TRUE
#     )
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 1, 1, 4),
#                      widths = c(4, 1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(2, 2, 2, 3))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 2,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 4, 2, 2))
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 2],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "MTM power",
#       xlab = "",
#       log = "x",
#       xlim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         x = 1 / MTM_res_1[, 1],
#         y = MTM_res_1[, i],
#         xlim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(mar = c(0, 4, 0, 2))
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 4],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Ar. conf lvl",
#       xlab = "",
#       log = "x",
#       ylim = c(80, 101),
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 3],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Har. conf lvl",
#       xlab = "",
#       log = "xy",
#       ylim = c(80, 101),
#       yaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(new = FALSE, mar = c(4, 4, 0, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = x_lab,
#       xaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#
#   if (add_data == FALSE & add_avg == TRUE & add_MTM == TRUE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(
#       c(2, 1, 3, 0, 4, 0, 5, 0, 6, 0),
#       nrow = 5,
#       ncol = 2,
#       byrow = TRUE
#     )
#     graphics::layout(
#       mat = layout.matrix,
#       heights = c(1, 1, 1, 1, 4),
#       widths = c(4, 1)
#     )
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(2, 2, 2, 3))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 2,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 4, 2, 2))
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 2],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "MTM power",
#       xlab = "",
#       log = "x",
#       xlim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         x = 1 / MTM_res_1[, 1],
#         y = MTM_res_1[, i],
#         xlim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(mar = c(0, 4, 0, 2))
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 4],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Ar. conf lvl",
#       xlab = "",
#       log = "x",
#       ylim = c(80, 101),
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 3],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Har. conf lvl",
#       xlab = "",
#       log = "xy",
#       ylim = c(80, 101),
#       yaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#
#     plot(
#       y = average,
#       x = 2^eha_log2$log2_period,
#       log = "x",
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Wt. power",
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     par(new = FALSE, mar = c(4, 4, 0, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = x_lab,
#       xaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#
#   if (add_data == TRUE & add_avg == TRUE & add_MTM == TRUE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(
#       c(1, 2, 0, 3, 0, 4, 0, 5, 6, 7),
#       nrow = 5,
#       ncol = 2,
#       byrow = TRUE
#     )
#     graphics::layout(
#       mat = layout.matrix,
#       heights = c(1, 1, 1, 1, 4),
#       widths = c(1, 4)
#     )
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 1, 2, 6), xpd = FALSE)
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 1,
#       at = mean(key.marks),
#       line = 4,
#       font = par()$font.axis,
#       cex = par()$cex.axis,
#       las = 1
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 0, 2, 2), xpd = FALSE)
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 2],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "MTM power",
#       xlab = "",
#       log = "x",
#       xlim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         x = 1 / MTM_res_1[, 1],
#         y = MTM_res_1[, i],
#         xlim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(mar = c(0, 0, 0, 2), xpd = FALSE)
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 4],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "",
#       xlab = "",
#       log = "x",
#       ylim = c(80, 101),
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     title(ylab = "Ar. conf lvl", xpd = NA)
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 3],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "",
#       xlab = "",
#       log = "xy",
#       ylim = c(80, 101),
#       yaxs = "i",
#       xlim = ylim_vals
#     )
#     title(ylab = "Har. conf lvl", xpd = NA)
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = average,
#       x = 2^eha_log2$log2_period,
#       log = "x",
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Wt. power",
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     title(ylab = "Wt. power", xpd = NA)
#     par(mar = c(4, 4, 0, 0), xpd = FALSE)
#     plot(
#       x = eha_log2$y,
#       y = eha_log2$x,
#       type = "l",
#       yaxs = "i",
#       xlab = "proxy value",
#       ylab = x_lab,
#       ylim = xlim_vals
#     )
#     par(new = FALSE, mar = c(4, 0, 0, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       xaxt = "n",
#       yaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$axis.2))
#     period.tick[period.tick < log2(eha_log2$Period[1])] = NA
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#   if (add_data == TRUE & add_avg == TRUE & add_MTM == FALSE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(c(1, 2, 3, 4),
#                             nrow = 2,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 4),
#                      widths = c(1, 4))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(2, 0.5, 2, 6), xpd = NA)
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 2,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 1,
#       at = 1,
#       line = 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis,
#       las = 1,
#       xpd = NA
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 0, 2, 2))
#     plot(
#       y = average,
#       x = 2^eha_log2$log2_period,
#       log = "x",
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Wt. power",
#       xlab = "",
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     title(ylab = "Wt. power", xpd = NA)
#     par(mar = c(4, 4, 0, 0), xpd = TRUE)
#     plot(
#       x = eha_log2$y,
#       y = eha_log2$x,
#       type = "l",
#       yaxs = "i",
#       xlab = "proxy value",
#       ylab = x_lab,
#       ylim = xlim_vals
#     )
#     par(new = FALSE,
#         mar = c(4, 0, 0, 2),
#         xpd = FALSE)
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = "",
#       xaxt = "n",
#       yaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#   if (add_data == TRUE & add_avg == FALSE & add_MTM == TRUE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(
#       c(1, 2, 0, 3, 0, 4, 5, 6),
#       nrow = 4,
#       ncol = 2,
#       byrow = TRUE
#     )
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 1, 1, 4),
#                      widths = c(1, 4))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 1, 2, 6), xpd = NA)
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 1,
#       at = mean(key.marks),
#       line = 4,
#       font = par()$font.axis,
#       cex = par()$cex.axis,
#       las = 1
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 0, 2, 2))
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 2],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "MTM power",
#       xlab = "",
#       log = "x",
#       xlim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         x = 1 / MTM_res_1[, 1],
#         y = MTM_res_1[, i],
#         xlim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(mar = c(0, 0, 0, 2), xpd = FALSE)
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 4],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "",
#       xlab = "",
#       log = "x",
#       ylim = c(80, 101),
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     title(ylab = "Ar. conf lvl", xpd = NA)
#     plot(
#       x = 1 / MTM_res_1[, 1],
#       y = MTM_res_1[, 3],
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "",
#       xlab = "",
#       log = "xy",
#       ylim = c(80, 101),
#       yaxs = "i",
#       xlim = ylim_vals
#     )
#     title(ylab = "Har. conf lvl", xpd = NA)
#     abline(
#       h = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(mar = c(4, 4, 0, 0))
#     plot(
#       x = eha_log2$y,
#       y = eha_log2$x,
#       type = "l",
#       yaxs = "i",
#       xlab = "proxy value",
#       ylab = x_lab,
#       ylim = xlim_vals
#     )
#     par(new = FALSE, mar = c(4, 0, 0, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       xaxt = "n",
#       yaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#   if (add_data == FALSE & add_avg == TRUE & add_MTM == FALSE &
#       plot_horizontal == FALSE) {
#     layout.matrix <- matrix(c(1, 2, 0, 3),
#                             nrow = 2,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1, 4),
#                      widths = c(1, 4))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(2, 2, 2, 3))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 1,
#       at = median(key.marks),
#       line = 3,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis,
#       xpd = NA
#     )
#     box(lwd = lwd.axis)
#     par(mar = c(0, 4, 2, 2))
#     plot(
#       y = average,
#       x = 2^eha_log2$log2_period,
#       log = "x",
#       type = "l",
#       xaxs = "i",
#       xaxt = "n",
#       ylab = "Wt. power",
#       xaxs = "i",
#       xlim = ylim_vals
#     )
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = (add_abline_v))
#     }
#     par(new = FALSE, mar = c(4, 4, 0, 2))
#     image(
#       y = eha_log2$depth,
#       x = eha_log2$log2_period,
#       z = (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       xlab = periodlab,
#       ylab = x_lab,
#       xaxt = "n",
#       main = main,
#       ylim = xlim_vals,
#       xlim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 1,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(y = add_lines[, 1], x = log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(y = add_points[, 1], x = log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(v = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = (add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = log2(add_abline_v))
#     }
#   }
#
#
#
#
#   if (add_data == TRUE & add_avg == FALSE & add_MTM == TRUE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(
#       c(1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 5, 6),
#       nrow = 2,
#       ncol = 9,
#       byrow = TRUE
#     )
#     graphics::layout(mat = layout.matrix,
#                      heights = c(0.25, 1),
#                      widths = c(1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 4, 2, 0))
#     plot(
#       y = eha_log2$y,
#       x = eha_log2$x,
#       type = "l",
#       xaxs = "i",
#       xlab = "",
#       ylab = "proxy value",
#       xaxt = "n",
#       xlim = xlim_vals
#     )
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE,
#         mar = c(3, 2, 2, 2),
#         mgp = c(2, 1, 0))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "Power",
#       ylab = ""
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.2
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.1,
#       las = 2,
#       cex = 0.75
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE,
#         mar = c(4, 4, 0, 0),
#         mgp = c(2, 1, 0))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE, mar = c(4, 0, 0, 0.5))
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 2],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "MTM power",
#       ylab = "",
#       log = "y",
#       ylim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         y = 1 / MTM_res_1[, 1],
#         x = MTM_res_1[, i],
#         ylim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 4],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Ar. conf lvl",
#       ylab = "",
#       log = "y",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 3],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Har. conf lvl",
#       ylab = "",
#       log = "xy",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#   }
#   if (add_data == TRUE & add_avg == TRUE & add_MTM == FALSE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(1, 2, 4, 3),
#                             nrow = 2,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(0.25, 1),
#                      widths = c(8, 2))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 4, 2, 0))
#     plot(
#       y = eha_log2$y,
#       x = eha_log2$x,
#       type = "l",
#       xaxs = "i",
#       xlab = "",
#       ylab = "proxy value",
#       xaxt = "n",
#       xlim = xlim_vals
#     )
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE,
#         mar = c(3, 2, 2, 2),
#         mgp = c(2, 1, 0))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "Power",
#       ylab = ""
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.2
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.1,
#       las = 2,
#       cex = 0.75
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 0, 0, 0.5))
#     plot(
#       x = average,
#       y = 2^eha_log2$log2_period,
#       log = "y",
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Wt. power",
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = add_abline_h)
#     }
#     par(new = FALSE,
#         mar = c(4, 4, 0, 0),
#         xpd = FALSE)
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_data == TRUE & add_avg == FALSE & add_MTM == FALSE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(1, 0, 3, 2),
#                             nrow = 2,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(0.25, 1),
#                      widths = c(8, 2))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 4, 2, 2))
#     plot(
#       y = eha_log2$y,
#       x = eha_log2$x,
#       type = "l",
#       xaxs = "i",
#       xlab = "",
#       ylab = "proxy value",
#       xaxt = "n",
#       xlim = xlim_vals
#     )
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE, mar = c(4, 0, 2, 5))
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 4,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 4,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 4, 0, 2))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_data == TRUE & add_avg == TRUE & add_MTM == TRUE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(
#       c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 5, 6, 7),
#       nrow = 2,
#       ncol = 10,
#       byrow = TRUE
#     )
#     graphics::layout(mat = layout.matrix,
#                      heights = c(0.25, 1),
#                      widths = c(1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(0, 4, 2, 0))
#     plot(
#       y = eha_log2$y,
#       x = eha_log2$x,
#       type = "l",
#       xaxs = "i",
#       xlab = "",
#       ylab = "proxy value",
#       xaxt = "n",
#       xlim = xlim_vals
#     )
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE,
#         mar = c(3, 2, 2, 2),
#         mgp = c(2, 1, 0))
#     image(
#       x = seq(from = 0, to = n.levels),
#       y = 1,
#       z = t(matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "Power",
#       ylab = ""
#     )
#     axis(
#       1,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.2
#     )
#     mtext(
#       key.labels,
#       side = 1,
#       at = key.marks,
#       line = 0.1,
#       las = 2,
#       cex = 0.75
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE,
#         mar = c(4, 4, 0, 0),
#         mgp = c(2, 1, 0))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#     par(new = FALSE, mar = c(4, 0, 0, 0.5))
#     plot(
#       x = average,
#       y = 2^eha_log2$log2_period,
#       log = "y",
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Wt. power",
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = add_abline_h)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 2],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "MTM power",
#       ylab = "",
#       log = "y",
#       ylim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         y = 1 / MTM_res_1[, 1],
#         x = MTM_res_1[, i],
#         ylim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 4],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Ar. conf lvl",
#       ylab = "",
#       log = "y",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 3],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Har. conf lvl",
#       ylab = "",
#       log = "xy",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#   }
#   if (add_data == FALSE & add_avg == TRUE & add_MTM == FALSE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(3, 2, 1),
#                             nrow = 1,
#                             ncol = 3,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1),
#                      widths = c(6, 1, 1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(4, 0, 2, 5))
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 4,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 4,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 0, 2, 0.5))
#     plot(
#       x = average,
#       y = 2^eha_log2$log2_period,
#       log = "y",
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Wt. power",
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = add_abline_h)
#     }
#     par(new = FALSE,
#         mar = c(4, 4, 2, 0),
#         mgp = c(2, 1, 0))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_data == FALSE & add_avg == TRUE & add_MTM == TRUE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(6, 5, 4, 3, 2, 1),
#                             nrow = 1,
#                             ncol = 6,
#                             byrow = TRUE)
#     graphics::layout(
#       mat = layout.matrix,
#       heights = c(1),
#       widths = c(6, 1, 1, 1, 1, 1)
#     )
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(4, 0, 2, 5))
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 4,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 4,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 0, 2, 0.5))
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 3],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Har. conf lvl",
#       ylab = "",
#       log = "xy",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 4],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Ar. conf lvl",
#       ylab = "",
#       log = "y",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 2],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "MTM power",
#       ylab = "",
#       log = "y",
#       ylim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         y = 1 / MTM_res_1[, 1],
#         x = MTM_res_1[, i],
#         ylim = c(min(y_axis), max(y_axis)),
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       x = average,
#       y = 2^eha_log2$log2_period,
#       log = "y",
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Wt. power",
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = add_abline_h)
#     }
#     par(new = FALSE,
#         mar = c(4, 4, 2, 0),
#         mgp = c(2, 1, 0))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_data == FALSE & add_avg == FALSE & add_MTM == TRUE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(5, 4, 3, 2, 1),
#                             nrow = 1,
#                             ncol = 5,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1),
#                      widths = c(7, 1, 1, 1, 1))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(4, 0, 2, 5))
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 4,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 4,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 0, 2, 0.5))
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 3],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Har. conf lvl",
#       ylab = "",
#       log = "xy",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 4],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "Ar. conf lvl",
#       ylab = "",
#       log = "y",
#       xlim = c(80, 101),
#       xaxs = "i",
#       ylim = ylim_vals
#     )
#     abline(
#       v = c(90, 95, 99),
#       lty = 3,
#       col = "grey",
#       lwd = 2
#     )
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     plot(
#       y = 1 / MTM_res_1[, 1],
#       x = MTM_res_1[, 2],
#       type = "l",
#       yaxs = "i",
#       yaxt = "n",
#       xlab = "MTM power",
#       ylab = "",
#       log = "y",
#       ylim = ylim_vals
#     )
#     for (i in 5:8) {
#       lines(
#         y = 1 / MTM_res_1[, 1],
#         x = MTM_res_1[, i],
#         ylim = ylim_vals,
#         lty = 3,
#         col = "grey",
#         lwd = 2
#       )
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = 1 / MTM_res_2[, 1],
#              col = "red",
#              lty = 3)
#     }
#     par(new = FALSE,
#         mar = c(4, 4, 2, 0),
#         mgp = c(2, 1, 0))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_data == FALSE & add_avg == FALSE & add_MTM == FALSE &
#       plot_horizontal == TRUE) {
#     layout.matrix <- matrix(c(2, 1),
#                             nrow = 1,
#                             ncol = 2,
#                             byrow = TRUE)
#     graphics::layout(mat = layout.matrix,
#                      heights = c(1),
#                      widths = c(10, 2.25))
#     power_max_mat.levels = quantile(mat, probs = seq(
#       from = 0,
#       to = 1,
#       length.out = n.levels + 1
#     ))
#     par(mar = c(4, 0, 2, 5))
#     image(
#       y = seq(from = 0, to = n.levels),
#       x = 1,
#       z = (matrix(power_max_mat.levels, nrow = 1)),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       yaxt = "n",
#       xaxt = "n",
#       xlab = "",
#       ylab = "",
#
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = key.marks,
#       labels = NA,
#       tck = 0.02,
#       tcl = 1.24
#     )
#     mtext(
#       key.labels,
#       side = 4,
#       at = key.marks,
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     mtext(
#       c("Power"),
#       side = 4,
#       at = mean(key.marks),
#       line = 0.5,
#       las = 2,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     box(lwd = lwd.axis)
#     par(new = FALSE, mar = c(4, 4, 2, 0.5))
#     image(
#       x = eha_log2$depth,
#       y = eha_log2$log2_period,
#       z =t (mat),
#       col = key.cols,
#       breaks = power_max_mat.levels,
#       useRaster = TRUE,
#       ylab = periodlab,
#       xlab = x_lab,
#       axes = TRUE,
#       yaxt = "n",
#       main = main,
#       xlim = xlim_vals,
#       ylim = log2(ylim_vals)
#     )
#     box(lwd = lwd.axis)
#     period.tick = unique(trunc(eha_log2$log2_period))
#     period.tick = na.omit(period.tick)
#     period.tick.label = 2^(period.tick)
#     axis(
#       2,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     axis(
#       4,
#       lwd = lwd.axis,
#       at = period.tick,
#       labels = NA,
#       tck = periodtck,
#       tcl = periodtcl
#     )
#     mtext(
#       period.tick.label,
#       side = 2,
#       at = period.tick,
#       las = 2,
#       line = par()$mgp[2] - 0.5,
#       font = par()$font.axis,
#       cex = par()$cex.axis
#     )
#     if (is.null(add_lines) != TRUE) {
#       for (i in 2:ncol(add_lines))
#         lines(add_lines[, 1], log2(add_lines[, i]))
#     }
#     if (is.null(add_points) != TRUE) {
#       for (i in 2:ncol(add_points))
#         points(add_points[, 1], log2(add_points[, i]))
#     }
#     if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) ==
#         FALSE) {
#       abline(h = log2(1 / MTM_res_2[, 1]),
#              col = "black",
#              lty = 3)
#     }
#     if (is.null(add_abline_h) != TRUE) {
#       abline(h = log2(add_abline_h))
#     }
#     if (is.null(add_abline_v) != TRUE) {
#       abline(v = add_abline_v)
#     }
#   }
#   if (add_MTM_peaks == TRUE & is.na(MTM_res_2[1, 2]) == FALSE) {
#     return(invisible(mtm_res))
#   }
# }

Try the WaverideR package in your browser

Any scripts or data that you put into this service are public.

WaverideR documentation built on April 6, 2026, 5:06 p.m.