
"spfmat"<-
function(x, y, np)
{
	n <- length(x)
	npar <- ((np + 1) * (np + 2))/2
	if(!is.loaded(symbol.C("VR_fmat"))) stop("Compiled code has not been dynamically loaded")
	z <- .C("VR_fmat",
		f = double(n * npar),
		as.double(x),
		as.double(y),
		as.integer(n),
		as.integer(np))
	z$f
}
"surf.ls"<-
function(np, x, y, z)
{
	if(!is.loaded(symbol.C("VR_frset"))) stop("Compiled code has not been dynamically loaded")
	if (np > 6) stop("np exceeds 6")
	if(is.data.frame(x)) {
		if(missing(y)) y <- x$y
		if(missing(z)) z <- x$z
		x <- x$x
	}
	rx <- range(x)
	ry <- range(y)
	.C("VR_frset",
		as.double(rx[1]),
		as.double(rx[2]),
		as.double(ry[1]),
		as.double(ry[2])
		)
	n <- length(x)
	npar <- ((np + 1) * (np + 2))/2
	f <- spfmat(x, y, np)
	Z <- .C("VR_ls",
		as.double(x),
		as.double(y),
		as.double(z),
		as.integer(n),
		as.integer(np),
		as.integer(npar),
		f = f,
		r = double((npar * (npar + 1))/2),
		beta = double(npar),
		wz = double(n),
		ifail = as.integer(0))
	structure(
	list(x=x, y=y, z=z, np = np, f=f, r = Z$r, beta=Z$beta, 
	wz=Z$wz, rx=rx, ry=ry, call=match.call()),
	class="trls")
}
"surf.gls"<-
function(np, covmod, x, y, z, nx=1000, ...)
{
	# if(nx > 10000) stop("nx is more than 10000")
	if (np > 6) stop("np exceeds 6")
	if(!is.loaded(symbol.C("VR_frset"))) stop("Compiled code has not been dynamically loaded")
	if(is.data.frame(x)) {
		if(missing(y)) y <- x$y
		if(missing(z)) z <- x$z
		x <- x$x
	}
	rx <- range(x)
	ry <- range(y)
	.C("VR_frset",
		as.double(rx[1]),
		as.double(rx[2]),
		as.double(ry[1]),
		as.double(ry[2])
		)
        covmod <- covmod
        arguments <- list(...)
        # covmod.argnames<-names(covmod)   should give the names of arguments
        # used in covmod(...), but it doesn't, so we have to do this "manually".
        # The following strips everthing before and after the parentheses of
        # the function declaration and splits its interior at ","s: 
        covmod.arg.tmp<-strsplit(strsplit(strsplit(deparse(args(covmod))[1],
                                                   "(")[[1]][2],
                                          ")")[[1]][1],
                                 ",")[[1]]
        # remove leading blanks, omit "=" and default values:
        f1<-function(x)strsplit(x," ")[[1]][1]
        covmod.argnames<-sapply(covmod.arg.tmp,f1)
        argnames.given<-names(arguments)
        if (length(arguments)) {
           pm <- pmatch(argnames.given, covmod.argnames, nomatch = 0)
           if (any(pm == 0)) 
           warning(paste("some of ... do not match"))
           # this doesn't work, use "..." later:
           # covmod[pm] <- unlist(arguments)
        }
        mm <- 1.5 * sqrt((rx[2] - rx[1])^2 + (ry[2] - ry[1])^2)
        # alph <- c(mm/nx, covmod(seq(0, mm, mm/nx)))
        alph <- c(mm/nx, covmod(seq(0, mm, mm/nx),...))
	.C("VR_alset",
		as.double(alph),
		as.integer(length(alph))
		)
	n <- length(x)
	npar <- ((np + 1) * (np + 2))/2
	f <- spfmat(x, y, np)
	Z <- .C("VR_gls",
		as.double(x),
		as.double(y),
		as.double(z),
		as.integer(n),
		as.integer(np),
		as.integer(npar),
		as.double(f),
		l = double((n * (n + 1))/2),
		r = double((npar * (npar + 1))/2),
		beta = double(npar),
		wz = double(n),
		yy = double(n),
		W = double(n),
		ifail = as.integer(0),
		l1f = double(n * npar)
		)
	if(Z$ifail > 0) stop("Rank failure in Choleski decomposition")
        if (nx > 1000) 
                alph <- alph[1]
	structure(
	list(x=x, y=y, z=z, np = np, f=f, alph=alph, l=Z$l, r = Z$r, 
	  beta=Z$beta, wz=Z$wz, yy=Z$yy, W=Z$W, l1f=Z$l1f, rx=rx, ry=ry,
	  covmod = covmod, call=match.call()), class=c("trgls", "trls"))
}
"trmat"<-
function(obj, xl, xu, yl, yu, n)
{
	if(!inherits(obj, "trls")) stop("object not a fitted trend surface")
	.C("VR_frset",
		as.double(obj$rx[1]),
		as.double(obj$rx[2]),
		as.double(obj$ry[1]),
		as.double(obj$ry[2])
		)
	dx <- (xu - xl)/n
	dy <- (yu - yl)/n
	x <- seq(xl, xu, dx)
	y <- seq(yl, yu, dy)
	z <- matrix(nrow = length(x), ncol = length(y))
	for(i in 1:length(y))
		z[, i] <- trval(obj, x, rep(y[i], length(x)))
	invisible(list(x = x, y = y, z = z))
}

"trval"<-
function(obj, x, y)
{
	n <- length(x)
	z <- .C("VR_valn",
		z = double(n),
		as.double(x),
		as.double(y),
		as.integer(n),
		as.double(obj$beta),
		as.integer(obj$np))
	z$z
}
"prmat"<-
function(obj, xl, xu, yl, yu, n)
{
	if(!inherits(obj, "trgls")) stop("object not from kriging")
	if(n > 999) stop("n is too large")
	.C("VR_frset",
		as.double(obj$rx[1]),
		as.double(obj$rx[2]),
		as.double(obj$ry[1]),
		as.double(obj$ry[2])
		)
         alph <- obj$alph
         if (length(alph) <= 1) {
                 mm <- 1.5 * sqrt((obj$rx[2] - obj$rx[1])^2 + 
                         (obj$ry[2] - obj$ry[1])^2)
                 alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
         }
	.C("VR_alset",
		as.double(obj$alph),
		as.integer(length(obj$alph))
		)
	dx <- (xu - xl)/n
	dy <- (yu - yl)/n
	xs <- seq(xl, xu, dx)
	ys <- seq(yl, yu, dy)
	z <- matrix(nrow = length(xs), ncol = length(ys))
	for(i in 1:length(ys))
		z[, i] <- trval(obj, xs, rep(ys[i], length(xs))) + 
		predval(obj, xs, rep(ys[i], length(xs)))
	invisible(list(x = xs, y = ys, z = z))
}
"predval"<-
function(obj, xp, yp)
{
	npt <- length(xp)
	z <- .C("VR_krpred",
		z = double(npt),
		as.double(xp),
		as.double(yp),
		as.double(obj$x),
		as.double(obj$y),
		as.integer(npt),
		as.integer(length(obj$x)),
		as.double(obj$yy)
	)
	z$z
}
"semat"<-
function(obj, xl, xu, yl, yu, n, se)
{
	if(!inherits(obj, "trgls")) stop("object not from kriging")
	if(n > 999) stop("n is too large")
	.C("VR_frset",
		as.double(obj$rx[1]),
		as.double(obj$rx[2]),
		as.double(obj$ry[1]),
		as.double(obj$ry[2])
		)
         alph <- obj$alph
         if (length(alph) <= 1) {
                 mm <- 1.5 * sqrt((obj$rx[2] - obj$rx[1])^2 + 
                         (obj$ry[2] - obj$ry[1])^2)
                 alph <- c(alph[1], obj$covmod(seq(0, mm, alph[1])))
         }
	.C("VR_alset",
		as.double(obj$alph),
		as.integer(length(obj$alph))
		)
	dx <- (xu - xl)/n
	dy <- (yu - yl)/n
	xs <- seq(xl, xu, dx)
	ys <- seq(yl, yu, dy)
	z <- matrix(nrow = length(xs), ncol = length(ys))
	np <- obj$np
	npar <- ((np + 1) * (np + 2))/2
	if(missing(se))
		se <- sqrt(sum(obj$W^2)/(length(obj$x) - npar))
	for(i in 1:length(ys))
		z[, i] <- se * sqrt(seval(obj, xs, rep(ys[i], length(xs))))
	invisible(list(x = xs, y = ys, z = z))
}
"seval"<-
function(obj, xp, yp)
{
	npt <- length(xp)
	np <- obj$np
	npar <- ((np + 1) * (np + 2))/2
	z <- .C("VR_prvar",
		z = double(npt),
		as.double(xp),
		as.double(yp),
		as.integer(npt),
		as.double(obj$x),
		as.double(obj$y),
		as.double(obj$l),
		as.double(obj$r),
		as.integer(length(obj$x)),
		as.integer(np),
		as.integer(npar),
		as.double(obj$l1f))
	z$z
}

"correlogram"<-
function(krig, nint, plotit=T, ...)
{
	if(!is.loaded(symbol.C("VR_correlogram"))) stop("Compiled code has not been dynamically loaded")
	z <- .C("VR_correlogram",
		xp = double(nint),
		yp = double(nint),
		nint = as.integer(nint),
		as.double(krig$x),
		as.double(krig$y),
		if(krig$np > 0) as.double(krig$wz) else as.double(krig$z),
		as.integer(length(krig$x)),
		cnt = integer(nint)
	        )
	xp <- z$xp[1:z$nint]
	yp <- z$yp[1:z$nint]
	z <- list(x = xp, y = yp, cnt = z$cnt[1:z$nint])
	if(plotit)
	   if(exists(".Device")) {
	      plot(xp, yp, type = "p", ylim = c(-1, 1), ...)
	      abline(0, 0)
	      invisible(z)
	   }
	   else {
	      warning("Device not active")
	      return(z)
	   }
	z
}

"variogram"<-
function(krig, nint, plotit=T, ...)
{
	if(!is.loaded(symbol.C("VR_variogram"))) stop("Compiled code has not been dynamically loaded")
	z <- .C("VR_variogram",
		xp = double(nint),
		yp = double(nint),
		nint = as.integer(nint),
		as.double(krig$x),
		as.double(krig$y),
		if(krig$np > 0) as.double(krig$wz) else as.double(krig$z),
		as.integer(length(krig$x)),
		cnt = integer(nint)
		)
	xp <- z$xp[1:z$nint]
	yp <- z$yp[1:z$nint]
	if(xp[1] > 0) {xp <- c(0, xp); yp <- c(0, yp)}
	z <- list(x = xp, y = yp, cnt = z$cnt[1:z$nint])
	if(plotit)
	   if(exists(".Device")) {
	      plot(xp, yp, type = "p", ...)
	      invisible(z)
	   }
	   else {
	      warning("Device not active")
	      return(z)
	   }
	z
}
expcov <- function(r, d, alpha=0, se=1)
{
	se^2*(alpha*(r < d/10000) + (1-alpha)*exp(-r/d))
}

gaucov <- function(r, d, alpha=0, se=1)
{
	se^2*(alpha*(r < d/10000) + (1-alpha)*exp(-(r/d)^2))
}

sphercov <- function(r, d, alpha=0, se=1, D=2)
{
	r <- r/d
	if(D==2) {
	t <- 1 - (2/pi)*(r*sqrt(1-r^2) + asin(r))
	} else {
	t <- 1 - 1.5*r + r^3/2
	}
	se^2*(alpha*(r < 1/10000) + (1-alpha)*t*(r < 1))
}
"ppinit"<- function(file)
{
  file.exists <- function(name) {
    if(machine() == "Unix") system(paste("test -r", name), intern = F) == 0
  }
  tfile <- file
  t1file <- paste(.Library,"spatial/data",file,sep="/")
  if(file.exists(t1file)) tfile <- t1file
  h <- scan(tfile, list(xl = 0, xu = 0, yl = 0, yu = 0, fac = 0),
	    nmax = 5, skip = 2)
  pp <- scan(tfile, list(x = 0, y = 0), skip = 3)
  pp$x <- pp$x/h$fac
  pp$y <- pp$y/h$fac
  pp$xl <- h$xl/h$fac
  pp$xu <- h$xu/h$fac
  pp$yl <- h$yl/h$fac
  pp$yu <- h$yu/h$fac
  if(!is.loaded(symbol.C("VR_ppset"))) 
    stop("Compiled code has not been dynamically loaded")
  ppregion(pp)
  invisible(pp)
}

"Kfn"<- function(pp, fs, k = 100)
{
  z <- .C("VR_sp_pp2",
	  as.double(pp$x),
	  as.double(pp$y),
	  as.integer(length(pp$x)),
	  k1 = as.integer(k),
	  h = double(k),
	  dmin = double(1),
	  lm = double(1),
	  as.double(fs))
  list(y = z$h[1:z$k1], x = (seq(1:z$k1) * fs)/k, k = k, 
       dmin = z$dmin, lm = max(z$dmin, z$lm),
       call=match.call())
}

"Kenvl"<- function(fs, nsim, ...)
{
  dot.expression <- as.expression(substitute(...))
  h <- Kfn(pp = eval(dot.expression), fs)
  hx <- h$x
  hu <- h$y
  hl <- h$y
  ha <- h$y^2
  for(i in 2:nsim) {
    h <- Kfn(pp = eval(dot.expression), fs)$y
    hu <- pmax(hu, h)
    hl <- pmin(hl, h)
    ha <- ha + h^2
  }
  list(x = hx, lower = hl, upper = hu, aver = sqrt(ha/nsim),
       call=match.call())
}

"Kaver"<- function(fs, nsim, ...)
{
  dot.expression <- as.expression(substitute(...))
  h <- Kfn(pp = eval(dot.expression), fs)
  hx <- h$x
  ha <- h$y^2
  for(i in 2:nsim) {
    h <- Kfn(pp = eval(dot.expression), fs)$y
    ha <- ha + h^2
  }
  list(x = hx, y = sqrt(ha/nsim), call=match.call())
}

"ppregion"<- function(xl = 0, xu = 1, yl = 0, yu = 1)
{
  if(is.list(xl))  .C("VR_ppset", as.double(xl$xl), as.double(xl$xu), 
		      as.double(xl$yl), as.double(xl$yu))
  else  .C("VR_ppset", as.double(xl), as.double(xu),
	   as.double(yl), as.double(yu))
  invisible()
}

"Psim"<- function(n)
{
  z <- .C("VR_pdata",
	  as.integer(n),
	  x = double(n),
	  y = double(n))
  invisible(list(x = z$x, y = z$y, call=match.call()))
}

"Strauss"<- function(n, c = 0, r)
{
  init <-  0
  if(!exists("ppx")) {
    init <-  1
    z <- .C("VR_pdata",
	    as.integer(n),
	    x = double(n),
	    y = double(n))
    assign("ppx", z$x)
    assign("ppy", z$y)
  }
  z <- .C("VR_simpat",
	  as.integer(n),
	  x = ppx,
	  y = ppy,
	  as.double(c),
	  as.double(r),
	  as.integer(init))
  assign("ppx", z$x)
  assign("ppy", z$y)
  invisible(list(x = z$x, y = z$y, call=match.call()))
}

"SSI"<- function(n, r)
{
  z <- .C("VR_simmat",
	  as.integer(n),
	  x = double(n),
	  y = double(n),
	  as.double(r))
  invisible(list(x = z$x, y = z$y, call=match.call()))
}

"pplikfn"<- function(cc, R, n, x, y, ng, target, trace=F)
{
  z <- .C("VR_plike",
	  as.double(x),
	  as.double(y),
	  as.integer(n),
	  as.double(cc),
	  as.double(R),
	  as.integer(ng),
	  as.double(target),
	  res=double(1)
	  )
  if(trace) print(c(cc, z$res))
  z$res
}

"pplik"<- function(pp, R, ng=50, trace=F)
{
  n <- length(pp$x)
  target <- n * (Kfn(pp, R,1)$y)^2 * pi /
    ((pp$xu - pp$xl) * (pp$yu - pp$yl))
  if(target == 0) return(0)
  tmp <- pplikfn(1, R, n, pp$x, pp$y, ng, target, F)
  if(tmp <= 0) return(1)
  uniroot(pplikfn, c(0,1), lower=-target, upper=tmp,
	  R=R, n=n, x=pp$x, y=pp$y, ng=ng, target=target,
	  trace=trace)$root
}

library.dynam("spatial")
