Circulant Embedding method for generating stationary Gaussian field

simulating stationary Gaussian field over an 'm' times 'n' grid
INPUT:
- 'm' and 'n' for evaluating the field over the m*n grid;
note that size of covariance matrix is m^2*n^2;
- scalar function rho(h), where 'h' is a two dimensional vector
input and cov(X_t,Y_s)=rho(t-s) is the cov. function of a
2-dimensional stationary Gaussian field; see reference below;
OUTPUT:
- two statistically independent fields 'field1' and 'field2'
over the m*n grid;
- vectors 'tx' and 'ty' so that the field is plotted via
imagesc(tx,ty,field1)
Example:
rho=@(h)((1-h(1)^2/50^2-h(1)*h(2)/(15*50)-h(2)^2/15^2)...
*exp(-(h(1)^2/50^2+h(2)^2/15^2))); % define covariance function
stationary_Gaussian_process(512,384,rho); % plot when no output wanted
Reference:
Kroese, D. P., & Botev, Z. I. (2015). Spatial Process Simulation.
In Stochastic Geometry, Spatial Statistics and Random Fields(pp. 369-404)
Springer International Publishing, DOI: 10.1007/978-3-319-10064-7_12