Back to the main page Back to the documentation 1. Overview: Grid assignment and interpolation 2. Grid assignment 3. Interpolation 4. Convolutions 5. Grid assignment using interlacing 6. Deconvolution

When working with particles we often want to assign particle to a grid to create a density field in either real or fourier space (useful to have different methods for this as we can save fourier transforms). We also want to be able to interpolate a grid to the particle positions. Here we provide methods for this.

The code provide methods for assigning particles to a grid, i.e. computing the convolution $\delta * W$ with $W$ being an arbitrary order B-spline kernel (or using any user provided kernel). The kernel is defined via $$H^{\rm ORDER} = H * H^{\rm ORDER-1}$$ where $*$ is a convolution and $H$ is the top-hat function $H = 1$ if a particle falls within $1/2$ of a cellsize from the center. In many dimensions the kernel is simply $H = \prod_{i=1}^{\rm NDIM} H_i$. For ORDER$=1$ we get NGP (Nearest Grid Point), ORDER$=2$ is CIC (Cloud in Cell), ORDER$=3$ is TSC (Triangular Shaped Cloud), ORDER$=4$ is PCS (Piecewis Cubic Spline) and ORDER$=5$ is PQS (Piecewise Quartic Spline). A nice property of this kernel is that its Fourier transform is simply $$\left(\prod_{i=1}^{\rm NDIM}\frac{\sin(\frac{\pi}{2}\frac{k}{k_{\rm ny}})}{ (\frac{\pi}{2}\frac{k}{k_{\rm ny}}) }\right)^{\rm ORDER}$$ where $k_{\rm ny} = \pi N_{\rm grid}$ is the Nyquist frequency of the grid (see Deconvolution below for where this is used). Note that we need to visit ${\rm ORDER}^{\rm NDIM}$ cells per particle to assign the density so the cost of the method goes up quite steeply with the order. For most applications CIC is what you want to use as it gives a good compromise between speed and accuracy.

// Assign particle to grid using NGP, CIC, TSC, ... etc.
// The result is the overdensity field rho/rho_mean - 1
template<int N, class T>
  void particles_to_grid(
      T *part,
      size_t NumPart,
      size_t NumPartTot,
      FFTWGrid<N> & overdensity,
      std::string density_assignment_method);

The particle class must have the method get_pos() which gives a pointer to the position of the particle and get_ndim() that tells us how many components the positions have. If your particles have different mass you should also have a get_mass() function. This mass will then be used when we assign particles to a grid. If we don't have this method then we assume that all particles have the same mass. For galaxy surveys the get_mass function should just return the weight of the galaxy. The masses can be arbitrarely normalized (though having them in units of the mean mass is a good convention to use) as we use the mass divided by the mean mass as the weight when assigning to the grid.

We also provide methods for interpolating a grid to particle positions using the same kernels as above (so interpolation_method takes the same values NGP, CIC, TSC etc. as above):

template<int N, class T>
  void interpolate_grid_to_particle_positions(
      const FFTWGrid<N> & grid,
      T *part,
      size_t NumPart,
      std::vector<FloatType> & interpolated_values,
      std::string interpolation_method);

And if you want to convolve the grid with an arbitary kernel you can use:

template<int N, int ORDER, class T>
  void convolve_grid_with_kernel(
      const FFTWGrid<N> & grid_in,
      FFTWGrid<N> & grid_out,
      std::function<FloatType(std::vector<double> &am; dx)> & convolution_kernel);

The argument to the kernel is how far away the point is from the cell in units of the cell-size so the NGP kernel would be $1$ if all $|dx_i| < 0.5$ and zero otherwisee.

Assign particles to grid using the original particles and particles displaced by $1/2$ grid-cell and combine this to get the fourier transform of the density field. This helps reduce aliasing and is useful for computing correlation functions in Fourier space from particles. NB: this method does not decolvolves the window function used for the density assignment so you will have to call deconvolve_window_function_fourier (see next section) if you want to do this:

// Returns density_grid_fourier with delta(k). This grid has to be allocated beforehand
template<int N, class T>
void particles_to_fourier_grid(T * part,
                               size_t NumPart,
                               size_t NumPartTot,
                               FFTWGrid<N> & density_grid_fourier,
                               std::string density_assignment_method,
                               bool interlacing);

When assigning to the grid we don't compute $\delta$, but rather $\delta * W$ where $W$ is the density assignment kernel. In fourier space this convolution becomes a simple product and we can get back the "true" $\delta$ if we want by dividing by the fourier transform of the kernel we use. The method below does this and is used by for example algorithms to compute fourier space correlation functions:

template <int N>
void deconvolve_window_function_fourier(density_grid_fourier, density_assignment_method);