bk.grid <- function(point.obj,
                    at,
                    prior,
                    var.mod.obj,
                    xsw=NULL,ysw=NULL,xne=NULL,yne=NULL,
                    dx=NULL,dy=NULL,nx=NULL,ny=NULL,
                    angle=0,
                    maxdist = NULL,
                    extrap = F,
                    border=NULL,
                    trend=0,
                    rsearch=0,
                    nsearch=0,
                    nsmin=-1,
                    nsmax=-1,
                    mode=3,
                    duplicate = "error",
                    dupfun = NULL,
                    method="gqr",
 		    get.lm=FALSE)
{
    tmp   <- check.gridparams(angle,xsw,xne,ysw,yne,
                              dx,dy,nx,ny)
    angle <- tmp$angle
    xsw   <- tmp$xsw
    xne   <- tmp$xne
    ysw   <- tmp$ysw
    ysw   <- tmp$ysw
    dx    <- tmp$dx 
    dy    <- tmp$dy 
    nx    <- tmp$nx 
    ny    <- tmp$ny 
    nz    <- tmp$nz 

    dog   <- check.border.grid(extrap,xsw,xne,ysw,yne,nx,ny,border,point.obj)

    
    point.obj <- remove.duplicates(point.obj,at,duplicate,dupfun)
    tmp       <- check.krigedata(point.obj,at,var.mod.obj,mode)
    n         <- tmp$n

    tmp     <- check.searchparams(maxdist,rsearch,nsearch,nsmin,nsmax,n)
    rsearch <- tmp$rsearch
    nsmin   <- tmp$nsmin
    nsmax   <- tmp$nsmax
    
    if(trend==0) ntrend<-1
    if(trend==1) ntrend<-3
    if(trend==2) ntrend<-6

    covtype<-switch(attr(var.mod.obj,"type"),
                    exponential=1,
                    gaussian=2,
                    spherical=3,
                    linear=4,
                    0)
    cov0<-0
    
#   determine optimum array sizes:
    if(!is.null(method)){
        if(method!="gqr" && method!="direct" && method!="ols")
            stop("method (used for glsfit) should be one of \"gqr\", \"ols\" or \"direct\"!")
    } else {
        method <-"gqr"
    }
    method<-switch(method,direct=2,gqr=1,ols=0)    
#    lwork <- glsfit.workquery(n,ntrend,method)
    lwork <- 3000

    if(prior$ntr!=ntrend)
        stop("model order of priors does not match!")
    npr<-prior$n
    typpr<-prior$info
    typpr[prior$type=="subjective"]<-typpr[prior$type=="subjective"]*(-1)


        
    snbbit<-rep(0,1+n*nz)
    snbbit[1]<-1
    if(get.lm){
      mu <- double(ntrend*nz)
      lambda <- double(n*nz)
      lambd0 <- double(nz)
    } else {
      mu <- double(ntrend)
      lambda <- double(n)
      lambd0 <- double(1)
    }
    

    ans<-.C("bk_grid",
            xsw     = as.double(xsw),
            ysw     = as.double(ysw),
	    xne     = as.double(xne),
            yne     = as.double(yne),
            angle   = as.double(angle),
	    nx      = as.integer(nx),
	    ny      = as.integer(ny),
            dx      = as.double(dx),
            dy      = as.double(dy),
	    xg      = double(nx),
	    yg      = double(ny),
	    zg      = double(nz),
	    varg    = double(nz),
	    dog     = as.integer(dog),
            lon     = as.double(point.obj$x),
            lat     = as.double(point.obj$y),
            z       = as.double(point.obj[[match(at, names(point.obj))]]),
	    extrap  = as.integer(extrap),
	    n       = as.integer(n),
	    covtype = as.integer(covtype),
            covpar  = as.double(var.mod.obj$parameters),
	    covmat  = double(n*n), 	
	    ldcov   = as.integer(n),
            extcov  = as.integer(0), # no external cov matrix
	    trend   = as.integer(trend),
	    ntrend  = as.integer(ntrend),
            mupr    = as.double(matlist.cbind(prior$mu)),
            ldmpr   = as.integer(ntrend),
            phipr   = as.double(matlist.cbind(prior$phi)),
            ldphpr  = as.integer(ntrend),
            lonpr   = as.double(prior$lon),
            latpr   = as.double(prior$lat),
            npr     = as.integer(npr),
            typpr   = as.integer(typpr),
            rsearch = as.double(rsearch),
	    nsearch = as.integer(nsearch),
	    nsmin   = as.integer(nsmin),
	    nsmax   = as.integer(nsmax),
            lwork   = as.integer(lwork),
            mode    = as.integer(mode),
            mu      = as.double(mu),
            lambda  = as.double(lambda),
            lambd0  = double(lambd0),
            bits    = as.integer(c(integer(nz),snbbit)),
            ierr    = integer(1),
            get.lm  = as.integer(get.lm),
            glsmth  = as.integer(method),
#           ,.Package= "baykrig"
            )
#browser()
    retval<-list(x=ans$xg,
                 y=ans$yg,
                 z=matrix(ans$zg,nx,ny),
                 var=matrix(ans$varg,nx,ny),
                 done=matrix(ans$dog, nx, ny),
                 snb=matrix(ans$bits[(nz+2):(nz+n*nz+1)],nrow=n,ncol=nz,byrow=F)
                 )
    if(get.lm){
      retval$lambda <- matrix(ans$lambda,nrow=n,ncol=nz,byrow=FALSE)
      retval$lambda0 <- matrix(ans$lambd0,nx,ny)
      retval$mu <- matrix(ans$mu,nrow=n,ncol=ntrend,byrow=FALSE)
    }
    retval$z[retval$done<=0] <- NA
    retval$var[retval$done<=0] <- NA
    retval$lambda0[retval$done<=0] <- NA
    retval$data<-point.obj
    retval$at<-at
    class(retval)<-"krige.map"
    retval    
  }