In this paper we will consider the problem of the numerical simulation of n
on-Gaussian, scalar random fields with a prescribed correlation structure e
ither provided by a theoretical model or computed on a set of observational
data. Although the numerical generation of a generic, non-Gaussian random
field is a trivial operation, the task becomes difficult when the field is
constrained with a prescribed correlation structure. In this regard, three
numerical methods, useful for astronomical applications, are presented. The
limits and capabilities of each method are discussed, and the pseudocodes
describing the numerical implementation are provided for two of them.