#' Plot \code{biwavelet} objects
#'
#' Plot \code{biwavelet} objects such as the cwt, cross-wavelet and wavelet coherence
#'
#' @param x \code{biwavelet} object generated by \code{wt}, \code{xwt}, or
#'   \code{wtc}.
#' @param ncol number of colors to use. Default is 64.
#' @param fill.cols Vector of fill colors to be used. Users can specify color
#'   vectors using \code{colorRampPalette} or \code{brewer.pal} from package
#'   \code{RColorBrewer}. Default is \code{NULL} and will generate MATLAB's jet
#'   color palette.
#' @param xlab xlabel of the figure. Default is "Time"
#' @param ylab ylabel of the figure. Default is "Period"
#' @param tol tolerance level for significance contours. Sigificance contours
#'   will be drawn around all regions of the spectrum where
#'   \code{spectrum/percentile >= tol}. Default is \code{1}. If strict
#'   \code{i^{th}} percentile regions are desired, then
#'   \code{tol} must be set to \code{1}.
#' @param plot.cb plot color bar if TRUE. Default is FALSE.
#' @param plot.phase Plot phases with black arrows. Default is FALSE.
#' @param type type of plot to create. Can be \code{power} to plot the power,
#'   \code{power.corr} to plot the bias-corrected power, \code{power.norm} to
#'   plot the power normalized by the variance, \code{power.corr.norm} to plot
#'   the bias-corrected power normalized by the variance, \code{wavelet} to plot
#'   thewavelet coefficients, or \code{phase} to plot the phase. Default is
#'   \code{power.corr.norm}.
#' @param plot.coi plot cone of influence (COI) as a semi-transparent polygon if
#'   TRUE. Default is TRUE. Areas that fall within the polygon can be affected
#'   by edge effects.
#' @param lwd.coi Line width of COI. Default is 1.
#' @param col.coi Color of COI. Default is \code{white}.
#' @param lty.coi Line type of COI. Default is 1 for solide lines.
#' @param alpha.coi Transparency of COI. Range is 0 (full transparency) to 1 (no
#'   transparency). Default is 0.5.
#' @param plot.sig plot contours for significance if TRUE. Default is TRUE.
#' @param lwd.sig Line width of significance contours. Default is 4.
#' @param col.sig Color of significance contours. Default is \code{black}.
#' @param lty.sig Line type of significance contours. Default is 1.
#' @param bw plot in black and white if TRUE. Default is FALSE.
#' @param legend.loc legend location coordinates as defined by
#'   \code{image.plot}. Default is \code{NULL}.
#' @param legend.horiz plot a horizontal legend if TRUE. Default is FALSE.
#' @param arrow.len size of the arrows. Default is based on plotting region
#'   (min(par()$pin[2]/30,par()$pin[1]/40).
#' @param arrow.lwd width/thickness of arrows. Default is arrow.len*0.3.
#' @param arrow.cutoff cutoff value for plotting phase arrows. Phase arrows will be
#'        be plotted in regions where the significance of the zvalues exceeds \code{arrow.cutoff}.
#'        If the object being plotted does not have a significance field, regions
#'        whose zvalues exceed the \code{arrow.cutoff} quantile will be plotted. Default is 1.
#' @param arrow.col Color of arrows. Default is \code{black}.
#' @param xlim the x limits. The default is \code{NULL}.
#' @param ylim the y limits. The default is \code{NULL}.
#' @param zlim the z limits. The default is \code{NULL}.
#' @param xaxt Add x-axis? The default is \code{s}; use \code{n} for none.
#' @param yaxt Add y-axis? The default is \code{s}; use \code{n} for none.
#' @param form format to use to display dates on the x-axis. Default is '\%Y'
#'   for 4-digit year. See \code{?Date} for other valid formats.
#' @param \dots other parameters.
#'
#' @details
#' Arrows pointing to the right mean that \code{x} and \code{y} are in phase.
#'
#' Arrows pointing to the left mean that \code{x} and \code{y} are in anti-phase.
#'
#' Arrows pointing up mean that \code{x} leads \code{y} by \eqn{\pi/2}.
#'
#' Arrows pointing down mean that \code{y} leads \code{x} by \eqn{\pi/2}.
#'
#' @references
#' Cazelles, B., M. Chavez, D. Berteaux, F. Menard, J. O. Vik, S. Jenouvrier, and
#' N. C. Stenseth. 2008. Wavelet analysis of ecological time series.
#' \emph{Oecologia} 156:287-304.
#'
#' Grinsted, A., J. C. Moore, and S. Jevrejeva. 2004. Application of the cross
#' wavelet transform and wavelet coherence to geophysical time series.
#' \emph{Nonlinear Processes in Geophysics} 11:561-566.
#'
#' Torrence, C., and G. P. Compo. 1998. A Practical Guide to Wavelet Analysis.
#' \emph{Bulletin of the American Meteorological Society} 79:61-78.
#'
#' Liu, Y., X. San Liang, and R. H. Weisberg. 2007. Rectification of the Bias in
#' the Wavelet Power Spectrum. \emph{Journal of Atmospheric and Oceanic Technology}
#' 24:2093-2102.
#'
#' @author Tarik C. Gouhier (tarik.gouhier@@gmail.com)
#' Code based on WTC MATLAB package written by Aslak Grinsted.
#'
#' @seealso \code{\link{image.plot}}
#'
#' @example vignettes/example-plot.biwavelet.R
#' @importFrom fields image.plot
#' @importFrom grDevices adjustcolor colorRampPalette
#' @importFrom graphics axTicks axis box contour image par polygon
#' @export
plot.biwavelet <- function(x, ncol = 64, fill.cols = NULL,
                           xlab = "Time", ylab = "Period",
                           tol = 1, plot.cb = FALSE, plot.phase = FALSE,
                           type = "power.corr.norm",
                           plot.coi = TRUE, lwd.coi = 1, col.coi = "white",
                           lty.coi = 1, alpha.coi = 0.5, plot.sig = TRUE,
                           lwd.sig = 4, col.sig = "black", lty.sig = 1,
                           bw = FALSE,
                           legend.loc = NULL,
                           legend.horiz = FALSE,
                           arrow.len = min(par()$pin[2] / 30,
                                           par()$pin[1] / 40),
                           arrow.lwd = arrow.len * 0.3,
                           arrow.cutoff = 1,
                           arrow.col = "black",
                           xlim = NULL, ylim = NULL, zlim = NULL,
                           xaxt = "s", yaxt = "s", form = "%Y", ...) {

  if (is.null(fill.cols)) {
    if (bw) {
      fill.cols <- c("black", "white")
    } else {
      fill.cols <- c("#00007F", "blue", "#007FFF",
                     "cyan","#7FFF7F", "yellow",
                     "#FF7F00", "red", "#7F0000")
    }
  }

  col.pal <- colorRampPalette(fill.cols)
  fill.colors <- col.pal(ncol)

  types <- c("power.corr.norm", "power.corr", "power.norm",
             "power", "wavelet", "phase")

  type <- match.arg(tolower(type), types)

  if (type == "power.corr" | type == "power.corr.norm") {
    if (x$type == "wtc" | x$type == "xwt") {
      x$power <- x$power.corr
      x$wave <- x$wave.corr
    } else {
      x$power <- x$power.corr
    }
  }

  if (type == "power.norm" | type == "power.corr.norm") {
    if (x$type == "xwt") {
      zvals <- log2(x$power) / (x$d1.sigma * x$d2.sigma)

      if (is.null(zlim)) {
        zlim <- range(c(-1, 1) * max(zvals))
      }

      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- 2 ^ locs

    } else if (x$type == "wtc" | x$type == "pwtc") {
      zvals <- x$rsq
      zvals[!is.finite(zvals)] <- NA
      if (is.null(zlim)) {
        zlim <- range(zvals, na.rm = TRUE)
      }
      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- locs
    } else {
      zvals <- log2(abs(x$power / x$sigma2))
      if (is.null(zlim)) {
        zlim <- range(c(-1, 1) * max(zvals))
      }
      zvals[zvals < zlim[1]] <- zlim[1]
      locs <- pretty(range(zlim), n = 5)
      leg.lab <- 2 ^ locs
    }
  } else if (type == "power" | type == "power.corr") {
    zvals <- log2(x$power)
    if (is.null(zlim)) {
      zlim <- range( c(-1, 1) * max(zvals) )
    }
    zvals[zvals < zlim[1]] <- zlim[1]
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- 2 ^ locs
  } else if (type == "wavelet") {
    zvals <- (Re(x$wave))
    if (is.null(zlim)) {
      zlim <- range(zvals)
    }
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- locs
  } else if (type == "phase") {
    zvals <- x$phase
    if (is.null(zlim)) {
      zlim <- c(-pi, pi)
    }
    locs <- pretty(range(zlim), n = 5)
    leg.lab <- locs
  }

  if (is.null(xlim)) {
    xlim <- range(x$t)
  }

  yvals <- log2(x$period)

  if (is.null(ylim)) {
    ylim <- range(yvals)
  } else {
    ylim <- log2(ylim)
  }

  image(x$t,
        yvals,
        t(zvals),
        zlim = zlim,
        xlim = xlim,
        ylim = rev(ylim),
        xlab = xlab,
        ylab = ylab,
        yaxt = "n",
        xaxt = "n",
        col = fill.colors, ...)

  box()
  if (class(x$xaxis)[1] == "Date" | class(x$xaxis)[1] == "POSIXct") {
    if (xaxt != "n") {
      xlocs <- pretty(x$t) + 1
      axis(side = 1, at = xlocs, labels = format(x$xaxis[xlocs], form, ...))
    }
  } else {
    if (xaxt != "n") {
      xlocs <- axTicks(1)
      axis(side = 1, at = xlocs, ...)
    }
  }

  if (yaxt != "n") {
    axis.locs <- axTicks(2)
    yticklab <- format(2 ^ axis.locs, dig = 1)
    axis(2, at = axis.locs, labels = yticklab, ...)
  }

  # COI
  if (plot.coi) {
    polygon(x = c(x$t, rev(x$t)), lty = lty.coi, lwd = lwd.coi,
            y = c(log2(x$coi),
                  rep(max(log2(x$coi), na.rm = TRUE), length(x$coi))),
            col = adjustcolor(col.coi, alpha.f = alpha.coi), border = col.coi)
  }

  # sig.level contour (default is 95%)
  if (plot.sig & length(x$signif) > 1) {
    if (x$type %in% c("wt", "xwt")) {
      contour(x$t, yvals, t(x$signif), level = tol, col = col.sig,
              lwd = lwd.sig, add = TRUE, drawlabels = FALSE)
    } else {
      tmp <- x$rsq / x$signif
      contour(x$t, yvals, t(tmp), level = tol, col = col.sig, lwd = lwd.sig,
              add = TRUE, drawlabels = FALSE)
    }
  }

  # Plot phases
  if (plot.phase) {
    a <- x$phase

    # Remove phases where power is weak
    if (!is.null(x$type)) {
      if (x$type %in% c("wt", "xwt")) {
        locs.phases <- which(x$signif <= arrow.cutoff)
      } else {
        v <- x$rsq / x$signif
        locs.phases <- which(v <= arrow.cutoff)
      }
    } else {
      locs.phases <- which(zvals < quantile(zvals, arrow.cutoff))
    }
    a[locs.phases] <- NA

    phase.plot(x$t, log2(x$period), a,
               arrow.len = arrow.len, arrow.lwd = arrow.lwd,
               arrow.col = arrow.col)
  }
  box()

  ## Add color bar: this must happen after everything, otherwise chaos ensues!
  if (plot.cb) {
    image.plot(x$t,
               yvals,
               t(zvals),
               zlim = zlim,
               ylim = rev(range(yvals)),
               xlab = xlab,
               ylab = ylab,
               col = fill.colors,
               smallplot = legend.loc,
               horizontal = legend.horiz,
               legend.only = TRUE,
               axis.args =
                 list(at = locs, labels = format(leg.lab, dig = 2), ...),
               xpd = NA)
  }
}