## Make RK for thin-plate splines
mkrk.tp <- function(dm,order,mesh,weight=1)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    if (xor(is.vector(mesh),dm==1)
        |xor(is.matrix(mesh),dm>=2)) {
        stop("gss error in mkrk.tp: mismatched inputs")
    }
    if ((min(weight)<0)|(max(weight)<=0)) {
        stop("gss error in mkrk.tp: negative weights")
    }
    ## Set weights
    if (is.vector(mesh)) N <- length(mesh)
    else N <- dim(mesh)[1]
    weight <- rep(weight,len=N)
    weight <- sqrt(weight/sum(weight))
    ## Obtain orthonormal basis
    phi.p <- mkphi.tp.p(dm,order)
    nnull <- choose(dm+order-1,dm)
    s <- NULL
    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
    s <- qr(weight*s)
    if (s$rank<nnull) {
        stop("gss error in mkrk.tp: insufficient normalizing mesh for thin-plate spline")
    }
    q <- qr.Q(s)
    r <- qr.R(s)
    ## Set Q^{T}E(|u_{i}-u_{j}|)Q
    rk.p <- mkrk.tp.p(dm,order)
    pep <- weight*t(weight*rk.p$fun(mesh,mesh,rk.p$env,out=TRUE))
    pep <- t(q)%*%pep%*%q
    ## Create the environment
    env <- list(dim=dm,order=order,weight=weight,
                phi.p=phi.p,rk.p=rk.p,q=q,r=r,mesh=mesh,pep=pep)
    ## Create the rk function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ## Check inputs
        if (env$dim==1) {
            if (!(is.vector(x)&is.vector(y))) {
                stop("gss error in rk: inputs are of wrong types")
            }
            nx <- length(x)
            ny <- length(y)
        }
        else {
            if (is.vector(x)) x <- t(as.matrix(x))
            if (env$dim!=dim(x)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            nx <- dim(x)[1]
            if (is.vector(y)) y <- t(as.matrix(y))
            if (env$dim!=dim(y)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            ny <- dim(y)[1]
        }
        ## Return the results
        nnull <- choose(env$dim+env$order-1,env$dim)
        if (outer.prod) {
            phix <- phiy <- NULL
            for (nu in 1:nnull) {
                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
            }
            phix <- backsolve(env$r,phix,transpose=TRUE)
            phiy <- backsolve(env$r,phiy,transpose=TRUE)
            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
            ex <- env$weight*ex
            ex <- t(env$q)%*%ex
            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
            ey <- env$weight*ey
            ey <- t(env$q)%*%ey
            env$rk.p$fun(x,y,env$rk.p$env,out=TRUE)-t(phix)%*%ey-
                t(ex)%*%phiy+t(phix)%*%env$pep%*%phiy
        }
        else {
            N <- max(nx,ny)
            phix <- phiy <- NULL
            for (nu in 1:nnull) {
                phix <- rbind(phix,env$phi.p$fun(x,nu,env$phi.p$env))
                phiy <- rbind(phiy,env$phi.p$fun(y,nu,env$phi.p$env))
            }
            phix <- backsolve(env$r,phix,transpose=TRUE)
            phix <- matrix(phix,nnull,N)
            phiy <- backsolve(env$r,phiy,transpose=TRUE)
            phiy <- matrix(phiy,nnull,N)
            ex <- env$rk.p$fun(env$mesh,x,env$rk.p$env,out=TRUE)
            ex <- env$weight*ex
            ex <- t(env$q)%*%ex
            ex <- matrix(ex,nnull,N)
            ey <- env$rk.p$fun(env$mesh,y,env$rk.p$env,out=TRUE)
            ey <- env$weight*ey
            ey <- t(env$q)%*%ey
            ey <- matrix(ey,nnull,N)
            fn1 <- function(x,n) x[1:n]%*%x[n+(1:n)]
            fn2 <- function(x,pep,n) t(x[1:n])%*%pep%*%x[n+(1:n)]
            env$rk.p$fun(x,y,env$rk.p$env)-apply(rbind(phix,ey),2,fn1,nnull)-
                apply(rbind(phiy,ex),2,fn1,nnull)+
                    apply(rbind(phix,phiy),2,fn2,env$pep,nnull)
        }
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make phi function for thin-plate splines
mkphi.tp <-  function(dm,order,mesh,weight)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    if (xor(is.vector(mesh),dm==1)
        |xor(is.matrix(mesh),dm>=2)) {
        stop("gss error in mkphi.tp: mismatched inputs")
    }
    if ((min(weight)<0)|(max(weight)<=0)) {
        stop("gss error in mkphi.tp: negative weights")
    }
    ## Set weights
    if (is.vector(mesh)) N <- length(mesh)
    else N <- dim(mesh)[1]
    weight <- rep(weight,len=N)
    weight <- sqrt(weight/sum(weight))
    ## Create the environment
    phi.p <- mkphi.tp.p(dm,order)
    nnull <- choose(dm+order-1,dm)
    s <- NULL
    for (nu in 1:nnull) s <- cbind(s,phi.p$fun(mesh,nu,phi.p$env))
    s <- qr(weight*s)
    if (s$rank<nnull) {
        stop("gss error in mkphi: insufficient normalizing mesh for thin-plate spline")
    }
    r <- qr.R(s)
    env <- list(dim=dm,order=order,phi.p=phi.p,r=r)
    ## Create the phi function
    fun <- function(x,nu,env) {
        nnull <- choose(env$dim+env$order-1,env$dim)
        phix <- NULL
        for(i in 1:nnull)
            phix <- rbind(phix,env$phi.p$fun(x,i,env$phi.p$env))
        t(backsolve(env$r,phix,transpose=TRUE))[,nu+1]
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make pseudo RK for thin-plate splines
mkrk.tp.p <- function(dm,order)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    ## Create the environment
    if (dm%%2) {                    
        theta <- gamma(dm/2-order)/2^(2*order)/pi^(dm/2)/gamma(order)
    }
    else {
        theta <- (-1)^(dm/2+order+1)/2^(2*order-1)/pi^(dm/2)/
            gamma(order)/gamma(order-dm/2+1)
    }
    env <- list(dim=dm,order=order,theta=theta)
    ## Create the rk.p function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ## Check inputs
        if (env$dim==1) {
            if (!(is.vector(x)&is.vector(y))) {
                stop("gss error in rk: inputs are of wrong types")
            }
        }
        else {
            if (is.vector(x)) x <- t(as.matrix(x))
            if (env$dim!=dim(x)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
            if (is.vector(y)) y <- t(as.matrix(y))
            if (env$dim!=dim(y)[2]) {
                stop("gss error in rk: inputs are of wrong dimensions")
            }
        }
        ## Return the results
        if (outer.prod) {               
            if (env$dim==1) {
                fn1 <- function(x,y) abs(x-y)
                d <- outer(x,y,fn1)
            }
            else {
                fn2 <- function(x,y) sqrt(sum((x-y)^2))
                d <- NULL
                for (i in 1:dim(y)[1]) d <- cbind(d,apply(x,1,fn2,y[i,]))
            }
        }
        else {
            if (env$dim==1) d <- abs(x-y)
            else {
                N <- max(dim(x)[1],dim(y)[1])
                x <- t(matrix(t(x),env$dim,N))
                y <- t(matrix(t(y),env$dim,N))
                fn <- function(x) sqrt(sum(x^2))
                d <- apply(x-y,1,fn)
            }
        }
        power <- 2*env$order-env$dim
        switch(1+env$dim%%2,
               env$theta*d^power*log(ifelse(d>0,d,1)),
               env$theta*d^power)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make pseudo phi function for thin-plate splines
mkphi.tp.p <- function(dm,order)
{
    ## Check inputs
    if (!((2*order>dm)&(dm>=1))) {
        stop("gss error: thin-plate spline undefined for the parameters")
    }
    ## Create the environment
    pol.code <- NULL
    for (i in 0:(order^dm-1)) {
        ind <- i; code <- NULL
        for (j in 1:dm) {
            code <- c(code,ind%%order)
            ind <- ind%/%order
        }
        if (sum(code)<order) pol.code <- cbind(pol.code,code)
    }
    env <- list(dim=dm,pol.code=pol.code)
    ## Create the phi function  
    fun <- function(x,nu,env) {
        if (env$dim==1) x <- as.matrix(x)
        if (env$dim!=dim(x)[2]) {
            stop("gss error in phi: inputs are of wrong dimensions")
        }
        apply(t(x)^env$pol.code[,nu],2,prod)
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}

## Make RK for spherical splines
mkrk.sphere <- function(order)
{
    ## Create the environment
    env <- list(order=order)
    ## Create the rk function
    fun <- function(x,y,env,outer.prod=FALSE) {
        ##% Check the inputs
        if (is.vector(x)) x <- t(as.matrix(x))
        if (is.vector(y)) y <- t(as.matrix(y))
        if (!(is.matrix(x)&is.matrix(y))) {
            stop("gss error in rk: inputs are of wrong types")
        }
        if ((dim(x)[2]!=2)|(dim(y)[2]!=2)) {
            stop("gss error in rk: inputs are of wrong dimensions")
        }
        if ((max(abs(x[,1]),abs(y[,1]))>90)|(max(abs(x[,2]),abs(y[,2]))>180)) {
            stop("gss error in rk: inputs are out of range")
        }
        ##% Convert to radian
        lat.x <- x[,1]/180*pi; lon.x <- x[,2]/180*pi
        lat.y <- y[,1]/180*pi; lon.y <- y[,2]/180*pi
        ##% Return the result
        rk <- function(lat.x,lon.x,lat.y,lon.y,order) {
            z <- cos(lat.x)*cos(lat.y)*cos(lon.x-lon.y)+sin(lat.x)*sin(lat.y)
            W <- ifelse(z<1-10^(-10),(1-z)/2,0)
            A <- ifelse(W>0,log(1+1/sqrt(W)),0)
            C <- ifelse(W>0,2*sqrt(W),0)
            switch(order-1,
                   (A*4*W*(3*W-1)+6*W*(1-C)+1)/2,
                   (W*W*(A*((840*W-720)*W+72)+420*W*(1-C)+220*C-150)-4*W+3)/12,
                   (W*W*W*(A*(((27720*W-37800)*W+12600)*W-600)+
                           (13860*(1-C)*W+14280*C-11970)*W-2772*C+1470)+
                    15*W*W-3*W+5)/30) - 1/(2*order-1)
        }
        if (outer.prod) {
            zz <- NULL
            for (i in 1:length(lat.y))
                zz <- cbind(zz,rk(lat.x,lon.x,lat.y[i],lon.y[i],env$order))
        }
        else zz <- rk(lat.x,lon.x,lat.y,lon.y,env$order)
        zz
    }
    ## Return the function and the environment
    list(fun=fun,env=env)
}