Conventions#

kspace makes a few specific choices about Fourier conventions and grid layout that matter if you’re combining its output with your own FFT code, or extending the package itself.

The \(2\pi\) wavenumber convention#

Wavenumbers are defined as \(k = 2\pi f\), where \(f\) is the ordinary (cycles-per-unit-length) frequency from numpy/scipy’s fftfreq. This is the convention used throughout physics for spatial power spectra. Concretely:

k = 2 * np.pi * np.fft.fftfreq(n, d=delta)

and correspondingly dk = 2π / width for the spacing between discrete wavenumbers in a box of a given width. If you generate your own arrays to feed into FourierAnalysis methods, make sure they use this same convention — mixing \(k\) and \(f\) conventions will silently give wrong amplitudes and mode locations.

Normal vs. modified wavenumbers#

Most of kspace — power spectra, fftn/ifftn, GRF generation — works with the “normal” wavenumbers described above: the exact FFT frequencies scaled by \(2\pi\), with no dependence on how a derivative might be discretized.

A few methods that compute spatial derivatives spectrally (divergence, curl, vector potentials — see generate_waves() and the diff_type argument below) instead offer modified wavenumbers: \(k\) values altered so that multiplying by \(ik\) in Fourier space reproduces a specific real-space finite-difference stencil exactly, rather than the exact spectral derivative. kspace calls this choice diff_type:

  • "continuum": the normal, unmodified FFT wavenumbers — an exact spectral derivative.

  • "central": wavenumbers rescaled to \(\sin(k\,dx)/dx\), matching a periodic central-difference stencil.

  • "forward": wavenumbers rescaled to a (complex-valued) periodic forward-difference stencil.

Modified wavenumbers exist so that a spectral derivative can be made to agree with a real-space finite-difference calculation done elsewhere in a pipeline (e.g. on a simulation grid); reach for "continuum" unless you specifically need that agreement. See Vector potentials and curl inversion below for the full set of caveats — in particular, "forward" wavenumbers are complex and change which dot product (k·k vs. k·k̄) is correct in downstream formulas.

FFT normalization#

FourierAnalysis.fftn() and ifftn() don’t use scipy.fft’s default normalization. They multiply/divide by the cell volume dV so that the discrete transform approximates a continuum Fourier transform:

\[ \hat{f}(\mathbf{k}) \approx \sum_n f(\mathbf{x}_n)\, e^{-i\mathbf{k}\cdot\mathbf{x}_n} \, dV \]

Don’t mix raw scipy.fft calls with these wrapped methods on the same data — the amplitudes won’t match.

FFTArray#

FFTArray is a thin ndarray subclass that carries the grid spacing (delta) as metadata, produced by fftn/make_powerspec/etc. Some FourierAnalysis methods (ifftn, integrate_kspace) require an FFTArray specifically and will raise TypeError on a plain array, so that operations aren’t silently applied to data on an incompatible grid.

Power, energy, and amplitude spectra#

A PowerSpectrum (PowerLaw, PowerLawBetaModel, DoublePowerLaw) parametrizes the isotropic power spectrum \(P(k)\) directly: calling power_spec(k) evaluates \(P(k)\). Two related spectra are also available, and both depend on ndim (\(d\)) because they involve integrating \(P(k)\) over the surface of a \(k\)-space sphere of radius \(k\):

  • E() — the energy spectrum, \(P(k)\) integrated over that shell, so that \(\int E(k)\,dk\) over some range of \(k\) gives the variance contributed by that range of scales:

    \[\begin{split} E(k) = \frac{1}{(2\pi)^d} \times \begin{cases} P(k) & d = 1 \\ 2\pi k\, P(k) & d = 2 \\ 4\pi k^2\, P(k) & d = 3 \end{cases} \end{split}\]

    The \(2\pi k\) and \(4\pi k^2\) factors are the circumference/surface area of the shell (in 1D the “shell” is just the two points \(\pm k\), so there’s no extra geometric factor). renormalize() uses this: it rescales \(P(k)\)’s normalization so that \(\int E(k)\,dk = f_{\rm rms}^2\) over the requested range.

  • A() — the amplitude spectrum, \(A(k) = \sqrt{E(k)\,k}\), the Fourier amplitude associated with wavenumber \(k\).

\(P(k)\) — not \(E(k)\) or \(A(k)\) — is the quantity that make_binned_powerspec() recovers from a field realization; see the next section.

Gridded vs. binned power spectra#

make_powerspec() and make_binned_powerspec() both compute a power spectrum from a field, but at different levels of reduction:

  • make_powerspec(data) returns \(P(\mathbf{k}) = |\hat{f}(\mathbf{k})|^2 / V\) on the full Fourier grid — one value per discrete mode \(\mathbf{k}\), as an FFTArray of the same shape as the (transformed) input. It makes no isotropy assumption, so it’s the right starting point if you need per-mode or directional power rather than a single curve — make_binned_powerspec() and integrate_kspace()-based calculations both build on it.

  • make_binned_powerspec(data, bins) calls make_powerspec internally, then histograms the gridded values by \(|\mathbf{k}|\) into 1-D bins, averaging over all modes in each shell — i.e. it assumes isotropy. The result is directly comparable to power_spec(k) on a PowerSpectrum instance (see Power, energy, and amplitude spectra above), which is exactly \(P(k)\), not \(E(k)\) or \(A(k)\) — as in the Quickstart power-spectrum-recovery example. Bins with no grid modes in range come back masked (numpy.ma.masked_invalid) rather than as 0 or NaN.

Divergence-free vector fields#

generate_vector_field_realization(divergence_free=True) projects out the longitudinal (along-\(\mathbf{k}\)) component of each Fourier mode:

\[ \mathbf{B}(\mathbf{k}) \mathrel{-}= \hat{\mathbf{k}}\,(\hat{\mathbf{k}} \cdot \mathbf{B}(\mathbf{k})) \]

then rescales the result by \(\sqrt{n/(n-1)}\) (\(n\) = 2 or 3 dimensions) to restore the per-component power lost by the projection, so the transverse field still has the same power spectrum \(P(k)\) as an unprojected component.

Vector potentials and curl inversion#

Given a divergence-free vector field \(\mathbf{B}\), potential_of_field() recovers a vector potential \(\mathbf{A}\) with \(\mathbf{B} = \nabla \times \mathbf{A}\) in the Coulomb gauge (\(\nabla \cdot \mathbf{A} = 0\)), using

\[ \hat{\mathbf{A}}(\mathbf{k}) = \frac{i\,\mathbf{k} \times \hat{\mathbf{B}}(\mathbf{k})}{\mathbf{k}\cdot\mathbf{k}} \]

In 2D, \(\mathbf{k}\) and \(\hat{\mathbf{B}}\) are both 2-vectors, so their cross product collapses to a scalar — \(\mathbf{A}\) is just that out-of-plane component.

This inversion, and the corresponding curl_of_field() / divergence_component(), accept a diff_type argument selecting which wavenumbers to use in the underlying formula (see generate_waves()):

  • "continuum" (default): the exact FFT wavenumbers. This is the exact algebraic inverse of a curl computed the same way — not of curl_of_field(), which differentiates in real space via numpy.gradient and only agrees with the spectral result at low \(k\).

  • "central": wavenumbers modified to match a real-space central-difference stencil (\(\sin(k\,dx)/dx\)). These remain real-valued, so the formulas above apply unchanged.

  • "forward": wavenumbers modified to match a forward-difference stencil. These are complex-valued, and the \(\mathbf{k}\cdot\mathbf{k}\) above is the plain (unconjugated) self-dot-product \(\sum k_i^2\) — not \(\sum |k_i|^2\) — since that’s what the underlying vector identity \(\mathbf{k}\times(\mathbf{k}\times\mathbf{B}) = \mathbf{k}(\mathbf{k}\cdot\mathbf{B}) - (\mathbf{k}\cdot\mathbf{k})\mathbf{B}\) requires. The two coincide for the real-valued "continuum"/"central" cases but not for "forward".