In this paper we will consider the problem of the numerical simulation of non-Gaussian, scalar random fields with a prescribed correlation structure provided either 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 tough when constraining the field with a prefixed correlation structure. At this regards, three numerical methods, useful for astronomical applications, are presented. The limits and capabilities of each method are discussed and the pseudo-codes describing the numerical implementation are provided for two of them.