# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"stationary.cov" <- function(x1, x2, Covariance = "Exponential", 
    Distance = "rdist", Dist.args = NULL, theta = 1, V = NULL, 
    C = NA, marginal = FALSE, derivative=0, ...) {
    # get covariance function arguments from call
    Cov.args <- list(...)
    # coerce x1 and x2 to matrices
    if (is.data.frame(x1)) 
        x1 <- as.matrix(x1)
    if (!is.matrix(x1)) 
        x1 <- matrix(c(x1), ncol = 1)
    if (missing(x2)) 
        x2 <- x1
    if (is.data.frame(x2)) 
        x2 <- as.matrix(x2)
    if (!is.matrix(x2)) 
        x2 <- matrix(c(x2), ncol = 1)
    d <- ncol(x1)
    n1 <- nrow(x1)
    n2 <- nrow(x2)
    #
    # separate out a single scalar transformation and a
    # more complicated scaling and rotation.
    # this is done partly to have the use of great circle distance make sense
    # by applying the scaling  _after_ finding the distance.
    #
    # tranform coordinates if theta not a scalar
    if (length(theta) > 1) {
        stop("theta as a vector matrix has been depreciated use the V argument")
    }
    #
    # following now treats V as a full matrix for scaling and rotation.
    #
    # try to catch incorrect conbination  of great circle distance and V
    if (Distance == "rdist.earth" & !is.null(V)) {
        stop("V not supported with great circle distance")
    }
    if (!is.null(V)) {
        if (theta != 1) {
            stop("can't specify both theta and V!")
        }
        x1 <- x1 %*% t(solve(V))
        x2 <- x2 %*% t(solve(V))
    }
    #
    # locations are now scaled and rotated correctly
    # now apply covariance function to pairwise distance matrix, or multiply
    # by C vector or just find marginal variance
    # this if block finds the cross covariance matrix
    if (is.na(C[1]) & !marginal) {
        # note overall scalling by theta (which is just theta under isotropic case)
        bigD <- do.call(Distance, c(list(x1 = x1, x2 = x2), Dist.args))/theta
        return(do.call(Covariance, c(list(d = bigD), Cov.args)))
    }
    # or multiply cross covariance by C
    # as coded below this is not particularly efficient of memory
    if (!is.na(C[1])) {
      bigD <- do.call(Distance, c(list(x1 = x1, x2 = x2), Dist.args))
      if( derivative==0){
        return(do.call(Covariance, c(list(d = bigD/theta), Cov.args)) %*%C) }
      else{
        # find partial derivatives
        tempW <-do.call(Covariance, c(list(d = bigD/theta),
                                      Cov.args, derivative=derivative))
        # loop over dimensions and knock out each partial accumulate these in
        # in temp
        temp <- matrix(NA, ncol = d, nrow = n1)
        for (kd in 1:d) {
          # Be careful if the distance (tempD) is close to zero.
          # Note that the x1 and x2 are in transformed ( V inverse) scale
          sM <- ifelse(bigD == 0, 0,
            (tempW * (outer(x1[, kd] , x2[, kd], "-") )/(theta * bigD)))
          # accumlate the new partial
          temp[, kd] <- sM %*% C }
        # transform back to original coordinates.
        if (!is.null(V)) {
          temp <- temp %*% t(solve(V)) }
        return(temp) }
    }
    # or find marginal variance and return  a vector.
    if (marginal) {
        sigma2 <- do.call(Covariance, c(list(d = 0), Cov.args))
        return(rep(sigma2, nrow(x1)))
    }

# should not get here based on sequence of conditional if statements above.
}