Quickstart

This walks through generating a Gaussian random field with a prescribed power spectrum, then recovering its power spectrum from the realization to check that the two match.

Define a power spectrum

from kspace import PowerLawBetaModel

# A power-law power spectrum with large- and small-scale cutoffs:
#   l_min, l_max set the cutoff scales; alpha is the power-law index.
power_spec = PowerLawBetaModel(l_min=10.0, l_max=200.0, alpha=-11.0 / 3.0)

# Normalize so the field has a given RMS amplitude
power_spec.renormalize(f_rms=10.0)

l_min/l_max and the grid coordinates below are in whatever length unit you choose (kpc is a common choice for the astrophysical use case this package targets, but nothing enforces that).

Generate a field realization

import numpy as np
from kspace import GaussianRandomField

left_edge = np.array([0.0, 0.0, 0.0])
right_edge = np.array([750.0, 750.0, 750.0])
ddims = [256, 256, 256]

grf = GaussianRandomField(left_edge, right_edge, ddims, power_spec, seed=10)
field = grf.generate_scalar_field_realization()

field is a real-valued NumPy array of shape (256, 256, 256).

Recover and check the power spectrum

from kspace import FourierAnalysis

fa = FourierAnalysis(right_edge - left_edge, ddims)
kbins, pk = fa.make_binned_powerspec(field, nbins=60)

# Compare to the input spectrum at the bin centers
k_centers = np.sqrt(kbins[1:] * kbins[:-1])
measured = pk
expected = power_spec(k_centers)

Plotting measured and expected against k_centers on a log-log axis should show good agreement away from the largest and smallest scales (where the finite grid and box size limit the achievable dynamic range).

Forward and inverse FFTs

fa above already wraps the FFT machinery used internally by make_binned_powerspec; you can also call it directly. fftn/ifftn aren’t drop-in replacements for scipy.fft’s functions of the same name — they use the \(k = 2\pi f\) convention and normalize by the cell volume dV (see Conventions), and ifftn requires its input to be the FFTArray that fftn returns, not a plain array:

# Forward transform: a real-space field -> an FFTArray in Fourier space
field_hat = fa.fftn(field)

# Inverse transform: exact inverse of fftn, back to real space
field_recovered = fa.ifftn(field_hat)

np.allclose(field, field_recovered)  # True

fftn also accepts vector fields, of shape (ndim, *ddims), transforming each component along the spatial axes.

Vector fields

# A 3-component vector field, e.g. a synthetic turbulent velocity field
vfield = grf.generate_vector_field_realization()

# Or a divergence-free ("solenoidal") one
vfield_solenoidal = grf.generate_vector_field_realization(divergence_free=True)

Divergence and curl

FourierAnalysis.divergence_of_field/curl_of_field compute these in real space via finite differences. Since vfield came from an FFT-based realization, its domain is periodic, so pass periodic=True to use wraparound (rather than one-sided) differences at the box edges:

div_v = fa.divergence_of_field(vfield, periodic=True)
curl_v = fa.curl_of_field(vfield, periodic=True)

# vfield_solenoidal was constructed to be divergence-free, so this should
# be close to zero everywhere (up to finite-difference truncation error)
div_v_solenoidal = fa.divergence_of_field(vfield_solenoidal, periodic=True)

There are also spectral counterparts, divergence_component/ solenoidal_component, which project a vector field into its compressional and solenoidal parts in Fourier space rather than differentiating in real space — see Conventions for how they relate to divergence_of_field/ curl_of_field and the choice of wavenumber (diff_type) they use.

See Conventions for the Fourier conventions this package uses, and Tutorials for complete, runnable examples.