Fits a trend surface by generalized least-squares.

Syntax:
       surf.gls(np, covmod, x, y, z, nx=1000, ...)

Input:
  np      degree of polynomial surface
  covmod  function to evaluate covariance or correlation function
  x       coordinates or a data frame with columns `x', `y', `z'
  y       coordinates
  z       coordinates. Will supersede `x$z'
  nx      Number of bins for table of the covariance.
          Increasing adds accuracy, and increases size of the object.
...       parameters for `covmod'

Output:
list with components the coefficients  and others for internal use
only.

See also:
`trmat',  `surf.ls',
`prmat',  `semat', `expcov', `gaucov', `sphercov'

Example:
topo.kr <-surf.gls(2, expcov, topo, d=0.7)
trsurf <- trmat(topo.kr, 0,6.5, 0, 6.5, 50)
contour(trsurf)
prsurf <- prmat(topo.kr, 0,6.5, 0, 6.5, 50)
contour(prsurf, levels=seq(700,  925,  25))
sesurf <- semat(topo.kr, 0, 6.5, 0, 6.5, 30)
contour(sesurf,levels=c(22,25))