A workshop dedicated to the Theoretical Virtual Observatory will take place at IAP on April 5-6th.
The goal is to bring together experts of the Virtual Observatory and theoreticians who would like to make results of their simulations (e.g. databases or catalogs) or numerical codes available to the worldwild astronomical community.
POWMES is a F90 program to measure very accurately the power spectrum in a N-body simulation, using Taylor expansion of some order on the cosine and sine transforms. It can read GADGET format and requires FFTW2 to be installed.
The code itself is found here
An extended version with an example is found there
The corresponding paper is found there
! Authors : Stephane Colombi (Institut d'Astrophysique de Paris,
! Dmitri Novikov (Imperial College, London,
! Publication : Colombi, Jaffe, Novikov, Pichon, MNRAS, in press
! The POWMES software is provided `as is', and the authors (SC & DN)
! make no representations or warranties, expressed or implied, of
! merchantability or fitness for any particular purpose, and disclaim
! all liability for direct or consequential damages resulting from
! anyone's use of the program.
! THE PROGRAM
! This programs measures the power spectrum in a 3D simulation using
! Fourier Taylor expansion of some order on the cosine and sine
! This program is still under development, and deserves further
! optimization to deal with large simulation samples.
! Taylor expansion :
! Approximation of zeroth order provides NGP interpolation (Nearest Grid
! Point) while higher order approximations take into account small
! displacements within cells of the mesh used to perform the Fourier
! transform increasingly accurately. Arbitrarily high order converges
! to the exact transform, namely to the sum
! delta_k=\Sum_i M_i exp(I.k.x_i/N), (1)
! where x_i is the position of the particle in [0,2.PI[, M_i its mass,
! (total mass of the system unity), k is an integer wavenumber, I^2=-1,
! and N is the size of the grid used to perform the calculations.
! The power spectrum P_k reads
! P_k=< |delta_k|^2 >_k,
! where <.>_k represents an angular average over the values of k. More
! exactly, function
! defines the size of the bin in k space (all k having the same E_k
! contribute to the same bin).
! Folding :
! To speed up the calculation, it is possible to perform foldings of the
! particle distribution. One folding along one axis consists in replacing
! all the positions veryfing x_i >= PI with x_i-PI, and then multiplying
! all the x_i's (including the unmodified ones) by a factor 2 to span
! again the range [0,2.PI[. This operation is performed on all the
! coordinates. One can then demonstrate that the modes given by
! equation (1) remain unchanged for even values of k.
! Assuming now that we have a periodic simulation box and that we
! use a grid of fixed size N to perform the calculations, successive
! foldings will allow us to estimate the power spectrum for
! 0 folding : k=0,1,2,3,4,...,N/2
! 1 folding : k=2,4,6,8,...,N
! 2 foldings : k=4,8,12,16,...,2N
! 4 foldings : k=8,16,24,32,...,4N
! and so on.
! Errors :
! There are 2 possible sources of errors in the calculation.
! (a) the errors due to the finite number of modes sampled in a given
! bin k. If N_k is the number of independent modes, then the relative
! error due to the finiteness of N_k is approximately
! Delta P_k/P_k=1/sqrt(N_k). (2)
! Obviously Delta P_k/P_k decreases while k increases, since it
! scales roughly like 1/k.
! In the code, the number of independent modes, N_k, is output.
! Another way of estimating the errors is also performed, by
! simple estimation of the scatter in each bin. It turns out in
! practice to give errors of the same order of equation (2), which
! gives the results for an underlying random Gaussian field.
! Actual calculation of the errors involve non Gaussian terms,
! so the errors provided by the code are only approximations that
! should be improved on for really accurate estimates of errorbars.
! (b) the biases due to the Taylor expansion and the discreteness of
! the particle distribution, namely
! (i) the so called 1/N_part shot noise bias where N_part is the number
! of particles. In the Fourier Taylor method, this bias is
! actually slightly different from 1/N_part, but this is estimated
! accurately by the code. Note that for a slightly perturbed grid
! pattern, such a bias does not arise and should not be corrected
! for. Only in a relaxed, locally Poissonian stage, the correction
! for shot noise makes really sense.
! (ii) the bias due to the interpolation involved in the calculation.
! It depends on the order Norder of the Taylor expansion.
! For example, for NGP, Norder=0, the bias is given by the
! famous sinc^2(k/2) term. This bias can be easily corrected
! for at any order. It is performed in the code.
! (iii) The unknown effects of aliasing of the contributions of the
! power-spectrum at scales smaller than the sampling imposed
! by the grid used to perform the calculations. In the limit
! of a locally Poissonian stationary random process, these effects
! of aliasing can be bounded. The code computes useful numbers
! to estimate the maximum effect of aliasing (see Colombi et
! al. paper for details). The effects of aliasing increase when
! approaching the nyquist frequency of the sampling grid, so
! it is wise to stay sufficiently away from the Nyquist frequency
! to avoid too much aliasing.
! Obviously, by combining foldings and using (a)+(b), one can find
! the best compromise that minimizes the error on the estimate of P_k for
! a given k at a minimum computational cost. However this is not performed
! at the present stage.
! In its current form, the code computes P(k) for each folding of the particle
! distribution. It computes as well shot noise contribution (b-i), and
! proposes a correction for the biase (b-ii) and usefull numbers for
! estimating (b-iii). Of course, the statistical errors are also estimated,
! (a). It outputs, for each folding, detailed data. However, for practical
! use, a global file contains all the combined information of all the
! foldings and this is what really matter for the user who does not
! necessarily want to have fully optimal results, but still fairly
! accurate ones. Right now, only that part of the code is fully explained
! in what follows.
! Goto the directory /bin and modify the Makefile according to the
! instructions in it, type make.
! Notes :
! - The program needs FFTW2, that has to be installed.
! - It is possible to use openMP acceleration on shared memory computers,
! in that case the thread part of FFTW2 has to be properly installed.
! - Possible problems of underscore while calling FFTW2, follow the
! instructions in the Makefile.
! - Be aware of the little/big endians problems for input file reading.
! Create a config file, following the example given in powmes.config, say
! myconfigfile.config. Then type the command
! whereispowmes/powmes myconfigfile.config, where whereispowmes
! is the location of the executable.
! Typing just the command powmes assumes that the name of the config file
! is powmes.config.
! PARAMETERS IN THE CONFIG FILE
! verbose (.true./.false.) : verbose mode
! megaverbose (.true./.false.) : detailed verbose mode
! filein ('myfile') : name of the input file
! nfile (1/2/3) : type of the input file
! 1 : GADGET file
! 2 : RAMSES file
! 3 : 4 columns ascii file :
! first line : npart where npart is the number of particles
! next lines : x(i),y(i),z(i),mass(i) with x,y,z in [0,1[
! nmpi (-1/positive integer) : number of MPI files for nfile=1
! (this parameter is thus ignored for nfile=2 or nfile=3)
! set nmpi=-1 if there is only one file
! read_mass (.true./.false.) : .true. for reading the mass in the
! file, .false. for ignoring it and giving the same mass to
! all the particles. The parameter read_mass is ignored
! if nfile=1 or nfile=2, where all the masses are supposed
! to be the same for all the particles. WARNING : the code
! has not been tested for unequal masses : this is still under
! nfoldpow (integer) : folding parameter
! - if nfoldpow >= 0, it corresponds to the number
! of foldings to be performed. Note that the largest value of k
! that will be probed will be kmax=k_nyquist*2^(nfoldpow-1), where
! k_nyquist is the nyquist frequency of the grid used to perform the
! - if nfoldpow < 0, then kmax=-nfoldpow is the maximum integer
! wavenumber up to which one wants to measure the power-spectrum.
! The number of foldings needed to achieve that is computed
! ngrid (even positive integer): resolution of the grid used to
! perform the calculations.
! ngrid must be an even number (but not necessarily a power of 2).
! norder (positive integer): order of the Taylor expansion.
! Norder=3 is the recommended value.
! shift(float,float,float) : vector corresponding to a constant shift to
! be applied to all the particles prior to the calculation. This
! can be handy to improve the accuracy of the measurements for a
! slightly perturbed grid. shift is expressed in units of the grid
! cell size.
! filepower ('powspec.dat') : the name of the main output file where all
! the informations of the measured power-spectrum are gathered
! in a handy (but not necessarily optimal in terms of signal
! to noise or biases) way.
! This is a 6 columns ascii file :
! column 1 : the integer wave number
! column 2 : the number of statistically independent modes C(k)
! contributing to the bin k. This number can be useful to
! compute errorbars : the statistical error
! on the rough power-spectrum
! is, under the Gaussian field hypothesis,
! DP_rough(k)/P_rough(k)=1/sqrt(C(k)) (3)
! column 3 : the estimated rough power spectrum PP_rough(k) before
! debiasing from Taylor interpolation.
! column 4 : the debiased estimated P_rough(k). To obtain the true
! power-spectrum, one can subtract the white noise contribution
! (except for a slightly deformed grid, e.g. initial conditions)
! P(k)=P_rough(k)-W(k)/N_part, (4)
! where W(k) is given by column 5 and is a number very
! close to unity.
! column 5 : the shot noise factor W(k) very close to unity, needed
! to correct for shot noise contribution the power spectrum
! given by equation (4).
! column 6 : the statistical error computed self-consistently from
! the dispersion inside each k bin.
! The number displayed is DP_rough(k)/P_rough(k) and
! should give in practice something very similar to
! equation (3).
! filepower_fold ('powspec'/'#junk') : root name of the output files where
! the measurements are stored separately for each folding. If the
! name of the file starts with a '#", these optional files are not
! Example : if filepower_fold='powspec', then 4 files are created :
! powspec.waven : the wavenumbers (1st column of filepower)
! powspec.nmodes : the number of modes file (second column of filepower)
! powspec.power : the rough power-spectrum, prior to debiasing
! (third column of filepower)
! powspec.powerdebiased : the debiased rough power-spectrum
! (fourth column of filepower)
! powspec.shotfac: the shot noise factor W(k), needed to correct
! for shot noise contribution when relevant
! (5th column of filepower)
! powspec.staterr: the statistical error (6th column of filepower)
! More specifically :
! powspec.waven : an ascii file with nfoldpow+1 columns, each column
! corresponding to a given folding (first column : no fold,
! second column : 1 fold, etc). In each column, the integer
! wavenumber probed by the calculation is written.
! Example : for ngrid=8 and nfoldpow=3 powspec.waven
! looks as follows :
! 0 0 0 0
! 1 2 4 8
! 2 4 8 16
! 3 6 12 24
! 4 8 16 32
! powspec.nmodes : similar as powspec.waven, but the number of modes
! in the corresponding k-bin is output
! powspec.power : similar as powspec.waven, but the measured value of
! the rough power-spectrum prior to debiasing is written in
! each column.
! Similarly for powspec.powerdebiased, powspec.shotfac,powspec.staterr.
! filetaylor ('powspec.taylor'/'#junk') : name of the file containing
! useful informations for the fourier-taylor calculation. This file
! is optional and is not output if the name starts with a '#'.
! This is a multicolumn ascii file, each line corresponds to an
! integer value of k, starting from 0.
! Column 1 : the number of independent modes C(k)
! Column 2 : an estimate of L(k)=W(k)*U(k) where U(k) is explained
! Column 3 : approximation for L(k) valid at small k
! Column 4 : the estimated bias U(k) on the rough power spectrum
! Basically, debiasing the power-spectrum (i.e. passing
! from column 3 to column 4 of filepower) consists in
! multiplying the rough power-spectrum by 1/U(k)
! Column 5 : The W(k) already mentionned earlier.
! For more details on filetaylor, ask to colombATiap.fr