function OUT = cov_stat(y,S,k)

%COV_STAT simulates covariance-stationary processes. 
%
%INPUTS: (i)   y = spherical parameters [range,sill,nug]
%        (ii)  S = (Nx2)-matrix of locations (Xi,Yi).
%        (iii) k = number of simulations desired.
%
%OUTPUT: OUT = (Nxk-matrix of with each column,(Zi), giving 
%              a vector of simulated values at all locations.

%Compute Covariogram

n = length(S) ;

D = distance_mat(S);

U = ones(n,1) ;

C = y(2)*U*U' - spher_mat(y,D) ;

%Correct C on diagonal (where nugget effect appears)

C = C + (y(3)*eye(n)) ; 

%Compute Cholesky Decompostion

L = chol(C)' ;

%simulate data

X = normrnd(0,1,n,k) ;

%convert to covariance stationary data

OUT = L * X ;