# 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
"Rad.cov" <- function(x1, x2, p = 1, m=NA,with.log = TRUE, 
    with.constant = TRUE, C = NA, marginal = FALSE, derivative = 0) {
    #
    # mth order thin plate spline radial basis functions
    # in d dimensions
    # usually called with p 2m-d
    #  marginal dummy argument
    #  this should only be called within predict.se.Krig
    #  and provides the correct calculation. Because this is
    #  a generalized covariance the marginal variance is not really
    #  defined.
    #
    if (marginal) {
        return(rep(0, nrow(x1)))
    }
    #
    # coerce locations to matrices, if x2 is missing use x1
    if (!is.matrix(x1)) 
        x1 <- as.matrix(x1)
    if (!is.matrix(x2)) 
        x2 <- as.matrix(x2)
    d <- ncol(x1)
    n1 <- nrow(x1)
    n2 <- nrow(x2)
    if( is.na(m)){
      m <- (d + p)/2}
    else{
      p <- 2*m -d}
    if( p < 0 ) {
      stop(" p is negative (m possibly too small)")}
    # parameter list to send to the FORTRAN
      par <- c(p/2, ifelse((d%%2 == 0) & (with.log), 1, 0))
    #
    # multiply by constant if requested
      rbf.constant <- ifelse(with.constant, radbas.constant(m, 
        d), 1)
    # compute matrix in FORTRAN
    if (is.na(C[1])) {
        temp <- .Fortran("radbas", nd = as.integer(d), x1 = as.double(x1), 
            n1 = as.integer(n1), x2 = as.double(x2), n2 = as.integer(n2), 
            par = as.double(par), k = as.double(rep(0, n1 * n2)))
        return(rbf.constant * matrix(temp$k, ncol = n2, nrow = n1))
    }
    else {
        #   do cross covariance matrix multiplication in FORTRAN
        if (derivative == 0) {
            #     evaluate function not partial derivatives.
            C <- as.matrix(C)
            n3 <- ncol(C)
            temp <- .Fortran("multrb", nd = as.integer(d), x1 = as.double(x1), 
                n1 = as.integer(n1), x2 = as.double(x2), n2 = as.integer(n2), 
                par = as.double(par), c = as.double(C), n3 = as.integer(n3), 
                h = as.double(rep(0, n1 * n3)), work = as.double(rep(0, 
                  n2)))$h
            return(rbf.constant * matrix(temp, nrow = n1, ncol = n3))
        }
        else {
            if (ncol(C) > 1) {
                stop("Can only evaluate derivatives on one spline fit")
            }
            temp <- .Fortran("mltdrb", nd = as.integer(d), x1 = as.double(x1), 
                n1 = as.integer(n1), x2 = as.double(x2), n2 = as.integer(n2), 
                par = as.double(par), c = as.double(C), h = as.double(rep(0, 
                  n1 * d)), work = as.double(rep(0, n2)))$h
            return(rbf.constant * matrix(temp, nrow = n1, ncol = d))
        }
    }
    stop("should not get here!")
}