"rad.cov" <-
function (x1, x2, p = 1, with.log = TRUE, with.constant = TRUE, 
            C = NA,marginal=FALSE) 
{
#
# mth order thin plate spline radial basis functions
# in d dimensions
# usually called with p 2m-d

#  marginal dummy argument 
    if( marginal){
    return( rep( 1, nrow( 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)
        m <- (d + p)/2
        par <- c(p/2, ifelse((d%%2 == 0) & (with.log), 1, 0))

# 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)), 
            PACKAGE="fields")

        temp<- matrix(temp$k, ncol = n2, nrow = n1)
    }

#
# do cross covariance matrix multiplication in FORTRAN
    else {

        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), h = as.double(rep(0, n1)), 
              work = as.double(rep(0, n2)), PACKAGE="fields")$h

    }

#
# multiply by constant if requested

        if (with.constant) {
            temp <- radbas.constant(m, d)*temp
        }

   return( temp)

}