# $Id: variogram.default.q,v 1.26 2008-04-08 12:20:43 edzer Exp $

"variogram.default" <-
function (object, locations, X, cutoff, width = cutoff/15.0, alpha = 0, 
    beta = 0, tol.hor = 90/length(alpha), tol.ver = 90/length(beta), 
    cressie = FALSE, dX = numeric(0), boundaries = numeric(0), 
    cloud = FALSE, trend.beta = NULL, debug.level = 1, cross = TRUE, 
	grid, map = FALSE, g = NULL, ..., projected = TRUE) 
{
    id1 = id2 = 0
    ret = NULL
    if (missing(cutoff)) {
		if (is.logical(map) && map == TRUE)
			stop("for variogram maps, supply at least a cutoff")
        cutoff = numeric(0)
	}
    else if (width <= 0) 
        stop("argument width should be positive")
    if (missing(width)) 
        width = numeric(0)
    if (cloud == TRUE) 
        width = 0
	if (is.logical(map) && map == TRUE) {
		require(sp)
		cells.dim = length(seq(-cutoff, cutoff, by = width))
		grid.topology = GridTopology(rep(-cutoff, 2), rep(width, 2), rep(cells.dim, 2))
		map = SpatialGrid(grid = grid.topology)
	}
    .Call("gstat_init", as.integer(debug.level))
	if (!is.null(g) && !is.null(g$set))
		gstat.load.set(g$set)
    id.names = NULL
    if (is.list(object) && is.list(locations)) {
        nvars = length(object)
        if (!is.null(names(object))) 
            id.names = names(object)
        for (i in 1:nvars) {
            if (missing(X)) 
                Xloc = rep(1, length(object[[i]]))
            else Xloc = X[[i]]
            t.beta = numeric(0)
            if (!is.null(trend.beta) && length(trend.beta) > 0) 
                t.beta = trend.beta[[i]]
            else t.beta = numeric(0)
			if (missing(grid) || !is.list(grid)) {
				if (!is.null(g) && gridded(g$data[i]$data))
					grd = unlist(gridparameters(g$data[i]$data))
				grd = numeric(0)
			} else
				grd = grid[[i]]
            .Call("gstat_new_data", as.double(object[[i]]), 
				as.double(locations[[i]]), as.double(Xloc), 
				as.integer(1), as.double(t.beta), as.integer(-1),
				as.integer(0), as.double(-1), as.integer(1), 
				double(0), grd, as.integer(0), as.integer(projected))
			if (!is.null(g) && !is.null(g$model[[id.names[i]]])) 
				load.variogram.model(g$model[[id.names[i]]], c(i - 1, i - 1))
        }
    } else 
		stop("argument object and locations should be lists")

	if (inherits(map, "SpatialGrid"))
		map = as.double(unlist(gridparameters(as(map, "SpatialGrid"))))

    pos = 0
    ids = NULL
	is.direct = NULL
	if (cross)
		id.range = nvars:1
	else
		id.range = 1:nvars
    for (id1 in id.range) {
        for (id2 in ifelse(cross, 1, id1):id1) {
            if (is.null(id.names)) 
                id = ifelse(id1 == id2, paste(id1), cross.name(id2, id1))
            else id = ifelse(id1 == id2, paste(id.names[id1]), 
                cross.name(id.names[id2], id.names[id1]))
            for (a in alpha) {
                for (b in beta) {
                  direction = as.numeric(c(a, b, tol.hor, tol.ver))
                  ret.call = .Call("gstat_variogram", 
				  		as.integer(c(id1 - 1, id2 - 1)), 
						as.numeric(cutoff), as.numeric(width), 
                    	as.numeric(direction), as.integer(cressie), 
                    	as.numeric(dX), as.numeric(boundaries), map)
				  boundaries = numeric(0)
                  if (is.logical(map) && map == FALSE) {
                    np = ret.call[[1]]
                    sel = np > 0
                    n.dir = length(sel[sel])
                    if (n.dir > 0) {
                      dist = ret.call[[2]]
                      gamma = ret.call[[3]]
                      dir.a = rep(a, n.dir)
                      dir.b = rep(b, n.dir)
                      ids = c(ids, rep(id, n.dir))
					  is.direct = c(is.direct, id1 == id2)
                      df = data.frame(np = np[sel], dist = dist[sel], 
                        gamma = gamma[sel], dir.hor = dir.a, dir.ver = dir.b)
                      if (pos > 0) 
                        ret[(pos + 1):(pos + n.dir), ] = df
                      else ret = df
                      pos = pos + n.dir
                    }
				  } else {
				  	if (is.null(ret)) {
					  ret = data.frame(ret.call[[1]], ret.call[[2]],
							ret.call[[4]], ret.call[[3]])
					  names = c("dx", "dy", id, paste("np", id, sep="."))
					} else { 
					  ret = data.frame(ret, ret.call[[4]], ret.call[[3]])
					  names = c(names, id, paste("np", id, sep="."))
					}
					names(ret) = names
				  }
                }
            }
        }
    }
    .Call("gstat_exit", NULL)
	if (is.logical(map) && map == FALSE) {
    	ret$id = factor(ids, levels = unique(ids))
		attr(ret, "direct") = data.frame(id = unique(ids), is.direct = is.direct)
    	if (cloud) {
        	class(ret) = c("variogramCloud", "data.frame")
			attr(ret, ".BigInt") = 2^(4 * .Machine$sizeof.long)
    	} else 
			class(ret) = c("gstatVariogram", "data.frame")
	} else {
		coordinates(ret) = c("dx", "dy")
		gridded(ret) = TRUE
		ret = list(map = ret)
		class(ret) = c("variogramMap", "list")
	}
    ret
}