"matern"<-function(x=seq(0,4*range,len=100),scale=1,range=1,smoothness=0.5)
{
#
#  usage:   out<-matern(x,scale,range,smoothness)
#
#  This function will compute the Matern covariance between
#       two points. The vector x represents the distances between pairs
#       of points. By default the correlation is calculated for uniformly
#       spaced distances from 0 to 4 times the range.
#
#
#  The other input parameters are:
#
#  scale:      the point variance of the covariance function (default 1).
#
#  range:      the spatial scale of correlation. The correlation is small 
#              for distances above about twice this value.
#
#  smoothness: the smoothness of the corresponding random field. The 
#              random field will have ceiling(smoothness) - 1 mean-square
#              derivatives. For Gaussian random fields this means that
#              realizations from this field will be ceiling(smoothness) - 1
#              times differentiable.
#              If the smoothness is above 50 the Squared Exponential
#              limiting covariance is used (default 0.5 i.e. Exponential).
#
#
#  The output contains, in addition to pertinent input:
#  
#  The output is a list whose elements are the pertinent input
#       listed above and 
#
#  y:        vector of Matern covariance at distances specified by the
#            vector x and parameters (scale, range and smoothness)
#
#            The names() function will list the names of the output.
#
#  For technical background see
#
#                        Handcock and Stein (1993),
#                     'A Bayesian Analysis of Kriging',
#                       Technometrics, 35, 4, 403-410.
#
#  or contact              mhandcock@stern.nyu.edu
#
# Has the Fortran subroutine rkmat has been loaded yet?
# If not, attempt to load it. You may need to alter the load fro your system.
#
	if(!is.loaded(symbol.For("rkmat"))) {
		temp <- dyn.load("rkmat.o", 2)	
#	Plan b
#		temp <- dyn.load(paste(getenv("SLOCAL"), ".Functions/.Dynload", 
#			"rkmat.o", sep = "/"), 2)	
		cat("rkmat.o  FORTRAN subroutine dynamically loaded:", temp, 
			fill = T)
	}
	y <- x
	length <- length(x)
	theta <- c(scale,range,smoothness)
	storage.mode(x) <- "double"
	storage.mode(y) <- "double"
	storage.mode(theta) <- "double"
	storage.mode(length) <- "integer"
	out <- .Fortran("rkmat",
	     theta=theta,
	     x=x,
	     y=y,
	     n=length
	     )
	list(x=out$x,y=out$y,scale=scale,range=range,smoothness=smoothness)
}