Covariance functions return the value of the covariance
  \(C(h)\) between a pair variables located at points separated by the
  distance \(h\).
  The covariance function can be written as a product of a variance
  parameter \(\sigma^2\) times a positive definite
  correlation function \(\rho(h)\):
  $$C(h) = \sigma^2 \rho(h).$$
  The expressions of the covariance functions available in geoR
  are given below. We recommend the LaTeX (and/or the corresponding
  .dvi, .pdf or .ps) version of this document for
  better visualization of the formulas.
  
Denote \(\phi\) the basic parameter of the correlation
  function and name it the range parameter.
  Some of the correlation functions will have an extra parameter
  \(\kappa\), the smoothness parameter.
  \(K_\kappa(x)\) denotes the modified Bessel
  function of the third kind of order \(\kappa\). See
  documentation of the function besselK for further details.
  In the equations below the functions are valid for  \(\phi>0\) and \(\kappa>0\), unless stated otherwise.
  
cauchy 
  $$\rho(h) = [1+(\frac{h}{\phi})^2]^{-\kappa}$$
gencauchy (generalised Cauchy) 
  $$\rho(h) = [1+(\frac{h}{\phi})^{\kappa_{2}}]^{-{\kappa_1}/{\kappa_2}},
    \kappa_1 > 0, 0 < \kappa_2 \leq 2 $$
    
circular
  Let \(\theta = \min(\frac{h}{\phi},1)\) and
  $$g(h)= 2\frac{(\theta\sqrt{1-\theta^2}+
      \sin^{-1}\sqrt{\theta})}{\pi}.$$
  Then, the circular model is given by:
  $$\rho(h) = \left\{ \begin{array}{ll}
    1 - g(h)  \mbox{ , if $h < \phi$}\cr
    0    \mbox{ , otherwise}
    \end{array} \right.$$
  
cubic
  $$\rho(h) = \left\{ \begin{array}{ll}
    1 - [7(\frac{h}{\phi})^2 - 8.75(\frac{h}{\phi})^3 +
    3.5(\frac{h}{\phi})^5-0.75(\frac{h}{\phi})^7] \mbox{ , if $h<\phi$} \cr
    0  \mbox{ , otherwise.}
    \end{array} \right.$$
  
gaussian
  $$\rho(h) = \exp[-(\frac{h}{\phi})^2]$$
exponential
  $$\rho(h) = \exp(-\frac{h}{\phi})$$
matern
  $$\rho(h) =
    \frac{1}{2^{\kappa-1}\Gamma(\kappa)}(\frac{h}{\phi})^\kappa
    K_{\kappa}(\frac{h}{\phi})$$
spherical
  $$\rho(h) = \left\{ \begin{array}{ll}
    1 - 1.5\frac{h}{\phi} + 0.5(\frac{h}{\phi})^3
    \mbox{ , if $h$ < $\phi$} \cr
    0    \mbox{ ,  otherwise}
    \end{array} \right.$$
    
power (and linear)
  The parameters of the this model
  \(\sigma^2\) and \(\phi\) can not be  
  interpreted as partial sill and range
  as for the other models.
  This model implies an unlimited dispersion and,
  therefore, has no sill and corresponds to a process which is only
  intrinsically stationary.
  The variogram function is given by:
  $$\gamma(h) = \sigma^2 {h}^{\phi} \mbox{ , } 0 < \phi < 2,
    \sigma^2 > 0$$
Since the corresponding process is not second order stationary the
  covariance and correlation functions are not defined.
  For internal calculations the geoR
  functions uses the fact the this model possesses locally
  stationary representations with covariance functions of the form:
  $$C_(h) = \sigma^2 (A - h^\phi)$$ ,
  where \(A\) is a suitable constant as given in
  Chiles & Delfiner (pag. 511, eq. 7.35).
  
The linear model corresponds a particular case with
  \(\phi = 1\).
  
powered.exponential (or stable)
  $$\rho(h) = \exp[-(\frac{h}{\phi})^\kappa] \mbox{ , } 0 < \kappa
    \leq 2$$
  
gneiting
  $$C(h)=\left(1 + 8 sh + 25 (sh)^2 + 32
    (sh)^3\right)(1-sh)^8 1_{[0,1]}(sh)$$
  where
  \(s=0.301187465825\).
  For further details see documentation of the function
  CovarianceFct 
  in the package
  RandomFields from where we extract the following :
  It is an alternative to the gaussian model since
  its graph is visually hardly distinguishable from the graph of
  the Gaussian model, but possesses neither the mathematical and nor the
  numerical disadvantages of the Gaussian model.
gneiting.matern
  Let \(\alpha=\phi\kappa_2\), \(\rho_m(\cdot)\) denotes the \(\mbox{Mat\'{e}rn}\) model
  and \(\rho_g(\cdot)\) the Gneiting model. Then the
  \(\mbox{Gneiting-Mat\'{e}rn}\) is given by
  $$\rho(h) = \rho_g(h|\phi=\alpha) \,
    \rho_m(h|\phi=\phi,\kappa=\kappa_1)$$
wave
  $$\rho(h) = \frac{\phi}{h}\sin(\frac{h}{\phi})$$ 
pure.nugget 
  $$\rho(h) = k$$
  where k is a constant value. This model corresponds to
  no spatial correlation.
Nested models
  Models with several structures
  usually called nested models
  in the geostatistical literature are also allowed.
  In this case the argument cov.pars takes a matrix and
  cov.model and lambda can either have length equal to
  the number of rows of this matrix or length 1.
  For the latter  cov.model and/or lambda are recycled, i.e. the same
  value is used for all structures.