Skip to content
MPI-parallelised code to compute polyspectrum (including bispectrum) from a 3D grid
C Makefile C++ Shell
Branch: master
Clone or download
Fetching latest commit…
Cannot retrieve the latest commit at this time.
Type Name Latest commit message Commit time
Failed to load latest commit information.
create_iniFile added option to create iniFile with phase option Nov 20, 2019
src removed test output Nov 20, 2019
README.rst edit to readme Jul 22, 2019
bispectra_tests.pdf theta values are also stored in output for bispectrum Aug 24, 2018
iniFile.ini added option of phase only bispectrum, adjusted README Jul 22, 2019



Code to compute polyspectrum from a 3D grid following the method described in Watkinson et al. (2017).

IMPORTANT:: The current input file only supports powerspectra and bispectra. Test results analogue to Fig. 7 in Watkinson et al. (2017) can be found here. Suggestions for improvements (particularly run time) are very welcome and can be sent to



Serial run

  1. fftw3 library: fftw3 >= 3.3.3

Parallel run

  1. MPI library
  2. fftw3 & fftw3-mpi library: fftw3 >= 3.3.3

Go to the FFTW webpage to install fftw3. Ensure to compile the library with the enable-mpi flag for parallel runs

$ ./configure --enable-mpi
$ make
$ make install

Note: To create the dynamic libraries, run configure with the --enable-shared flag.

Download & Build

$ git clone
$ make

This will download the code and first test case from the github directory and compile the source code.


The first test case can then be run by

$ ./polyspectrum iniFile.ini

iniFile.ini contains all input parameters that are needed for any runs. For a different simulation the code does not need to be recompiled but just this parameter file iniFile.ini to be adapted.

Parameter file

The following parameters are specified in iniFile.ini:

  • gridsize: size of the 3D grid along one axis
  • boxsize: comoving boxsize in Mpc/h
  • gasInputsInDoublePrecision: set to 0 for single, 1 for double precision for density fields being read
  • ionInputsInDoublePrecision: set to 0 for single, 1 for double precision for ionization fields being read
  • densityFile: path to 3D density grid
  • ionFile: path to 3D ionization grid
  • hubble_h: H = 100*hubble_h km/s/Mpc
  • omega_b: baryon density parameter
  • omega_m: matter density parameter
  • omega_l: lambda density parameter
  • Y: mass fraction of Helium in the primordial gas (assumed to consist of H and He)
  • whichField: field from which the polyspectrum is calculated; options are: DENS for density, XHII for ionization fraction and XHI_DENS for neutral gas density
  • useOnlyPhase: 1 yields the phase-only bispectrum
  • n: number of vectors, i.e. n=2 yields powerspectrum, n=3 yields bispectrum
  • equilateral: 1 yields bispectrum for equilateral triangles with numValues providing the number of k-values depending on the box and gridsize of the simulation box, for all values not equal to 1 non-equilateral triangles are computed
  • k1: length of first vector in Mpc/h
  • k2: length of second vector in Mpc/h (if n>2)
  • numValues: set to 1, if only single value for polyspectrum should be calculated; for >1 and n>2 number of k bins along theta = 0 to 180
  • theta: angle between k1 and k2 in rad
  • kbinwidth: bin width in k-space in units of grid cells for computing the bispectra
  • kbinningCase: DEFAULT assumes for all k vectors a binwidth in k-space of kbinwidth, GIVEN_BINNING assumes for k1 and k2 a binwidth in k-space of kbinwidth and for k3 a width that corresponds to cosThetaBinwidth, DERIVED_BINNING assumes for k1 and k2 a binwidth in k-space of kbinwidth and for k3 it derives the binwidth from the allowed ranges given by the uncertainties of k1 and k2
  • cosThetaBinwidth: binwidth in cos(theta)
  • output_dir: directory where output should be written
  • output_basename: output name of the run
  • write_numpolygons: 1 if output files should contain number of polygons, otherwise they are not included in the output
You can’t perform that action at this time.