Algorithm 1013

A continuous spectral algorithm and computer routines in the R programming environment that enable the simulation of second-order stationary and intrinsic (i.e., with second-order stationary increments or generalized increments) vector Gaussian random fields in Euclidean spaces are presented. The simulation is obtained by computing a weighted sum of cosine and sine waves, with weights that depend on the matrix-valued spectral density associated with the spatial correlation structure of the random field to simulate. The computational cost is proportional to the number of locations targeted for simulation, below that of sequential, matrix decomposition and discrete spectral algorithms. Also, the implementation is versatile, as there is no restriction on the number of vector components, workspace dimension, number and geometrical configuration of the target locations. The computer routines are illustrated with synthetic examples and statistical testing is proposed to check the normality of the distribution of the simulated random field or of its generalized increments. A by-product of this work is a spectral representation of spherical, cubic, penta, Askey, J-Bessel, Cauchy, Laguerre, hypergeometric, iterated exponential, gamma, and stable covariance models in the d-dimensional Euclidean space.