R/plot_RadialPlot.R

In Luminescence: Comprehensive Luminescence Dating Data Analysis

#' @title Function to create a Radial Plot
#'
#' @description A Galbraith's radial plot is produced on a logarithmic or a linear scale.
#'
#' @details Details and the theoretical background of the radial plot are given in the
#' cited literature. This function is based on an S script of Rex Galbraith. To
#' reduce the manual adjustments, the function has been rewritten. Thanks to
#' Rex Galbraith for useful comments on this function. \cr
#' Plotting can be disabled by adding the argument `plot = "FALSE"`, e.g.
#' to return only numeric plot output.
#'
#' Earlier versions of the Radial Plot in this package had the 2-sigma-bar
#' drawn onto the z-axis. However, this might have caused misunderstanding in
#' that the 2-sigma range may also refer to the z-scale, which it does not!
#' Rather it applies only to the x-y-coordinate system (standardised error vs.
#' precision). A spread in doses or ages must be drawn as lines originating at
#' zero precision (x0) and zero standardised estimate (y0). Such a range may be
#' drawn by adding lines to the radial plot ( `line`, `line.col`,
#' `line.label`, cf. examples).
#'
#' A statistic summary, i.e. a collection of statistic measures of
#' centrality and dispersion (and further measures) can be added by specifying
#' one or more of the following keywords:
#' - `"n"` (number of samples),
#' - `"mean"` (mean De value),
#' - `"mean.weighted"` (error-weighted mean),
#' - `"median"` (median of the De values),
#' - `"sdrel"` (relative standard deviation in percent),
#' - `"sdrel.weighted"` (error-weighted relative standard deviation in percent),
#' - `"sdabs"` (absolute standard deviation),
#' - `"sdabs.weighted"` (error-weighted absolute standard deviation),
#' - `"serel"` (relative standard error),
#' - `"serel.weighted"` (error-weighted relative standard error),
#' - `"seabs"` (absolute standard error),
#' - `"seabs.weighted"` (error-weighted absolute standard error),
#' - `"in.2s"` (percent of samples in 2-sigma range),
#' - `"kurtosis"` (kurtosis) and
#' - `"skewness"` (skewness).
#'
#' @param data [data.frame] or [RLum.Results-class] object (**required**):
#' for `data.frame` two columns: De (`data[,1]`) and De error (`data[,2]`).
#' To plot several data sets in one plot, the data sets must be provided as
#' `list`, e.g. `list(data.1, data.2)`.
#'
#' @param na.rm [logical] (*with default*):
#' excludes `NA` values from the data set prior to any further operations.
#'
#' @param log.z [logical] (*with default*):
#' Option to display the z-axis in logarithmic scale. Default is `TRUE`.
#'
#' @param central.value [numeric]:
#' User-defined central value, primarily used for horizontal centring
#' of the z-axis.
#'
#' @param centrality [character] or [numeric] (*with default*):
#' measure of centrality, used for automatically centring the plot and drawing
#' the central line. Can either be one out of
#' - `"mean"`,
#' - `"median"`,
#' - `"mean.weighted"` and
#' - `"median.weighted"` or a
#' - numeric value used for the standardisation.
#'
#' @param mtext [character]:
#' additional text below the plot title.
#'
#' @param summary [character] (*optional*):
#' add statistic measures of centrality and dispersion to the plot.
#' Can be one or more of several keywords. See details for available keywords.
#'
#' @param summary.pos [numeric] or [character] (*with default*):
#' optional position coordinates or keyword (e.g. `"topright"`)
#' for the statistical summary. Alternatively, the keyword `"sub"` may be
#' specified to place the summary below the plot header. However, this latter
#' option is only possible if `mtext` is not used.
#'
#' @param legend [character] vector (*optional*):
#' legend content to be added to the plot.
#'
#' @param legend.pos [numeric] or [character] (with
#' default): optional position coordinates or keyword (e.g. `"topright"`)
#' for the legend to be plotted.
#'
#' @param stats [character]: additional labels of statistically
#' important values in the plot. One or more out of the following:
#' - `"min"`,
#' - `"max"`,
#' - `"median"`.
#'
#' @param rug [logical]:
#' Option to add a rug to the z-scale, to indicate the location of individual values
#'
#' @param plot.ratio [numeric]:
#' User-defined plot area ratio (i.e. curvature of the z-axis). If omitted,
#' the default value (`4.5/5.5`) is used and modified automatically to optimise
#' the z-axis curvature. The parameter should be decreased when data points
#' are plotted outside the z-axis or when the z-axis gets too elliptic.
#'
#' @param bar.col [character] or [numeric] (*with default*):
#' colour of the bar showing the 2-sigma range around the central
#' value. To disable the bar, use `"none"`. Default is `"grey"`.
#'
#' @param y.ticks [logical]:
#' Option to hide y-axis labels. Useful for data with small scatter.
#'
#' @param grid.col [character] or [numeric] (*with default*):
#' colour of the grid lines (originating at `[0,0]` and stretching to
#' the z-scale). To disable grid lines, use `"none"`. Default is `"grey"`.
#'
#' @param line [numeric]:
#' numeric values of the additional lines to be added.
#'
#' @param line.col [character] or [numeric]:
#' colour of the additional lines.
#'
#' @param line.label [character]:
#' labels for the additional lines.
#'
#' @param output [logical]:
#' Optional output of numerical plot parameters. These can be useful to
#' reproduce similar plots. Default is `FALSE`.
#'
#' @param ... Further plot arguments to pass. `xlab` must be a vector of
#' length 2, specifying the upper and lower x-axes labels.
#'
#' @return Returns a plot object.
#'
#' @section Function version: 0.5.9
#'
#' @author
#' Michael Dietze, GFZ Potsdam (Germany)\cr
#' Sebastian Kreutzer, Institute of Geography, Heidelberg University (Germany)\cr
#' Based on a rewritten S script of Rex Galbraith, 2010
#'
#' @seealso [plot], [plot_KDE], [plot_Histogram], [plot_AbanicoPlot]
#'
#' @references
#' Galbraith, R.F., 1988. Graphical Display of Estimates Having
#' Differing Standard Errors. Technometrics, 30 (3), 271-281.
#'
#' Galbraith, R.F., 1990. The radial plot: Graphical assessment of spread in
#' ages. International Journal of Radiation Applications and Instrumentation.
#' Part D. Nuclear Tracks and Radiation Measurements, 17 (3), 207-214.
#'
#' Galbraith, R. & Green, P., 1990. Estimating the component ages in a finite
#' mixture. International Journal of Radiation Applications and
#' Instrumentation. Part D. Nuclear Tracks and Radiation Measurements, 17 (3)
#' 197-206.
#'
#' Galbraith, R.F. & Laslett, G.M., 1993. Statistical models for mixed fission
#' track ages. Nuclear Tracks And Radiation Measurements, 21 (4), 459-470.
#'
#' Galbraith, R.F., 1994. Some Applications of Radial Plots. Journal of the
#' American Statistical Association, 89 (428), 1232-1242.
#'
#' Galbraith, R.F., 2010. On plotting OSL equivalent doses. Ancient TL, 28 (1),
#' 1-10.
#'
#' Galbraith, R.F. & Roberts, R.G., 2012. Statistical aspects of equivalent
#' dose and error calculation and display in OSL dating: An overview and some
#' recommendations. Quaternary Geochronology, 11, 1-27.
#'
#' @examples
#'
#' ## load example data
#' data(ExampleData.DeValues, envir = environment())
#' ExampleData.DeValues <- Second2Gray(
#'   ExampleData.DeValues$BT998, c(0.0438,0.0019))
#'
#' ## plot the example data straightforward
#' plot_RadialPlot(data = ExampleData.DeValues)
#'
#' ## now with linear z-scale
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   log.z = FALSE)
#'
#' ## now with output of the plot parameters
#' plot1 <- plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   log.z = FALSE,
#'   output = TRUE)
#' plot1
#' plot1$zlim
#'
#' ## now with adjusted z-scale limits
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   log.z = FALSE,
#'   zlim = c(100, 200))
#'
#' ## now the two plots with serious but seasonally changing fun
#' #plot_RadialPlot(data = data.3, fun = TRUE)
#'
#' ## now with user-defined central value, in log-scale again
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   central.value = 150)
#'
#' ## now with a rug, indicating individual De values at the z-scale
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   rug = TRUE)
#'
#' ## now with legend, colour, different points and smaller scale
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   legend.text = "Sample 1",
#'   col = "tomato4",
#'   bar.col = "peachpuff",
#'   pch = "R",
#'   cex = 0.8)
#'
#' ## now without 2-sigma bar, y-axis, grid lines and central value line
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   bar.col = "none",
#'   grid.col = "none",
#'   y.ticks = FALSE,
#'   lwd = 0)
#'
#' ## now with user-defined axes labels
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   xlab = c("Data error (%)", "Data precision"),
#'   ylab = "Scatter",
#'   zlab = "Equivalent dose [Gy]")
#'
#' ## now with minimum, maximum and median value indicated
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   central.value = 150,
#'   stats = c("min", "max", "median"))
#'
#' ## now with a brief statistical summary
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   summary = c("n", "in.2s"))
#'
#' ## now with another statistical summary as subheader
#' plot_RadialPlot(
#'   data = ExampleData.DeValues,
#'   summary = c("mean.weighted", "median"),
#'   summary.pos = "sub")
#'
#' ## now the data set is split into sub-groups, one is manipulated
#' data.1 <- ExampleData.DeValues[1:15,]
#' data.2 <- ExampleData.DeValues[16:25,] * 1.3
#'
#' ## now a common dataset is created from the two subgroups
#' data.3 <- list(data.1, data.2)
#'
#' ## now the two data sets are plotted in one plot
#' plot_RadialPlot(data = data.3)
#'
#' ## now with some graphical modification
#' plot_RadialPlot(
#'   data = data.3,
#'   col = c("darkblue", "darkgreen"),
#'   bar.col = c("lightblue", "lightgreen"),
#'   pch = c(2, 6),
#'   summary = c("n", "in.2s"),
#'   summary.pos = "sub",
#'   legend = c("Sample 1", "Sample 2"))
#'
#' @md
#' @export
plot_RadialPlot <- function(
  data,
  na.rm = TRUE,
  log.z = TRUE,
  central.value,
  centrality = "mean.weighted",
  mtext,
  summary,
  summary.pos,
  legend,
  legend.pos,
  stats,
  rug = FALSE,
  plot.ratio,
  bar.col,
  y.ticks = TRUE,
  grid.col,
  line,
  line.col,
  line.label,
  output = FALSE,
  ...
) {
  ## Homogenise input data format
  if(is(data, "list") == FALSE) {data <- list(data)}

  ## Check input data
  for(i in 1:length(data)) {
    if(is(data[[i]], "RLum.Results") == FALSE &
         is(data[[i]], "data.frame") == FALSE) {
      stop(paste("[plot_RadialPlot] Error: Input data format is neither",
                 "'data.frame' nor 'RLum.Results'"), call. = FALSE)
    } else {
      if(is(data[[i]], "RLum.Results") == TRUE) {
        data[[i]] <- get_RLum(data[[i]], "data")
      }

      ##use only the first two columns
      data[[i]] <- data[[i]][,1:2]
    }
  }

  ## check data and parameter consistency--------------------------------------
  if(missing(stats) == TRUE) {stats <- numeric(0)}
  if(missing(summary) == TRUE) {
    summary <- c("n", "in.2s")
  }

  if(missing(summary.pos) == TRUE) {
    summary.pos <- "sub"
  }
  if(missing(bar.col) == TRUE) {
    bar.col <- rep("grey80", length(data))
  }

  if(missing(grid.col) == TRUE) {
    grid.col <- rep("grey70", length(data))
  }

  if(missing(summary) == TRUE) {
    summary <- NULL
  }

  if(missing(summary.pos) == TRUE) {
    summary.pos <- "topleft"
  }

  if(missing(mtext) == TRUE) {
    mtext <- ""
  }

  ## check z-axis log-option for grouped data sets
  if(is(data, "list") == TRUE & length(data) > 1 & log.z == FALSE) {
    warning(paste("Option 'log.z' is not set to 'TRUE' altough more than one",
                  "data set (group) is provided."))
  }

  ## optionally, remove NA-values
  if(na.rm == TRUE) {
    for(i in 1:length(data)) {
      data[[i]] <- na.exclude(data[[i]])
    }
  }

  ## create preliminary global data set
  De.global <- data[[1]][,1]
  if(length(data) > 1) {
    for(i in 2:length(data)) {
      De.global <- c(De.global, data[[i]][,1])
    }
  }

  ## calculate major preliminary tick values and tick difference
  extraArgs <- list(...)
  if("zlim" %in% names(extraArgs)) {
    limits.z <- extraArgs$zlim
  } else {
    z.span <- (mean(De.global) * 0.5) / (sd(De.global) * 100)
    z.span <- ifelse(z.span > 1, 0.9, z.span)
    limits.z <- c((ifelse(test = min(De.global) <= 0,
                          yes = 1.1,
                          no =  0.9) - z.span) * min(De.global),
                  (1.1 + z.span) * max(De.global))
  }

  ## calculate correction dose to shift negative values
  if(min(De.global) < 0 && log.z) {
    if("zlim" %in% names(extraArgs)) {
      De.add <- abs(extraArgs$zlim[1])
    } else {
      ## estimate delta De to add to all data
      De.add <-  min(10^ceiling(log10(abs(De.global))) * 10)

      ## optionally readjust delta De for extreme values
      if(De.add <= abs(min(De.global))) {
        De.add <- De.add * 10
      }
    }
  } else {
    De.add <- 0
  }

  ticks <- round(pretty(limits.z, n = 5), 3)
  De.delta <- ticks[2] - ticks[1]

  ## optionally add correction dose to data set and adjust error
  if(log.z) {
    for(i in 1:length(data))
      data[[i]][,1] <- data[[i]][,1] + De.add

    De.global <- De.global + De.add

  }

  ## calculate major preliminary tick values and tick difference
  extraArgs <- list(...)
  if("zlim" %in% names(extraArgs)) {
    limits.z <- extraArgs$zlim
  } else {
    z.span <- (mean(De.global) * 0.5) / (sd(De.global) * 100)
    z.span <- ifelse(z.span > 1, 0.9, z.span)
    limits.z <- c((ifelse(min(De.global) <= 0, 1.1, 0.9) - z.span) * min(De.global),
                  (1.1 + z.span) * max(De.global))

  }
  ticks <- round(pretty(limits.z, n = 5), 3)
  De.delta <- ticks[2] - ticks[1]

  ## calculate and append statistical measures --------------------------------
  ## z-values based on log-option
  z <- lapply(1:length(data), function(x){
    if(log.z == TRUE) {log(data[[x]][,1])} else {data[[x]][,1]}})

  if(is(z, "list") == FALSE) {z <- list(z)}
  data <- lapply(1:length(data), function(x) {
     cbind(data[[x]], z[[x]])})
  rm(z)

  ## calculate se-values based on log-option
  se <- lapply(1:length(data), function(x, De.add){
    if(log.z == TRUE) {
      if(De.add != 0) {
        data[[x]][,2] <- data[[x]][,2] / (data[[x]][,1] + De.add)
      } else {
        data[[x]][,2] / data[[x]][,1]
      }
    } else {
      data[[x]][,2]
    }}, De.add = De.add)

  if(is(se, "list") == FALSE) {se <- list(se)}
  data <- lapply(1:length(data), function(x) {
    cbind(data[[x]], se[[x]])})
  rm(se)

  ## calculate central values
  if(centrality[1] == "mean") {
    z.central <- lapply(1:length(data), function(x){
      rep(mean(data[[x]][,3], na.rm = TRUE), length(data[[x]][,3]))})
  } else if(centrality[1] == "median") {
    z.central <- lapply(1:length(data), function(x){
      rep(median(data[[x]][,3], na.rm = TRUE), length(data[[x]][,3]))})
  } else  if(centrality[1] == "mean.weighted") {
    z.central <- lapply(1:length(data), function(x){
      sum(data[[x]][,3] / data[[x]][,4]^2) /
        sum(1 / data[[x]][,4]^2)})
  } else if(centrality[1] == "median.weighted") {
    ## define function after isotone::weighted.median
    median.w <- function (y, w)
    {
      ox <- order(y)
      y <- y[ox]
      w <- w[ox]
      k <- 1
      low <- cumsum(c(0, w))
      up <- sum(w) - low
      df <- low - up
      repeat {
        if (df[k] < 0)
          k <- k + 1
        else if (df[k] == 0)
          return((w[k] * y[k] + w[k - 1] * y[k - 1]) / (w[k] + w[k - 1]))
        else return(y[k - 1])
      }
    }
    z.central <- lapply(1:length(data), function(x){
      rep(median.w(y = data[[x]][,3],
                   w = data[[x]][,4]), length(data[[x]][,3]))})
  } else if(is.numeric(centrality) & length(centrality) == length(data)) {
    z.central.raw <- if(log.z == TRUE) {
      log(centrality + De.add)
    } else {
      centrality + De.add
    }
    z.central <- lapply(1:length(data), function(x){
      rep(z.central.raw[x], length(data[[x]][,3]))})
  } else if(is.numeric(centrality) == TRUE &
              length(centrality) > length(data)) {
    z.central <- lapply(1:length(data), function(x){
      rep(median(data[[x]][,3], na.rm = TRUE), length(data[[x]][,3]))})
  } else {
    stop("[plot_RadialPlot()] Measure of centrality not supported!", call. = FALSE)
  }

  data <- lapply(1:length(data), function(x) {
    cbind(data[[x]], z.central[[x]])})
  rm(z.central)

  ## calculate precision
  precision <- lapply(1:length(data), function(x){
    1 / data[[x]][,4]})
  if(is(precision, "list") == FALSE) {precision <- list(precision)}
  data <- lapply(1:length(data), function(x) {
    cbind(data[[x]], precision[[x]])})
  rm(precision)

  ## calculate standard estimate
  std.estimate <- lapply(1:length(data), function(x){
    (data[[x]][,3] - data[[x]][,5]) / data[[x]][,4]})
  if(is(std.estimate, "list") == FALSE) {std.estimate <- list(std.estimate)}
  data <- lapply(1:length(data), function(x) {
    cbind(data[[x]], std.estimate[[x]])})

  ## append empty standard estimate for plotting
  data <- lapply(1:length(data), function(x) {
    cbind(data[[x]], std.estimate[[x]])})
  rm(std.estimate)

  ## generate global data set
  data.global <- cbind(data[[1]],
                       rep(x = 1,
                           times = nrow(data[[1]])))

  colnames(data.global) <- rep("", 9)

  if(length(data) > 1) {
    for(i in 2:length(data)) {
      data.add <- cbind(data[[i]],
                        rep(x = i, times = nrow(data[[i]])))
      colnames(data.add) <- rep("", 9)
      data.global <- rbind(data.global,
                           data.add)
    }
  }

  ## create column names
  colnames(data.global) <- c(
    "De", "error", "z", "se", "z.central", "precision", "std.estimate",
    "std.estimate.plot")

## calculate global central value
if(centrality[1] == "mean") {
  z.central.global <- mean(data.global[,3], na.rm = TRUE)
} else if(centrality[1] == "median") {
  z.central.global <- median(data.global[,3], na.rm = TRUE)
} else  if(centrality[1] == "mean.weighted") {
  z.central.global <- sum(data.global[,3] / data.global[,4]^2) /
    sum(1 / data.global[,4]^2)
} else if(centrality[1] == "median.weighted") {
  ## define function after isotone::weighted.mean
  median.w <- function (y, w)
  {
    ox <- order(y)
    y <- y[ox]
    w <- w[ox]
    k <- 1
    low <- cumsum(c(0, w))
    up <- sum(w) - low
    df <- low - up
    repeat {
      if (df[k] < 0)
        k <- k + 1
      else if (df[k] == 0)
        return((w[k] * y[k] + w[k - 1] * y[k - 1])/(w[k] +
                                                      w[k - 1]))
      else return(y[k - 1])
    }
  }
  z.central.global <- median.w(y = data.global[,3], w = data.global[,4])
} else if(is.numeric(centrality) == TRUE &
            length(centrality == length(data))) {
  z.central.global <- mean(data.global[,3], na.rm = TRUE)
}

  ## optionally adjust central value by user-defined value
  if(missing(central.value) == FALSE) {
    # ## adjust central value for De.add
    central.value <- central.value + De.add

    z.central.global <- ifelse(log.z == TRUE,
                               log(central.value),
                               central.value)
  }

  ## create column names
  for(i in 1:length(data)) {
    colnames(data[[i]]) <- c("De",
                             "error",
                             "z",
                             "se",
                             "z.central",
                             "precision",
                             "std.estimate",
                             "std.estimate.plot")
  }

  ## re-calculate standardised estimate for plotting
  for(i in 1:length(data)) {
    data[[i]][,8] <- (data[[i]][,3] - z.central.global) / data[[i]][,4]
  }

  data.global.plot <- data[[1]][,8]
  if(length(data) > 1) {
    for(i in 2:length(data)) {
      data.global.plot <- c(data.global.plot, data[[i]][,8])
    }
  }
  data.global[,8] <- data.global.plot

  ## print warning for too small scatter
  if(max(abs(1 / data.global[6])) < 0.02) {
    small.sigma <- TRUE
    message(paste("Attention, small standardised estimate scatter.",
                "Toggle off y.ticks?"))
}

  ## read out additional arguments---------------------------------------------
  extraArgs <- list(...)

  main <- if("main" %in% names(extraArgs)) {extraArgs$main} else
    {expression(paste(D[e], " distribution"))}

  sub <- if("sub" %in% names(extraArgs)) {extraArgs$sub} else {""}

  if("xlab" %in% names(extraArgs)) {
    if(length(extraArgs$xlab) != 2) {
      stop("Argmuent xlab is not of length 2!")
    } else {xlab <- extraArgs$xlab}
  } else {
    xlab <- c(if(log.z == TRUE) {
      "Relative standard error (%)"
      } else {
        "Standard error"
        },
      "Precision")
  }

  ylab <- if("ylab" %in% names(extraArgs)) {
    extraArgs$ylab
    } else {
      "Standardised estimate"
    }

  zlab <- if("zlab" %in% names(extraArgs)) {
    extraArgs$zlab
    } else {
      expression(paste(D[e], " [Gy]"))
    }

  if("zlim" %in% names(extraArgs)) {
    limits.z <- extraArgs$zlim
  } else {
    z.span <- (mean(data.global[,1]) * 0.5) / (sd(data.global[,1]) * 100)
    z.span <- ifelse(z.span > 1, 0.9, z.span)
    limits.z <- c((0.9 - z.span) * min(data.global[[1]]),
                  (1.1 + z.span) * max(data.global[[1]]))

  }

  if("xlim" %in% names(extraArgs)) {
    limits.x <- extraArgs$xlim
  } else {
    limits.x <- c(0, max(data.global[,6]))
  }

  if(limits.x[1] != 0) {
    limits.x[1] <- 0
    warning("Lower x-axis limit not set to zero, issue corrected!")
  }

  if("ylim" %in% names(extraArgs)) {
    limits.y <- extraArgs$ylim
  } else {
    y.span <- (mean(data.global[,1]) * 10) / (sd(data.global[,1]) * 100)
    y.span <- ifelse(y.span > 1, 0.98, y.span)
    limits.y <- c(-(1 + y.span) * max(abs(data.global[,7])),
                   (0.8 + y.span) * max(abs(data.global[,7])))
  }

  cex <- if("cex" %in% names(extraArgs)) {
    extraArgs$cex
  } else {
    1
  }

  lty <- if("lty" %in% names(extraArgs)) {
    extraArgs$lty
    } else {
      rep(2, length(data))
    }

  lwd <- if("lwd" %in% names(extraArgs)) {
    extraArgs$lwd
    } else {
      rep(1, length(data))
    }

  pch <- if("pch" %in% names(extraArgs)) {
    extraArgs$pch
    } else {
      rep(1, length(data))
    }

  col <- if("col" %in% names(extraArgs)) {
    extraArgs$col
    } else {
      1:length(data)
    }

  tck <- if("tck" %in% names(extraArgs)) {
    extraArgs$tck
  } else {
    NA
  }

  tcl <- if("tcl" %in% names(extraArgs)) {
    extraArgs$tcl
  } else {
    -0.5
  }

  show <- if("show" %in% names(extraArgs)) {extraArgs$show} else {TRUE}
  if(show != TRUE) {show <- FALSE}

  fun <- if("fun" %in% names(extraArgs)) {
    extraArgs$fun
  } else {
    FALSE
  }

  ## define auxiliary plot parameters -----------------------------------------

  ## optionally adjust plot ratio
  if(missing(plot.ratio)) {
    if(log.z) {
      plot.ratio <- 1 /  (1 * ((max(data.global[,6]) - min(data.global[,6])) /
        (max(data.global[,7]) - min(data.global[,7]))))
    } else {
      plot.ratio <- 4.5 / 5.5
    }
  }

  ##limit plot ratio
  plot.ratio <- min(c(1e+06, plot.ratio))

  ## calculate conversion factor for plot coordinates
  f <- (max(data.global[,6]) - min(data.global[,6])) /
       (max(data.global[,7]) - min(data.global[,7])) * plot.ratio

  ## calculate major and minor z-tick values
  tick.values.major <- signif(c(limits.z, pretty(limits.z, n = 5)))
  tick.values.minor <- signif(pretty(limits.z, n = 25), 3)

  tick.values.major <- tick.values.major[tick.values.major >=
    min(tick.values.minor)]
  tick.values.major <- tick.values.major[tick.values.major <=
    max(tick.values.minor)]
  tick.values.major <- tick.values.major[tick.values.major >=
    limits.z[1]]
  tick.values.major <- tick.values.major[tick.values.major <=
    limits.z[2]]
  tick.values.minor <- tick.values.minor[tick.values.minor >=
    limits.z[1]]
  tick.values.minor <- tick.values.minor[tick.values.minor <=
    limits.z[2]]

  if(log.z == TRUE) {
    tick.values.major <- log(tick.values.major)
    tick.values.minor <- log(tick.values.minor)
  }

  ## calculate z-axis radius
  r.x <- limits.x[2] / max(data.global[,6]) + 0.05
  r <- max(sqrt((data.global[,6])^2+(data.global[,7] * f)^2)) * r.x

  ## calculate major z-tick coordinates
  tick.x1.major <- r / sqrt(1 + f^2 * (
    tick.values.major - z.central.global)^2)

  tick.y1.major <- (tick.values.major - z.central.global) * tick.x1.major
  tick.x2.major <- (1 + 0.015 * cex) * r / sqrt(
    1 + f^2 * (tick.values.major - z.central.global)^2)
  tick.y2.major <- (tick.values.major - z.central.global) * tick.x2.major
  ticks.major <- cbind(0,
    tick.x1.major, tick.x2.major, tick.y1.major, tick.y2.major)


  ## calculate minor z-tick coordinates
  tick.x1.minor <- r / sqrt(1 + f^2 * (
    tick.values.minor - z.central.global)^2)
  tick.y1.minor <- (tick.values.minor - z.central.global) * tick.x1.minor
  tick.x2.minor <- (1 + 0.007 * cex) * r / sqrt(
    1 + f^2 * (tick.values.minor - z.central.global)^2)
  tick.y2.minor <- (tick.values.minor - z.central.global) * tick.x2.minor
  ticks.minor <- cbind(tick.x1.minor,
                       tick.x2.minor,
                       tick.y1.minor,
                       tick.y2.minor)

  ## calculate z-label positions
  label.x <- 1.03 * r / sqrt(1 + f^2 *
    (tick.values.major - z.central.global)^2)
  label.y <- (tick.values.major - z.central.global) * tick.x2.major

  ## create z-axes labels
  if(log.z) {
    label.z.text <- signif(exp(tick.values.major), 3)

  } else {
    label.z.text <- signif(tick.values.major, 3)

  }

  ## subtract De.add from label values
  if(De.add != 0)
    label.z.text <- label.z.text - De.add

  labels <- cbind(label.x, label.y, label.z.text)

  ## calculate coordinates for 2-sigma-polygon overlay
  polygons <- matrix(nrow = length(data), ncol = 8)

  for(i in 1:length(data)) {
    polygons[i,1:4] <- c(limits.x[1],
                         limits.x[1],
                         max(data.global[,6]),
                         max(data.global[,6]))
    polygons[i,5:8] <- c(-2,
                         2,
                         (data[[i]][1,5] - z.central.global) *
                           polygons[i,3] + 2,
                         (data[[i]][1,5] - z.central.global) *
                           polygons[i,4] - 2)
  }

  ## calculate node coordinates for semi-circle
  user.limits <- if(log.z) log(limits.z) else limits.z

  ellipse.values <- seq(
    from = min(c(tick.values.major, tick.values.minor, user.limits[1])),
    to = max(c(tick.values.major,tick.values.minor, user.limits[2])),
    length.out = 500)
  ellipse.x <- r / sqrt(1 + f^2 * (ellipse.values - z.central.global)^2)
  ellipse.y <- (ellipse.values - z.central.global) * ellipse.x
  ellipse <- cbind(ellipse.x, ellipse.y)
  ellipse.lims <- rbind(range(ellipse[,1]), range(ellipse[,2]))

  ## check if z-axis overlaps with 2s-polygon
  polygon_y_max <- max(polygons[,7])
  polygon_y_min <- min(polygons[,7])

  z_2s_upper <- ellipse.x[abs(ellipse.y - polygon_y_max) ==
                            min(abs(ellipse.y - polygon_y_max))]

  z_2s_lower <- ellipse.x[abs(ellipse.y - polygon_y_min) ==
                            min(abs(ellipse.y - polygon_y_min))]

  if(max(polygons[,3]) >= z_2s_upper | max(polygons[,3]) >= z_2s_lower) {
    warning("[plot_RadialPlot()] z-scale touches 2s-polygon. Decrease plot ratio.",
            call. = FALSE)
  }

  ## calculate statistical labels
  if(length(stats == 1)) {stats <- rep(stats, 2)}
  stats.data <- matrix(nrow = 3, ncol = 3)
  data.stats <- as.numeric(data.global[,1])

  if("min" %in% stats == TRUE) {
    stats.data[1, 3] <- data.stats[data.stats == min(data.stats)][1]
    stats.data[1, 1] <- data.global[data.stats == stats.data[1, 3], 6][1]
    stats.data[1, 2] <- data.global[data.stats == stats.data[1, 3], 8][1]
  }

  if("max" %in% stats == TRUE) {
    stats.data[2, 3] <- data.stats[data.stats == max(data.stats)][1]
    stats.data[2, 1] <- data.global[data.stats == stats.data[2, 3], 6][1]
    stats.data[2, 2] <- data.global[data.stats == stats.data[2, 3], 8][1]
  }

  if("median" %in% stats == TRUE) {
    stats.data[3, 3] <- data.stats[data.stats == quantile(data.stats, 0.5, type = 3)]
    stats.data[3, 1] <- data.global[data.stats == stats.data[3, 3], 6][1]
    stats.data[3, 2] <- data.global[data.stats == stats.data[3, 3], 8][1]
  }

  ## recalculate axes limits if necessary
  limits.z.x <- range(ellipse[,1])
  limits.z.y <- range(ellipse[,2])
  if(!("ylim" %in% names(extraArgs))) {
    if(limits.z.y[1] < 0.66 * limits.y[1]) {
      limits.y[1] <- 1.8 * limits.z.y[1]
    }
    if(limits.z.y[2] > 0.77 * limits.y[2]) {
      limits.y[2] <- 1.3 * limits.z.y[2]
    }
  }
  if(!("xlim" %in% names(extraArgs))) {
    if(limits.z.x[2] > 1.1 * limits.x[2]) {
      limits.x[2] <- limits.z.x[2]
    }
  }

  ## calculate and paste statistical summary
  De.stats <- matrix(nrow = length(data), ncol = 18)
  colnames(De.stats) <- c("n", "mean", "mean.weighted", "median", "median.weighted",
    "kde.max", "sd.abs", "sd.rel", "se.abs", "se.rel", "q25", "q75", "skewness",
    "kurtosis", "sd.abs.weighted", "sd.rel.weighted", "se.abs.weighted",
    "se.rel.weighted")

  for(i in 1:length(data)) {
    data_to_stats <- data[[i]][,1:2]

    ## remove added De
    if(log.z) data_to_stats$De <- data_to_stats$De - De.add

    statistics <- calc_Statistics(data = data_to_stats)
    De.stats[i,1] <- statistics$weighted$n
    De.stats[i,2] <- statistics$unweighted$mean
    De.stats[i,3] <- statistics$weighted$mean
    De.stats[i,4] <- statistics$unweighted$median
    De.stats[i,5] <- statistics$unweighted$median
    De.stats[i,7] <- statistics$unweighted$sd.abs
    De.stats[i,8] <- statistics$unweighted$sd.rel
    De.stats[i,9] <- statistics$unweighted$se.abs
    De.stats[i,10] <- statistics$weighted$se.rel
    De.stats[i,11] <- quantile(data[[i]][,1], 0.25)
    De.stats[i,12] <- quantile(data[[i]][,1], 0.75)
    De.stats[i,13] <- statistics$unweighted$skewness
    De.stats[i,14] <- statistics$unweighted$kurtosis
    De.stats[i,15] <- statistics$weighted$sd.abs
    De.stats[i,16] <- statistics$weighted$sd.rel
    De.stats[i,17] <- statistics$weighted$se.abs
    De.stats[i,18] <- statistics$weighted$se.rel

    ## kdemax - here a little doubled as it appears below again
    De.density <- try(density(x = data[[i]][,1],
                              kernel = "gaussian",
                              from = limits.z[1],
                              to = limits.z[2]),
                      silent = TRUE)

    if(!inherits(De.density, "try-error")) {
      De.stats[i,6] <- NA

    } else {
      De.stats[i,6] <- De.density$x[which.max(De.density$y)]
    }
  }

  label.text = list(NA)

  if(summary.pos[1] != "sub") {
    n.rows <- length(summary)

    for(i in 1:length(data)) {
      stops <- paste(rep("\n", (i - 1) * n.rows), collapse = "")

      summary.text <- character(0)

      for(j in 1:length(summary)) {
        summary.text <- c(summary.text,
                          paste(
                            "",
                            ifelse("n" %in% summary[j] == TRUE,
                                   paste("n = ",
                                         De.stats[i,1],
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("mean" %in% summary[j] == TRUE,
                                   paste("mean = ",
                                         round(De.stats[i,2], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("mean.weighted" %in% summary[j] == TRUE,
                                   paste("weighted mean = ",
                                         round(De.stats[i,3], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("median" %in% summary[j] == TRUE,
                                   paste("median = ",
                                         round(De.stats[i,4], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("median.weighted" %in% summary[j] == TRUE,
                                   paste("weighted median = ",
                                         round(De.stats[i,5], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("kdemax" %in% summary[j] == TRUE,
                                   paste("kdemax = ",
                                         round(De.stats[i,6], 2),
                                         " \n ",
                                         sep = ""),
                                   ""),
                            ifelse("sdabs" %in% summary[j] == TRUE,
                                   paste("sd = ",
                                         round(De.stats[i,7], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("sdrel" %in% summary[j] == TRUE,
                                   paste("rel. sd = ",
                                         round(De.stats[i,8], 2), " %",
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("seabs" %in% summary[j] == TRUE,
                                   paste("se = ",
                                         round(De.stats[i,9], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("serel" %in% summary[j] == TRUE,
                                   paste("rel. se = ",
                                         round(De.stats[i,10], 2), " %",
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("skewness" %in% summary[j] == TRUE,
                                   paste("skewness = ",
                                         round(De.stats[i,13], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("kurtosis" %in% summary[j] == TRUE,
                                   paste("kurtosis = ",
                                         round(De.stats[i,14], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("in.2s" %in% summary[j] == TRUE,
                                   paste("in 2 sigma = ",
                                         round(sum(data[[i]][,7] > -2 &
                                                     data[[i]][,7] < 2) /
                                                 nrow(data[[i]]) * 100 , 1),
                                         " %",
                                         sep = ""),
                                   ""),
                            ifelse("sdabs.weighted" %in% summary[j] == TRUE,
                                   paste("abs. weighted sd = ",
                                         round(De.stats[i,15], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("sdrel.weighted" %in% summary[j] == TRUE,
                                   paste("rel. weighted sd = ",
                                         round(De.stats[i,16], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("seabs.weighted" %in% summary[j] == TRUE,
                                   paste("abs. weighted se = ",
                                         round(De.stats[i,17], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            ifelse("serel.weighted" %in% summary[j] == TRUE,
                                   paste("rel. weighted se = ",
                                         round(De.stats[i,18], 2),
                                         "\n",
                                         sep = ""),
                                   ""),
                            sep = ""))
      }

      summary.text <- paste(summary.text, collapse = "")

      label.text[[length(label.text) + 1]] <- paste(stops,
                                                    summary.text,
                                                    stops,
                                                    sep = "")
    }
  } else {
    for(i in 1:length(data)) {

      summary.text <- character(0)

      for(j in 1:length(summary)) {
        summary.text <- c(summary.text,
                          ifelse("n" %in% summary[j] == TRUE,
                                 paste("n = ",
                                       De.stats[i,1],
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("mean" %in% summary[j] == TRUE,
                                 paste("mean = ",
                                       round(De.stats[i,2], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("mean.weighted" %in% summary[j] == TRUE,
                                 paste("weighted mean = ",
                                       round(De.stats[i,3], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("median" %in% summary[j] == TRUE,
                                 paste("median = ",
                                       round(De.stats[i,4], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("median.weighted" %in% summary[j] == TRUE,
                                 paste("weighted median = ",
                                       round(De.stats[i,5], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("kdemax" %in% summary[j] == TRUE,
                                 paste("kdemax = ",
                                       round(De.stats[i,6], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("sdrel" %in% summary[j] == TRUE,
                                 paste("rel. sd = ",
                                       round(De.stats[i,8], 2), " %",
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("sdabs" %in% summary[j] == TRUE,
                                 paste("abs. sd = ",
                                       round(De.stats[i,7], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("serel" %in% summary[j] == TRUE,
                                 paste("rel. se = ",
                                       round(De.stats[i,10], 2), " %",
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("seabs" %in% summary[j] == TRUE,
                                 paste("abs. se = ",
                                       round(De.stats[i,9], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("skewness" %in% summary[j] == TRUE,
                                 paste("skewness = ",
                                       round(De.stats[i,13], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("kurtosis" %in% summary[j] == TRUE,
                                 paste("kurtosis = ",
                                       round(De.stats[i,14], 2),
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("in.2s" %in% summary[j] == TRUE,
                                 paste("in 2 sigma = ",
                                       round(sum(data[[i]][,7] > -2 &
                                                   data[[i]][,7] < 2) /
                                               nrow(data[[i]]) * 100 , 1),
                                       " %   ",
                                       sep = ""),
                                 ""),
                          ifelse("sdabs.weighted" %in% summary[j] == TRUE,
                                 paste("abs. weighted sd = ",
                                       round(De.stats[i,15], 2), " %",
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("sdrel.weighted" %in% summary[j] == TRUE,
                                 paste("rel. weighted sd = ",
                                       round(De.stats[i,16], 2), " %",
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("seabs.weighted" %in% summary[j] == TRUE,
                                 paste("abs. weighted se = ",
                                       round(De.stats[i,17], 2), " %",
                                       " | ",
                                       sep = ""),
                                 ""),
                          ifelse("serel.weighted" %in% summary[j] == TRUE,
                                 paste("rel. weighted se = ",
                                       round(De.stats[i,18], 2), " %",
                                       " | ",
                                       sep = ""),
                                 "")
        )
      }

      summary.text <- paste(summary.text, collapse = "")

      label.text[[length(label.text) + 1]]  <- paste(
        "  ",
        summary.text,
        sep = "")
    }

    ## remove outer vertical lines from string
    for(i in 2:length(label.text)) {
      label.text[[i]] <- substr(x = label.text[[i]],
                                start = 3,
                                stop = nchar(label.text[[i]]) - 3)
    }
  }

## remove dummy list element
label.text[[1]] <- NULL
  ## convert keywords into summary placement coordinates
  if(missing(summary.pos) == TRUE) {
    summary.pos <- c(limits.x[1], limits.y[2])
    summary.adj <- c(0, 1)
  } else if(length(summary.pos) == 2) {
    summary.pos <- summary.pos
    summary.adj <- c(0, 1)
  } else if(summary.pos[1] == "topleft") {
    summary.pos <- c(limits.x[1], limits.y[2])
    summary.adj <- c(0, 1)
  } else if(summary.pos[1] == "top") {
    summary.pos <- c(mean(limits.x), limits.y[2])
    summary.adj <- c(0.5, 1)
  } else if(summary.pos[1] == "topright") {
    summary.pos <- c(limits.x[2], limits.y[2])
    summary.adj <- c(1, 1)
  }  else if(summary.pos[1] == "left") {
    summary.pos <- c(limits.x[1], mean(limits.y))
    summary.adj <- c(0, 0.5)
  } else if(summary.pos[1] == "center") {
    summary.pos <- c(mean(limits.x), mean(limits.y))
    summary.adj <- c(0.5, 0.5)
  } else if(summary.pos[1] == "right") {
    summary.pos <- c(limits.x[2], mean(limits.y))
    summary.adj <- c(1, 0.5)
  }else if(summary.pos[1] == "bottomleft") {
    summary.pos <- c(limits.x[1], limits.y[1])
    summary.adj <- c(0, 0)
  } else if(summary.pos[1] == "bottom") {
    summary.pos <- c(mean(limits.x), limits.y[1])
    summary.adj <- c(0.5, 0)
  } else if(summary.pos[1] == "bottomright") {
    summary.pos <- c(limits.x[2], limits.y[1])
    summary.adj <- c(1, 0)
  }

  ## convert keywords into legend placement coordinates
  if(missing(legend.pos) == TRUE) {
    legend.pos <- c(limits.x[1], limits.y[2])
    legend.adj <- c(0, 1)
  } else if(length(legend.pos) == 2) {
    legend.pos <- legend.pos
    legend.adj <- c(0, 1)
  } else if(legend.pos[1] == "topleft") {
    legend.pos <- c(limits.x[1], limits.y[2])
    legend.adj <- c(0, 1)
  } else if(legend.pos[1] == "top") {
    legend.pos <- c(mean(limits.x), limits.y[2])
    legend.adj <- c(0.5, 1)
  } else if(legend.pos[1] == "topright") {
    legend.pos <- c(limits.x[2], limits.y[2])
    legend.adj <- c(1, 1)
  } else if(legend.pos[1] == "left") {
    legend.pos <- c(limits.x[1], mean(limits.y))
    legend.adj <- c(0, 0.5)
  } else if(legend.pos[1] == "center") {
    legend.pos <- c(mean(limits.x), mean(limits.y))
    legend.adj <- c(0.5, 0.5)
  } else if(legend.pos[1] == "right") {
    legend.pos <- c(limits.x[2], mean(limits.y))
    legend.adj <- c(1, 0.5)
  } else if(legend.pos[1] == "bottomleft") {
    legend.pos <- c(limits.x[1], limits.y[1])
    legend.adj <- c(0, 0)
  } else if(legend.pos[1] == "bottom") {
    legend.pos <- c(mean(limits.x), limits.y[1])
    legend.adj <- c(0.5, 0)
  } else if(legend.pos[1] == "bottomright") {
    legend.pos <- c(limits.x[2], limits.y[1])
    legend.adj <- c(1, 0)
  }

  ## calculate line coordinates and further parameters
  if(!missing(line)) {
    #line = line + De.add

    if(log.z == TRUE) {line <- log(line)}

    line.coords <- list(NA)

    for(i in 1:length(line)) {
      line.x <- c(limits.x[1],
                  r / sqrt(1 + f^2 * (line[i] - z.central.global)^2))
      line.y <- c(0, (line[i] - z.central.global) * line.x[2])

      line.coords[[length(line.coords) + 1]] <- rbind(line.x, line.y)
    }

    line.coords[1] <- NULL

    if(missing(line.col) == TRUE) {
      line.col <- seq(from = 1, to = length(line.coords))
    }

    if(missing(line.label) == TRUE) {
      line.label <- rep("", length(line.coords))
    }
  }

  ## calculate rug coordinates
  if(missing(rug) == FALSE) {
    if(log.z == TRUE) {
      rug.values <- log(De.global)
    } else {
      rug.values <- De.global
    }

    rug.coords <- list(NA)

    for(i in 1:length(rug.values)) {
      rug.x <- c(r / sqrt(1 + f^2 * (rug.values[i] - z.central.global)^2) * 0.988,
                 r / sqrt(1 + f^2 * (rug.values[i] - z.central.global)^2) * 0.995)
      rug.y <- c((rug.values[i] - z.central.global) * rug.x[1],
                 (rug.values[i] - z.central.global) * rug.x[2])
      rug.coords[[length(rug.coords) + 1]] <- rbind(rug.x, rug.y)
    }

    rug.coords[1] <- NULL
  }

  ## Generate plot ------------------------------------------------------------

  ## check if plotting is enabled
  if(show) {

    ## determine number of subheader lines to shift the plot
    if(length(summary) > 0 & summary.pos[1] == "sub") {
      shift.lines <- length(data) + 1
    } else {shift.lines <- 1}

    ## setup plot area
    default <- par(mar = c(4, 4, shift.lines + 1.5, 7),
        xpd = TRUE,
        cex = cex)

    ## reset on exit
    on.exit(par(default))

    ## create empty plot
    plot(NA,
         xlim = limits.x,
         ylim = limits.y,
         main = "",
         sub = sub,
         xlab = "",
         ylab = "",
         xaxs = "i",
         yaxs = "i",
         frame.plot = FALSE,
         axes = FALSE)

    ## add y-axis label
    mtext(side = 2,
          line = 2.5,
          at = 0,
          adj = 0.5,
          cex = cex,
          text = ylab)

    ## calculate upper x-axis label values
    label.x.upper <- if(log.z == TRUE) {
      as.character(round(1/axTicks(side = 1)[-1] * 100, 1))
    } else {
      as.character(round(1/axTicks(side = 1)[-1], 1))
    }

    ## optionally, plot 2-sigma-bar
    if(bar.col[1] != "none") {
      for(i in 1:length(data)) {
        polygon(x = polygons[i,1:4],
                y = polygons[i,5:8],
                lty = "blank",
                col = bar.col[i])
      }
    }

    ## optionally, add grid lines
    if(grid.col[1] != "none") {
      for(i in 1:length(tick.x1.major)) {
        lines(x = c(limits.x[1], tick.x1.major[i]),
              y = c(0, tick.y1.major[i]),
              col = grid.col)
      }
    }

    ## optionally, plot central value lines
    if(lwd[1] > 0 & lty[1] > 0) {
      for(i in 1:length(data)) {
        x2 <- r / sqrt(1 + f^2 * (
          data[[i]][1,5] - z.central.global)^2)
        y2 <- (data[[i]][1,5] - z.central.global) * x2
        lines(x = c(limits.x[1], x2),
              y = c(0, y2),
              lty = lty[i],
              lwd = lwd[i],
              col = col[i])
      }
    }

    ## optionally add further lines
    if(missing(line) == FALSE) {
      for(i in 1:length(line)) {
        lines(x = line.coords[[i]][1,],
              y = line.coords[[i]][2,],
              col = line.col[i])
        text(x = line.coords[[i]][1,2],
             y = line.coords[[i]][2,2] + par()$cxy[2] * 0.3,
             labels = line.label[i],
             pos = 2,
             col = line.col[i],
             cex = cex * 0.9)
      }
    }

    ## overplot unwanted parts
    polygon(x = c(ellipse[,1], limits.x[2] * 2, limits.x[2] * 2),
            y = c(ellipse[,2], max(ellipse[,2]), min(ellipse[,2])),
            col = "white",
            lty = 0)

    ## add plot title
    title(main = main, line = shift.lines, font = 2)

    ## plot lower x-axis (precision)
    x.axis.ticks <- axTicks(side = 1)
    x.axis.ticks <- x.axis.ticks[c(TRUE, x.axis.ticks <= limits.x[2])]
    x.axis.ticks <- x.axis.ticks[x.axis.ticks <= limits.x[2]]

    ## axis with lables and ticks
    axis(side = 1,
         at = x.axis.ticks,
         lwd = 1,
         xlab = "")

    ## extend axis line to right side of the plot
    lines(x = c(max(x.axis.ticks, na.rm = TRUE), limits.x[2]),
          y = c(limits.y[1], limits.y[1]))

    ## draw closing tick on right hand side
    axis(side = 1, tcl = 0.5, lwd = 0, lwd.ticks = 1, at = limits.x[2],
         labels = FALSE)
    axis(side = 1, tcl = -0.5, lwd = 0, lwd.ticks = 1, at = limits.x[2],
         labels = FALSE)

    ## add upper axis label
    mtext(text = xlab[1],
          at = (limits.x[1] + limits.x[2]) / 2,
          side = 1,
          line = -3.5,
          cex = cex)

    ## add lower axis label
    mtext(text = xlab[2],
          at = (limits.x[1] + limits.x[2]) / 2,
          side = 1,
          line = 2.5,
          cex = cex)

    ## plot upper x-axis
    axis(side = 1,
         tcl = 0.5,
         lwd = 0,
         lwd.ticks = 1,
         at = x.axis.ticks[-1],
         labels = FALSE)

    ## remove first tick label (infinity)
    label.x.upper <- label.x.upper[1:(length(x.axis.ticks) - 1)]

    ## add tick labels
    axis(side = 1,
         lwd = 0,
         labels = label.x.upper,
         at = x.axis.ticks[-1],
         line = -3)

    ## plot minor z-ticks
    for(i in 1:length(tick.values.minor)) {
      lines(x = c(tick.x1.minor[i], tick.x2.minor[i]),
            y = c(tick.y1.minor[i], tick.y2.minor[i]))
    }

    ## plot major z-ticks
    for(i in 1:length(tick.values.major)) {
      lines(x = c(tick.x1.major[i], tick.x2.major[i]),
            y = c(tick.y1.major[i], tick.y2.major[i]))
    }

    ## plot z-axis
    lines(ellipse)

    ## plot z-values
    text(x = label.x,
         y = label.y,
         label = label.z.text, 0)

    ## plot z-label
    mtext(side = 4,
          at = 0,
          line = 5,
          las = 3,
          adj = 0.5,
          cex = cex,
          text = zlab)

    ## optionally add rug
    if(rug == TRUE) {
      for(i in 1:length(rug.coords)) {
        lines(x = rug.coords[[i]][1,],
              y = rug.coords[[i]][2,],
              col = col[data.global[i,9]])
      }
    }

    ## plot values
    for(i in 1:length(data)) {
      points(data[[i]][,6][data[[i]][,6] <= limits.x[2]],
             data[[i]][,8][data[[i]][,6] <= limits.x[2]],
             col = col[i],
             pch = pch[i])
    }

    ## optionally add min, max, median sample text
    if(length(stats) > 0) {
      text(x = stats.data[,1],
           y = stats.data[,2],
           labels = round(stats.data[,3], 1),
           pos = 2,
           cex = 0.85)
    }

    ## optionally add legend content
    if(missing(legend) == FALSE) {
      legend(x = legend.pos[1],
             y = 0.8 * legend.pos[2],
             xjust = legend.adj[1],
             yjust = legend.adj[2],
             legend = legend,
             pch = pch,
             col = col,
             text.col = col,
             cex = 0.8 * cex,
             bty = "n")
    }

    ## plot y-axis
    if(y.ticks == TRUE) {
      char.height <- par()$cxy[2]
      tick.space <- axisTicks(usr = limits.y, log = FALSE)
      tick.space <- (max(tick.space) - min(tick.space)) / length(tick.space)
      if(tick.space < char.height * 1.5) {
        axis(side = 2, at = c(-2, 2), labels = c("", ""), las = 1)
        axis(side = 2, at = 0, tcl = 0, labels = paste("\u00B1", "2"), las = 1)
      } else {
        axis(side = 2, at = seq(-2, 2, by = 2), las = 2)
      }
    } else {
      axis(side = 2, at = 0)
    }

    ## optionally add subheader text
    mtext(side = 3,
          line = shift.lines - 2,
          text = mtext,
          cex = 0.8 * cex)

    ## add summary content
    for(i in 1:length(data)) {
      if(summary.pos[1] != "sub") {
        text(x = summary.pos[1],
             y = 0.8 * summary.pos[2],
             adj = summary.adj,
             labels = label.text[[i]],
             cex = 0.8 * cex,
             col = col[i])
      } else {
        if(mtext == "") {
          mtext(side = 3,
                line = shift.lines - 1 - i,
                text = label.text[[i]],
                col = col[i],
                cex = 0.8 * cex)
        }
      }
    }

    ##FUN by R Luminescence Team
    if(fun==TRUE){sTeve()}
  }

  if(output) {
    return(list(data = data,
                data.global = data.global,
                xlim = limits.x,
                ylim = limits.y,
                zlim = limits.z,
                r = r,
                plot.ratio = plot.ratio,
                ticks.major = ticks.major,
                ticks.minor = ticks.minor,
                labels = labels,
                polygons = polygons,
                ellipse.lims = ellipse.lims))
  }

}