randomfieldgenerator.C

#include "random.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <limits>

/*
    GenRan: Non-Gaussian Random field generator.
    Copyright (C) 2006 Peter Grassl
    Distributed under the MIT licence.

    Top-level entry points:
      generateRandomField    — orchestrates target-spectrum computation,
                               optional SDB-2011 iteration, Gaussian field
                               generation, and translation to the prescribed
                               marginal CDF.
      mapGaussianToRandom    — applies the memoryless translation
                               X = F_NG^{-1}(Φ(g/σ_G)) to a Gaussian field.

    The optional iteration mode implements
        Shields, Deodatis & Bocchini, "A simple and efficient methodology to
        approximate a general non-Gaussian stationary stochastic process by
        a translation process," Probab. Eng. Mech. 26 (2011) 511–519.
    It finds an underlying Gaussian PSDF whose translation image best
    matches the prescribed target spectrum, in 1D, 2D, and 3D.
*/

#define MYPI 3.14159265


// ===========================================================================
//  Helpers shared by the Gaussian-to-non-Gaussian mapping and the SDB-2011
//  iterative scheme.
// ===========================================================================

// Returns F_NG^{-1}(p): the non-Gaussian quantile corresponding to a
// probability p in [0, 1]. typeOfCDF and randomParameter follow the same
// conventions as the rest of the program (1 = Gaussian, 2 = Weibull,
// 3 = grafted Weibull-Gaussian).
static double inverseNonGaussianCdf(double p, double *params, int typeOfCDF)
{
  if (typeOfCDF == 1) {
    const double std = params[0] * params[1];
    return normal_cdf_inv(p, params[0], std);
  }
  if (typeOfCDF == 2) {
    const double scaling = params[0];
    const double modulus = params[1];
    return scaling * std::pow(-std::log(1.0 - p), 1.0 / modulus);
  }
  if (typeOfCDF == 3) {
    const double scaling = params[0];
    const double modulus = params[1];
    const double graftingProbability = params[2];
    const double rescaling = params[4];
    const double gaussianMean = params[5];
    const double gaussianStd  = params[6];
    const double pp = p / rescaling;
    if (pp < graftingProbability) {
      return scaling * std::pow(-std::log(1.0 - pp), 1.0 / modulus);
    }
    const double length = std::pow(-std::log(1.0 - graftingProbability / rescaling),
                                   1.0 / modulus);
    const double helpVal = normal_01_cdf((length - gaussianMean) / gaussianStd);
    double helpProbability = pp - graftingProbability / rescaling + helpVal;
    if (helpProbability >= 1.0) helpProbability = 0.99999;
    return normal_cdf_inv(helpProbability, gaussianMean, gaussianStd);
  }
  std::fprintf(stderr, "inverseNonGaussianCdf: unknown typeOfCDF %d\n", typeOfCDF);
  std::exit(1);
}


// Compute the unconditional mean and variance of the non-Gaussian
// distribution X = F_NG^{-1}(Φ(z)) by trapezoidal integration of
// h(z) and h(z)² against the standard normal density φ(z). Used by the
// SDB-2011 path for distributions whose moments do not have a clean
// analytical form (currently the grafted Weibull-Gaussian, typeOfCDF=3).
static void computeNonGaussianMomentsByQuadrature(double *params, int typeOfCDF,
                                                  double *meanOut, double *varOut)
{
  const int M = 2001;
  const double zMax = 8.0;
  const double dz = (2.0 * zMax) / (M - 1);
  const double invSqrt2Pi = 1.0 / std::sqrt(2.0 * MYPI);

  double sumX = 0.0, sumXX = 0.0;
  for (int j = 0; j < M; ++j) {
    const double z = -zMax + j * dz;
    const double w = (j == 0 || j == M - 1) ? 0.5 : 1.0;   // trapezoidal endpoints
    const double phi = invSqrt2Pi * std::exp(-0.5 * z * z);
    const double p = normal_01_cdf(z);
    const double pClamped = std::max(1.0e-12, std::min(1.0 - 1.0e-12, p));
    const double X = inverseNonGaussianCdf(pClamped, params, typeOfCDF);
    sumX  += w * X * phi;
    sumXX += w * X * X * phi;
  }
  *meanOut = sumX * dz;
  const double E_X2 = sumXX * dz;
  *varOut = E_X2 - (*meanOut) * (*meanOut);
}


// Compute the Hermite expansion coefficients of the centred translation
// kernel h(z) = F_NG^{-1}(Φ(z)) - μ_NG, namely
//
//     c_k = ∫ h(z) He_k(z) φ(z) dz,    k = 0, 1, ..., Kmax
//
// where He_k is the probabilist's Hermite polynomial and φ is the standard
// normal PDF. By Mehler's formula the autocorrelation of the translated
// process is
//
//     R_Y(τ) = Σ_{k≥1} c_k² ξ(τ)^k / k!
//
// with ξ(τ) the normalised Gaussian correlation. c_0 should evaluate to 0
// since the kernel is centred; we keep the slot for indexing convenience.
//
// Integration uses the composite trapezoidal rule on z ∈ [-zMax, zMax] with
// M points; this converges spectrally for our Gaussian-weighted integrand.
static void computeHermiteCoefficients(double *c, int Kmax,
                                       double *params, int typeOfCDF, double meanNG)
{
  const int M = 2001;        // odd, so z = 0 is a node
  const double zMax = 8.0;   // ±8σ exhausts the Gaussian weight
  const double dz = (2.0 * zMax) / (M - 1);
  const double invSqrt2Pi = 1.0 / std::sqrt(2.0 * MYPI);

  std::vector<double> hVals(M);
  std::vector<double> phi(M);
  for (int j = 0; j < M; ++j) {
    const double z = -zMax + j * dz;
    const double p = normal_01_cdf(z);
    // Clamp p away from the open boundaries so the inverse CDF stays finite
    // for the extreme tails of the Gauss-Hermite domain.
    const double pClamped = std::max(1.0e-12, std::min(1.0 - 1.0e-12, p));
    hVals[j] = inverseNonGaussianCdf(pClamped, params, typeOfCDF) - meanNG;
    phi[j] = invSqrt2Pi * std::exp(-0.5 * z * z);
  }

  // Trapezoidal integration of h(z) · He_k(z) · φ(z) over the grid.
  // Hermite polynomials are computed by the standard recurrence
  //     He_{k+1}(z) = z · He_k(z) - k · He_{k-1}(z)
  // Maintain rolling buffers to avoid storing all polynomials.
  std::vector<double> hePrev(M, 1.0);   // He_0
  std::vector<double> heCurr(M);
  for (int j = 0; j < M; ++j) heCurr[j] = -zMax + j * dz;  // He_1(z) = z

  auto trapz = [&](const std::vector<double> &he) {
    double s = 0.5 * (hVals[0]   * he[0]   * phi[0]
                    + hVals[M-1] * he[M-1] * phi[M-1]);
    for (int j = 1; j < M - 1; ++j) s += hVals[j] * he[j] * phi[j];
    return s * dz;
  };

  c[0] = trapz(hePrev);
  if (Kmax >= 1) c[1] = trapz(heCurr);

  std::vector<double> heNext(M);
  for (int k = 1; k < Kmax; ++k) {
    for (int j = 0; j < M; ++j) {
      const double z = -zMax + j * dz;
      heNext[j] = z * heCurr[j] - k * hePrev[j];
    }
    c[k + 1] = trapz(heNext);
    std::swap(hePrev, heCurr);
    std::swap(heCurr, heNext);
  }
}


// Multi-axis 1D FFT applied to a complex array of shape `N[nDim]` stored in
// row-major order with two doubles (real, imag) per element. Each axis is
// transformed in place by NR's `four1` without any normalisation, so the
// composite is an unnormalised nD FFT.
static void multidimFft(double *data, const int N[], int nDim, int isign)
{
  if (nDim == 1) {
    four1(data - 1, N[0], isign);
    return;
  }

  int total = 1;
  for (int i = 0; i < nDim; ++i) total *= N[i];

  for (int axis = 0; axis < nDim; ++axis) {
    int stride = 1;
    for (int i = axis + 1; i < nDim; ++i) stride *= N[i];
    const int outer  = total / (N[axis] * stride);  // # blocks of size N[axis]*stride
    std::vector<double> slice(2 * N[axis]);

    for (int block = 0; block < outer; ++block) {
      for (int inner = 0; inner < stride; ++inner) {
        const int baseOff = block * (N[axis] * stride) + inner;
        for (int t = 0; t < N[axis]; ++t) {
          const int off = baseOff + t * stride;
          slice[2 * t]     = data[2 * off];
          slice[2 * t + 1] = data[2 * off + 1];
        }
        four1(slice.data() - 1, N[axis], isign);
        for (int t = 0; t < N[axis]; ++t) {
          const int off = baseOff + t * stride;
          data[2 * off]     = slice[2 * t];
          data[2 * off + 1] = slice[2 * t + 1];
        }
      }
    }
  }
}


// True iff every component of the row-major linear index `lin` lies
// strictly inside (1, N[i]/2) — the "interior octant" used by the rest of
// this program (boundary-zero convention shared by computeTargetSpectralDensity).
static bool linearIndexIsInterior(int lin, const int N[], int nDim)
{
  int rem = lin;
  for (int i = nDim - 1; i >= 0; --i) {
    const int idx_i = rem % N[i];
    rem /= N[i];
    if (idx_i == 0 || idx_i >= N[i] / 2) return false;
  }
  return true;
}


// Pack a one-sided real PSDF (interior-only nonzero, row-major shape `N`)
// into a 2*total-double NR-style complex buffer with full symmetric
// extension. Imaginary parts are zero.
static void packSymmetricPsdfND(const double *S, const int N[], int nDim,
                                double *fftBuf)
{
  int total = 1;
  for (int i = 0; i < nDim; ++i) total *= N[i];
  for (int i = 0; i < 2 * total; ++i) fftBuf[i] = 0.0;

  // Iterate over the interior box [1, N[i]/2-1]^nDim and place each value at
  // all 2^nDim symmetric reflections (k_i ↔ N[i] - k_i).
  std::vector<int> idx(nDim, 1);
  while (true) {
    int linOneSided = 0;
    for (int i = 0; i < nDim; ++i) linOneSided = linOneSided * N[i] + idx[i];
    const double v = S[linOneSided];

    for (int s = 0; s < (1 << nDim); ++s) {
      int linFull = 0;
      for (int i = 0; i < nDim; ++i) {
        const int ki = (s & (1 << i)) ? (N[i] - idx[i]) : idx[i];
        linFull = linFull * N[i] + ki;
      }
      fftBuf[2 * linFull]     = v;
      fftBuf[2 * linFull + 1] = 0.0;
    }

    // Advance multi-index over the interior box.
    int axis = nDim - 1;
    while (axis >= 0) {
      idx[axis]++;
      if (idx[axis] < N[axis] / 2) break;
      idx[axis] = 1;
      axis--;
    }
    if (axis < 0) break;
  }
}


// SDB-2011 iteration generalised to nDim ∈ {1, 2, 3}. Given the user's
// prescribed target non-Gaussian PSDF S_T (interior-only nonzero in the
// same convention as computeTargetSpectralDensity) and the Hermite
// coefficients of the centred translation kernel, fills S_G_out with the
// converged underlying Gaussian PSDF and returns the std of that field.
//
// dKappa[i] must match the spectral discretisation used elsewhere in the
// program: kappaUpper / (N[i] / 2) for the hard-coded kappaUpper = 1.8π.
//
// The algorithm is identical to the 1D case described in SDB-2011 §6; only
// the FFT, the symmetric extension, and the interior-bin iteration become
// dimension-aware.
static double iterateUnderlyingGaussianPsdfND(double *S_G_out,
                                              const double *S_T,
                                              const int N[], int nDim,
                                              const double dKappa[],
                                              const double *cHermite, int Kmax,
                                              int maxIter, double beta,
                                              double tolerance)
{
  int total = 1;
  for (int i = 0; i < nDim; ++i) total *= N[i];

  double dKappaProd = 1.0;
  for (int i = 0; i < nDim; ++i) dKappaProd *= dKappa[i];

  std::vector<double> invFact(Kmax + 1);
  invFact[0] = 1.0;
  for (int k = 1; k <= Kmax; ++k) invFact[k] = invFact[k - 1] / k;

  // ||S_T||² over the interior octant (Eq. 19 denominator).
  double targetNormSq = 0.0;
  for (int j = 0; j < total; ++j) {
    if (linearIndexIsInterior(j, N, nDim)) targetNormSq += S_T[j] * S_T[j];
  }
  if (targetNormSq <= 0.0) {
    std::fprintf(stderr, "iterateUnderlyingGaussianPsdfND: target PSDF is identically zero\n");
    std::exit(1);
  }

  std::vector<double> S_G(total);
  for (int j = 0; j < total; ++j) S_G[j] = S_T[j];

  std::vector<double> fftBuf(2 * total);
  std::vector<double> R_G(total), xi(total), R_NG(total), S_NG(total);

  const double fftToDensity = 1.0 / (static_cast<double>(total) * dKappaProd);

  double prevError = std::numeric_limits<double>::infinity();
  double finalError = prevError;
  int iter = 0;

  for (iter = 0; iter < maxIter; ++iter) {
    // -- Step 2: IFFT(S_G_full) → R_G (unnormalised, NR convention).
    packSymmetricPsdfND(S_G.data(), N, nDim, fftBuf.data());
    multidimFft(fftBuf.data(), N, nDim, +1);
    for (int j = 0; j < total; ++j) R_G[j] = fftBuf[2 * j];

    const double R0 = R_G[0];
    if (R0 <= 0.0) {
      std::fprintf(stderr,
                   "iterateUnderlyingGaussianPsdfND: non-positive variance at iter %d\n",
                   iter);
      std::exit(1);
    }

    // -- Step 3: ξ = R_G / R_G[0], clamped to [-1, 1].
    for (int j = 0; j < total; ++j) {
      double v = R_G[j] / R0;
      if (v >  1.0) v =  1.0;
      if (v < -1.0) v = -1.0;
      xi[j] = v;
    }

    // -- Step 4: Mehler series.
    for (int j = 0; j < total; ++j) {
      double sum = 0.0;
      double xiPow = xi[j];
      for (int k = 1; k <= Kmax; ++k) {
        sum += cHermite[k] * cHermite[k] * xiPow * invFact[k];
        xiPow *= xi[j];
      }
      R_NG[j] = sum;
    }

    // -- Step 5: FFT(R_NG) → S_NG and zero boundary bins.
    for (int j = 0; j < total; ++j) {
      fftBuf[2 * j]     = R_NG[j];
      fftBuf[2 * j + 1] = 0.0;
    }
    multidimFft(fftBuf.data(), N, nDim, -1);
    for (int j = 0; j < total; ++j) {
      S_NG[j] = linearIndexIsInterior(j, N, nDim) ? fftBuf[2 * j] * fftToDensity : 0.0;
    }

    // Renormalise S_NG so its interior integral matches S_T's. Without this,
    // the discrete DC truncation in S_T (set to zero by convention while the
    // Mehler-derived S_NG retains its full κ=0 contribution) causes a global
    // amplitude bias of order ~10% that swamps the per-frequency shape
    // difference the iteration is meant to correct.
    double sumST = 0.0, sumSNG = 0.0;
    for (int j = 0; j < total; ++j) {
      if (linearIndexIsInterior(j, N, nDim)) {
        sumST  += S_T[j];
        sumSNG += S_NG[j];
      }
    }
    if (sumSNG > 0.0) {
      const double rescale = sumST / sumSNG;
      for (int j = 0; j < total; ++j) S_NG[j] *= rescale;
    }

    // -- Step 6: relative error.
    double diffSq = 0.0;
    for (int j = 0; j < total; ++j) {
      if (linearIndexIsInterior(j, N, nDim)) {
        const double d = S_NG[j] - S_T[j];
        diffSq += d * d;
      }
    }
    finalError = 100.0 * std::sqrt(diffSq / targetNormSq);
    std::printf("SDB-2011 iter %3d: ε = %.4f%%\n", iter, finalError);

    if (iter > 0) {
      if (finalError > prevError) break;
      if (prevError - finalError < tolerance) break;
    }
    prevError = finalError;

    // -- Step 7: update S_G over the interior octant.
    for (int j = 0; j < total; ++j) {
      if (!linearIndexIsInterior(j, N, nDim)) continue;
      if (S_NG[j] <= 0.0 || S_T[j] <= 0.0) continue;
      S_G[j] = std::pow(S_T[j] / S_NG[j], beta) * S_G[j];
      if (S_G[j] < 0.0) S_G[j] = 0.0;
    }
  }

  std::printf("SDB-2011 finished: %d iterations, final ε = %.4f%%\n",
              iter + 1, finalError);

  for (int j = 0; j < total; ++j) S_G_out[j] = S_G[j];

  // σ_G² = ∫ S_G dκ = (∏ 2·dκ_i) · Σ over interior octant.
  double sumSG = 0.0;
  for (int j = 0; j < total; ++j) {
    if (linearIndexIsInterior(j, N, nDim)) sumSG += S_G[j];
  }
  double octantFactor = 1.0;
  for (int i = 0; i < nDim; ++i) octantFactor *= 2.0 * dKappa[i];
  return std::sqrt(sumSG * octantFactor);
}


// Copy the corner [0..destShape) of a `srcShape`-sized field into a
// `destShape`-sized buffer. Used to crop a padded simulation grid down to
// the user's requested region.
static void cropField(double *dest, int destShape[],
                      const double *src, int srcShape[], int nDim)
{
  if (nDim == 1) {
    for (int k = 0; k < destShape[0]; ++k) dest[k] = src[k];
  } else if (nDim == 2) {
    for (int k = 0; k < destShape[0]; ++k) {
      for (int l = 0; l < destShape[1]; ++l) {
        dest[k * destShape[1] + l] = src[k * srcShape[1] + l];
      }
    }
  } else if (nDim == 3) {
    for (int k = 0; k < destShape[0]; ++k) {
      for (int l = 0; l < destShape[1]; ++l) {
        for (int m = 0; m < destShape[2]; ++m) {
          const int dIdx = (k * destShape[1] + l) * destShape[2] + m;
          const int sIdx = (k * srcShape[1] + l) * srcShape[2] + m;
          dest[dIdx] = src[sIdx];
        }
      }
    }
  }
}


void generateRandomField(double field[], int numberReal[], double autoLength[],
			 double randomParameter[], int nDim, int typeOfCDF, long randomInteger,
			 int padding, int iterate, int iterMaxIter, double iterTolerance,
			 double iterBeta, int iterHermiteOrder)
{
  // Compute mean and variance from the user's parameters.
  double variance, meanNG;

  // Build the padded shape that the FFT actually runs on. With padding == 1
  // the padded shape equals the user-requested shape and the original code
  // path is taken (the cheap memcpy at the end is a no-op).
  int paddedReal[3] = {1, 1, 1};
  int paddedTotal = 1;
  for (int i = 0; i < nDim; ++i) {
    paddedReal[i] = numberReal[i] * padding;
    paddedTotal *= paddedReal[i];
  }

  double *targetSpectralDensity = new double[paddedTotal];

  if (typeOfCDF == 1) {       // Gaussian
    meanNG   = randomParameter[0];
    variance = std::pow(randomParameter[1] * randomParameter[0], 2.0);
  } else if (typeOfCDF == 2) { // Weibull (analytical mean and variance)
    meanNG   = randomParameter[2];
    variance = randomParameter[3];
  } else if (typeOfCDF == 3) { // Grafted Weibull-Gaussian
    // The grafted CDF has no closed-form unconditional mean or variance
    // (the formula previously used in this slot, (cov · gaussianMean)², is
    // the variance of the Gaussian core only). Compute both numerically by
    // trapezoidal integration of the inverse CDF against the standard
    // normal density.
    computeNonGaussianMomentsByQuadrature(randomParameter, typeOfCDF,
                                          &meanNG, &variance);
  }

  computeTargetSpectralDensity(targetSpectralDensity, paddedReal, autoLength,
                               variance, nDim);

  // Underlying-Gaussian PSDF used by generateGaussianField. Without
  // iteration, this is identical to the user's target PSDF.
  double *underlyingPsdf = targetSpectralDensity;
  std::vector<double> iteratedPsdf;
  double sigmaG = std::sqrt(variance);  // default: σ_G = σ_NG (no iteration)

  // SDB-2011 iteration: search for the underlying Gaussian PSDF whose
  // translation image best matches the user's prescribed target spectrum.
  // Only useful for non-Gaussian distributions; the Gaussian path
  // (typeOfCDF == 1) is the trivial identity transformation. Works in 1D,
  // 2D and 3D using the same Mehler kernel.
  if (iterate && typeOfCDF != 1) {
    double dKappaArr[3] = {0., 0., 0.};
    for (int i = 0; i < nDim; ++i) dKappaArr[i] = (1.8 * MYPI) / (paddedReal[i] / 2);

    std::vector<double> cHerm(iterHermiteOrder + 1, 0.0);
    computeHermiteCoefficients(cHerm.data(), iterHermiteOrder,
                               randomParameter, typeOfCDF, meanNG);
    {
      double sigmaCheck = 0.0; double f = 1.0;
      for (int k = 1; k <= iterHermiteOrder; ++k) {
        f /= k;
        sigmaCheck += cHerm[k] * cHerm[k] * f;
      }
      std::printf("SDB-2011: Hermite c_0 = %.3e (should be near 0); "
                  "Σ c_k²/k! = %.6e (target σ² = %.6e)\n",
                  cHerm[0], sigmaCheck, variance);
    }

    iteratedPsdf.resize(paddedTotal);
    sigmaG = iterateUnderlyingGaussianPsdfND(iteratedPsdf.data(),
                                             targetSpectralDensity,
                                             paddedReal, nDim, dKappaArr,
                                             cHerm.data(), iterHermiteOrder,
                                             iterMaxIter, iterBeta, iterTolerance);
    underlyingPsdf = iteratedPsdf.data();
  } else if (iterate && typeOfCDF == 1) {
    std::printf("SDB-2011 iteration not needed for Gaussian (typeOfCDF=1); ignored.\n");
  }

  // Generate the Gaussian field on the (possibly padded) grid. With
  // padding == 1 we write directly into the caller's buffer; otherwise we
  // use a temporary, then crop the corner the user actually asked for.
  double *gaussianBuffer = (padding == 1) ? field : new double[paddedTotal];
  generateGaussianField(gaussianBuffer, underlyingPsdf, paddedReal,
                        autoLength, variance, randomInteger, nDim);

  if (padding > 1) {
    cropField(field, numberReal, gaussianBuffer, paddedReal, nDim);
    delete[] gaussianBuffer;
  }

  // Translate Gaussian → non-Gaussian. The std passed to mapGaussianToRandom
  // must equal the actual std of the underlying Gaussian field: σ_G after
  // SDB-2011 iteration, otherwise √variance (the value used to scale the
  // target spectrum).
  const double stdOverride = (iterate && typeOfCDF != 1) ? sigmaG : std::sqrt(variance);
  mapGaussianToRandom(field, numberReal, randomParameter, nDim, typeOfCDF, stdOverride);

  printf("Mapping done\n");

  delete[] targetSpectralDensity;
  
  return;
}


// Apply the memoryless translation X = F_NG^{-1}(Φ(g/σ_G)) to a Gaussian
// field in place. The caller specifies σ_G via gaussianStdOverride; when
// that is non-positive (default), the function infers σ from
// randomParameter under the assumption σ_G == σ_NG (correct only without
// SDB-2011 iteration).
void mapGaussianToRandom(double field[], int numberReal[], double randomParameter[],
                         int nDim, int typeOfCDF, double gaussianStdOverride)
{
  int numberRealTotal = numberReal[0];
  if (nDim >= 2) numberRealTotal *= numberReal[1];
  if (nDim >= 3) numberRealTotal *= numberReal[2];

  double standardDeviation;
  if (gaussianStdOverride > 0.0) {
    standardDeviation = gaussianStdOverride;
  } else if (typeOfCDF == 1) {       // Gaussian
    standardDeviation = randomParameter[0] * randomParameter[1];
  } else if (typeOfCDF == 2) {       // Weibull
    standardDeviation = std::sqrt(randomParameter[3]);
  } else if (typeOfCDF == 3) {       // Grafted Weibull-Gaussian
    standardDeviation = std::sqrt(std::pow(randomParameter[3] * randomParameter[5], 2.0));
  } else {
    std::fprintf(stderr, "mapGaussianToRandom: unknown typeOfCDF %d\n", typeOfCDF);
    std::exit(1);
  }

  // For the grafted CDF, precompute the constant grafting length so the
  // inner loop only touches per-sample work.
  double graftingProbability = 0., rescaling = 0., modulus = 0., length = 0.;
  double graftedGaussianMean = 0., graftedGaussianStd = 0.;
  if (typeOfCDF == 3) {
    modulus            = randomParameter[1];
    graftingProbability = randomParameter[2];
    rescaling          = randomParameter[4];
    graftedGaussianMean = randomParameter[5];
    graftedGaussianStd  = randomParameter[6];
    length = pow(-log(1. - graftingProbability / rescaling), 1. / modulus);
  }

  for (int k = 0; k < numberRealTotal; ++k) {
    double gaussianProbability = normal_cdf(field[k], 0., standardDeviation);
    if (typeOfCDF == 1) {           // Gaussian: shift and scale to target.
      field[k] = normal_cdf_inv(gaussianProbability, randomParameter[0], standardDeviation);
    } else if (typeOfCDF == 2) {    // Weibull inverse CDF.
      field[k] = randomParameter[0] * pow(-log(1. - gaussianProbability),
                                          1. / randomParameter[1]);
    } else if (typeOfCDF == 3) {    // Grafted Weibull-Gaussian inverse CDF.
      gaussianProbability /= rescaling;
      if (gaussianProbability < graftingProbability) {
        // Weibull tail.
        field[k] = randomParameter[0] * pow(-log(1. - gaussianProbability),
                                            1. / randomParameter[1]);
      } else {
        // Rescaled Gaussian core.
        const double help = normal_01_cdf((length - graftedGaussianMean) / graftedGaussianStd);
        double helpProbability = gaussianProbability - graftingProbability / rescaling + help;
        if (helpProbability >= 1.) helpProbability = 0.99999;
        field[k] = normal_cdf_inv(helpProbability, graftedGaussianMean, graftedGaussianStd);
        if (field[k] < 0.) {
          std::fprintf(stderr, "mapGaussianToRandom: grafted core produced a negative value\n");
          std::exit(1);
        }
      }
    }
  }
}



// Tabulate the target two-sided spectral density that corresponds to a
// Gaussian autocorrelation kernel R(τ) = σ² ∏_i exp(-(τ_i / b_i)²) with
// b_i = autoLength[i]. Bins where any axis index is 0 or ≥ N_i/2 are zeroed
// — this is the simulator-wide convention and excludes the DC and Nyquist
// components from the field synthesis.
void computeTargetSpectralDensity(double spectralDensityFunction[], int numberReal[],
                                  double autoLength[], double variance, int nDim)
{
  int numberAnal[nDim];
  double deltaKappa[nDim];
  const double kappaUpper = 1.8 * MYPI;

  for (int m = 0; m < nDim; ++m) {
    numberAnal[m] = numberReal[m] / 2;
    deltaKappa[m] = kappaUpper / numberAnal[m];
  }

  if (nDim == 1) {
    for (int k = 0; k < numberReal[0]; ++k) {
      if (k == 0 || k >= numberReal[0] / 2) {
        spectralDensityFunction[k] = 0.;
      } else {
        spectralDensityFunction[k] = variance * autoLength[0] / (2. * sqrt(MYPI))
            * exp(-pow(autoLength[0] * k * deltaKappa[0] / 2., 2.));
      }
    }
  } else if (nDim == 2) {
    for (int k = 0; k < numberReal[0]; ++k) {
      for (int l = 0; l < numberReal[1]; ++l) {
        if (k == 0 || l == 0 || k >= numberReal[0] / 2 || l >= numberReal[1] / 2) {
          spectralDensityFunction[k * numberReal[1] + l] = 0.;
        } else {
          spectralDensityFunction[k * numberReal[1] + l] =
              variance * autoLength[0] * autoLength[1] / (4. * MYPI)
              * exp(-pow(autoLength[0] * k * deltaKappa[0] / 2., 2.)
                    - pow(autoLength[1] * l * deltaKappa[1] / 2., 2.));
        }
      }
    }
  } else if (nDim == 3) {
    for (int k = 0; k < numberReal[0]; ++k) {
      for (int l = 0; l < numberReal[1]; ++l) {
        for (int t = 0; t < numberReal[2]; ++t) {
          const int idx = k * numberReal[1] * numberReal[2] + l * numberReal[2] + t;
          if (k == 0 || l == 0 || t == 0
              || k >= numberReal[0] / 2 || l >= numberReal[1] / 2
              || t >= numberReal[2] / 2) {
            spectralDensityFunction[idx] = 0.;
          } else {
            spectralDensityFunction[idx] =
                variance * autoLength[0] * autoLength[1] * autoLength[2]
                / (8. * pow(MYPI, 1.5))
                * exp(-pow(autoLength[0] * k * deltaKappa[0] / 2., 2.)
                      - pow(autoLength[1] * l * deltaKappa[1] / 2., 2.)
                      - pow(autoLength[2] * t * deltaKappa[2] / 2., 2.));
          }
        }
      }
    }
  } else {
    std::fprintf(stderr, "computeTargetSpectralDensity: nDim must be 1, 2, or 3 (got %d)\n", nDim);
    std::exit(1);
  }
}