# $Id: fit.variogram.reml.q,v 1.12 2006-12-10 14:17:17 edzer Exp $

"fit.variogram.reml" <-
function (formula, locations, data, model, debug.level = 1, set, degree = 0)
{
    if (missing(formula)) 
        stop("nothing to fit to")
    if (!inherits(formula, "formula"))
		stop("first argument should be of class formula")
	if (!missing(locations)) {
    	if (inherits(locations, "formula"))
			coordinates(data) = locations
		else if (is(locations, "Spatial"))
			data = as(locations, "SpatialPointsDataFrame")
	}
    if (!is(data, "SpatialPointsDataFrame"))
        stop("data should (now) be of class SpatialPointsDataFrame")
    if (missing(model)) 
        stop("no model to fit")
    if (!inherits(model, "variogramModel"))
        stop("model should be of class variogramModel (use vgm)")
    fit.sills = rep(TRUE, length(model$model))
    fit.ranges = rep(FALSE, length(model$model))
    .Call("gstat_init", as.integer(debug.level))
    ret = gstat.formula(formula, data)
    ret$y <- residuals(lm(formula, data))
    .Call("gstat_new_data", as.double(ret$y), as.double(ret$locations),
		as.double(ret$X), as.integer(1), double(0), as.integer(-1), 
		as.integer(0), as.double(-1), as.integer(1), 
		double(0), double(0), as.integer(degree),
		as.integer(is.projected(data)))
    load.variogram.model(model)
    if (!missing(set))
    	gstat.load.set(set)
    ret = .Call("gstat_fit_variogram", as.integer(5), 
        as.integer(fit.sills), as.integer(fit.ranges))
    .Call("gstat_exit", 0) 
    model$psill = ret[[1]]
    model$range = ret[[2]]
    model
}