Power Spectra
=============

Though SWIFT includes functionality for calculating power spectra, this tool
runs on-the-fly, and as such after a run has completed you may wish to create
a number of more non-standard power spectra.

These tools are available as part of the :mod:`swiftsimio.visualisation.power_spectrum`
module. Making a power spectrum consists of two major steps: depositing the particles
on grid(s), and then binning their fourier transform to get the one-dimensional power.


Depositing on a grid
--------------------

Depositing your particles on a grid is performed using
:func:`swiftsimio.visualisation.power_spectrum.render_to_deposit`. This function
performs a nearest-grid-point (NGP) of all particles in the provided particle
dataset. For example:

.. code-block:: python

    from swiftsimio import load
    from swiftsimio.visualisation.power_spectrum import render_to_deposit

    data = load("cosmo_volume_example.hdf5")

    gas_mass_deposit = render_to_deposit(
        data.gas,
        resolution=512,
        project="masses",
        parallel=True,
    )

The specific field being depositied can be controlled with the ``project``
keyword argument. The ``resolution``` argument gives the one-dimensional
resolution of the 3D grid, so in this case you would recieve a ``512x512x512``
grid. Note that the ``gas_mass_deposit`` is a :obj:`~swiftsimio.objects.cosmo_array`,
and as such includes cosmological and unit information that is used later
in the process.


Generating a power spectrum
---------------------------

Once you have your grid deposited, you can easily generate a power spectrum
using the
:func:`~swiftsimio.visualisation.power_spectrum.deposition_to_power_spectrum`
function. For example, using the above deposit:

.. code-block:: python

    from swiftsimio.visualisation.power_spectrum import deposition_to_power_spectrum

    wavenumbers, power_spectrum, _ = deposition_to_power_spectrum(
        deposition=gas_mass_deposit,
        boxsize=data.metadata.boxsize,
    )

This power spectrum can then be plotted. Units are included on both the wavenumbers
and the power spectrum. Cross-spectra are also supported through the
``cross_deposition`` keyword, but by default this generates the auto power.
Wavenumbers are calculated to be at the weighted mean of the k-values in each
bin, rather than representing the center of the bin.


More complex scenarios
----------------------

In a realistic simualted power spectrum, you will need to perform 'folding'
to achieve a viable dynamic range within an achievable memory footprint.
Consider, for instance, a 1 Gpc simulation volume with a 1 kpc resolution
limit. In that case, you would need a deposition grid with :math:`10^{18}`
cells, amounting to a 4 EB (yes, exabyte!) memory footprint. That is,
of course, not realistic!

As in a power spectrum we are only interested in the periodicity of the
system, we can fold it back in on itself during the rendering process.
The position of the particle in the box is set to be:

:math:`x_i' := \left( x_i \% \frac{L}{2^{n}} \right) \frac{2^{n}}{L}`

where :math:`L` is the box-size and :math:`n` is some integer greater
than or equal to zero. This allows you to probe modes in the reduced
box-length :math:`L / 2^{n}` with the same fixed resolution deposition
buffer.

The ``folding`` parameter is available for both ``render_to_deposit``
and ``deposition_to_power_spectrum``, but it may be easier to use the
utility functions provided for automatically stitching together
the folded spectra. The function
:func:`~swiftsimio.visualsation.power_spectrum.folded_depositions_to_power_spectrum`
allows you to do this easily:

.. code-block:: python

   import unyt as u
   from swiftsimio.visualisation.power_spectrum import folded_depositions_to_power_spectrum
   from swiftsimio.objects import cosmo_array

   folded_depositions = {}

   for folding in [x * 2 for x in range(5)]:
       folded_depositions[folding] = render_to_deposit(
           data.gas,
           resolution=512,
           project="masses",
           parallel=True,
           folding=folding,
       )

   bins, centers, power_spectrum, foldings = folded_depositions_to_power_spectrum(
       depositions=folded_depositions,
       boxsize=data.metadata.boxsize,
       number_of_wavenumber_bins=128,
       wavenumber_range=cosmo_array(
           [1e-2, 1e2],
           u.Mpc**-1,
           comoving=True,
           scale_factor=data.metadata.a,
           scale_exponent=-1,
       ),
       log_wavenumber_bins=True,
       workers=4,
       minimal_sample_modes=8192,
       cutoff_above_wavenumber_fraction=0.75,
       shot_noise_norm=len(gas_mass_deposit),
   )

The 'used' foldings of the power spectrum are shown in the
``foldings`` return vaule, which is an array containing the folding
that was used for each given bin. This is useful for debugging and
visualisation.

There are a few crucial parameters to this function:

1. ``workers`` is the number of threads to use for the calculation of
   the Fourier transforms.
2. ``minimal_sample_modes`` is the minimum number of modes that must be
   present in a bin for it to be included in the final power spectrum.
   Generally for a big simulation you want to set this to around 10,000,
   and this number is ignored for the lowest wavenumber bin.
3. ``cutoff_above_wavenumber_fraction`` is the fraction of the
   individual fold's (as represented by the FFT itself) maximally sampled
   wavenumber. Ignored for the last fold, and we always cap the maximal
   wavenumber to the Nyquist frequency.
4. ``shot_noise_norm`` is the number of particles in the simulation
   that contribute to the power spectrum. This is used to normalise
   the power spectrum to the shot noise level. This is very
   important in this case because of the use of NGP deposition.

Foldings are stitched using a simple method where the 'better sampled'
foldings are used preferentially, up to the cutoff value.