diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 2086d52cd46df3d8ee41daaafe58b6e1fc0cd50e..830f676a4f00c11a447c7eb87ad541be2a588a3a 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -1280,6 +1280,7 @@ The options are:
* The number of grid foldings to use: ``num_folds``.
* The factor by which to fold at each iteration: ``fold_factor`` (default: 4)
* The order of the window function: ``window_order`` (default: 3)
+ * Whether or not to correct the placement of the centre of the k-bins for small k values: ``shift_centre_small_k_bins`` (default: 1)
The window order sets the way the particle properties get assigned to the mesh.
Order 1 corresponds to the nearest-grid-point (NGP), order 2 to cloud-in-cell
@@ -1319,6 +1320,20 @@ are mutually exclusive. The particles selected in each half are different in
each output. Note that neutrino PS can only be computed when neutrinos are
simulated using particles.
+SWIFT uses bins of integer :math:`k`, with bins :math:`[0.5,1.5]`, :math:`[1.5,2.5]` etc. The
+representative :math:`k` values used to be assigned to the bin centres (so k=1, 2, etc), which are
+then transformed to physical :math:`k` by a factor :math:`kL/(2*pi)`. For the first few bins, only
+few modes contribute to each bin. It is then advantageous to move the "centre" of the bin to the
+actual location correponding to the mean of the contributing modes. The :math:`k` label of the bin
+is thus shifted by a small amount. The way to calculate these shifts is to consider a 3D cube of
+:math:`(kx,ky,kz)` cells and check which cells fall inside a spherical shell with boundaries
+:math:`(i+0.5,i+1.5)`, then calculate the average :math:`k=\sqrt{kx^2+ky^2+kz^2}`. So for
+:math:`i=0` there cells :math:`k=1` and 12 cells :math:`k=\sqrt(2)`, so the weighted k becomes
+:math:`(6 * 1 + 12 * \sqrt(2)) / 18 = 1.2761424`. Note that only the first 7 (22) bins require a
+correction larger than 1 (0.1) percent. We apply a correction to the first 128 terms. This
+correction is activated when ``shift_centre_small_k_bins`` is switched on (that's the default
+behaviour).
+
An example of a valid power-spectrum section of the parameter file looks like:
.. code:: YAML
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index f379b938a18b370a6246bac5ca7d3a6944a4bdff..232eed341614884e687a62ddcc8fecdd42b24b7a 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -932,6 +932,7 @@ PowerSpectrum:
num_folds: 6 # Number of foldings (1 means no foldings), determines the max k
fold_factor: 4 # (Optional) factor by which to reduce the box along each side each folding (default: 4)
window_order: 3 # (Optional) order of the mass assignment scheme (default: 3, TSC)
+ shift_centre_small_k_bins: 1 # (Optional) Correct the centre of the bins with a small k to account for the small number of modes entering the bin.
output_list_on: 0 # (Optional) Enable the output list
output_list: ./output_list_ps.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
requested_spectra: ["matter-matter","cdm-cdm","starBH-starBH","gas-matter","pressure-pressure","matter-pressure", "neutrino0-neutrino1"] # Array of strings indicating which components should be correlated for power spectra
diff --git a/src/power_spectrum.c b/src/power_spectrum.c
index 08189a25250179174668f5fe08cda4ea2764904d..37c6cd0025b8bb8fc5ac727e7e9b4e27d80be452 100644
--- a/src/power_spectrum.c
+++ b/src/power_spectrum.c
@@ -51,6 +51,37 @@
#define power_data_default_fold_factor 4
#define power_data_default_window_order 3
+/* The way to calculate these shifts is to consider a 3D cube of (kx,ky,kz)
+ * cells and check which cells fall inside a spherical shell with boundaries
+ * (i+0.5,i+1.5), then calculate the average k=sqrt(kx^2+ky^2+kz^2). So for i=0
+ * you'd find 6 cells k=1 and 12 cells k=sqrt(2), so the weighted k becomes
+ * (6 * 1 + 12 * sqrt(2)) / 18 = 1.2761424 – etc.
+ * Note that beyond the 7th term, the correction is < 1%. */
+#define number_of_corrected_bins 128
+static const float correction_shift_k_values[number_of_corrected_bins] = {
+ 1.2761424f, 1.1154015f, 1.0447197f, 1.0151449f, 1.0195166f, 1.0203214f,
+ 1.0102490f, 1.0031348f, 1.0063766f, 1.0093355f, 1.0055681f, 1.0024279f,
+ 1.0034435f, 1.0038386f, 1.0011069f, 1.0002888f, 1.0018693f, 1.0029172f,
+ 1.0019128f, 1.0009282f, 1.0015312f, 1.0016361f, 1.0009436f, 1.0003777f,
+ 1.0005931f, 1.0010948f, 1.0010581f, 1.0009779f, 1.0010282f, 1.0008224f,
+ 1.0006637f, 1.0004002f, 1.0002419f, 1.0005172f, 1.0005523f, 1.0004342f,
+ 1.0005183f, 1.0005357f, 1.0003162f, 1.0001836f, 1.0003737f, 1.0004792f,
+ 1.0004169f, 1.0003660f, 1.0004468f, 1.0004218f, 1.0001436f, 1.0000479f,
+ 1.0002012f, 1.0003710f, 1.0003234f, 1.0002661f, 1.0003446f, 1.0003313f,
+ 1.0001844f, 1.0000630f, 1.0001714f, 1.0002382f, 1.0001507f, 1.0001663f,
+ 1.0002199f, 1.0002403f, 1.0000911f, 0.9999714f, 1.0001136f, 1.0001907f,
+ 1.0001917f, 1.0001684f, 1.0001875f, 1.0002158f, 1.0000941f, 1.0000646f,
+ 1.0000930f, 1.0001497f, 1.0001589f, 1.0001215f, 1.0001563f, 1.0001254f,
+ 1.0000557f, 1.0000220f, 1.0000517f, 1.0001039f, 1.0001185f, 1.0000778f,
+ 1.0000848f, 1.0001415f, 1.0001108f, 1.0000709f, 1.0000724f, 1.0001201f,
+ 1.0001480f, 1.0001204f, 1.0001185f, 1.0000844f, 1.0000224f, 0.9999752f,
+ 0.9999997f, 1.0000969f, 1.0001076f, 1.0000756f, 1.0000700f, 1.0000854f,
+ 1.0001067f, 1.0000390f, 1.0000443f, 1.0000863f, 1.0000585f, 1.0000352f,
+ 1.0000677f, 1.0001081f, 1.0000537f, 1.0000199f, 1.0000308f, 1.0000585f,
+ 1.0000479f, 1.0000304f, 1.0000751f, 1.0000710f, 1.0000152f, 1.0000083f,
+ 1.0000342f, 1.0000530f, 1.0000543f, 1.0000442f, 1.0000680f, 1.0000753f,
+ 1.0000369f, 1.0000117f};
+
#ifdef HAVE_FFTW
/**
@@ -1154,7 +1185,9 @@ void power_spectrum(const enum power_type type1, const enum power_type type2,
else
sprintf(powunits, "Mpc^3 eV cm^(-3)");
- fprintf(outputfile, "# Folding %d, all lengths/volumes are comoving\n",
+ fprintf(outputfile,
+ "# Folding %d, all lengths/volumes are comoving. k-bin centres "
+ "are not corrected for the weights of the modes.\n",
i);
fprintf(outputfile, "# Shotnoise [%s]\n", powunits);
fprintf(outputfile, "%g\n", shot);
@@ -1208,8 +1241,17 @@ void power_spectrum(const enum power_type type1, const enum power_type type2,
power_init_output_file(outputfile, type1, type2, us, phys_const);
for (int j = 0; j < numtot; ++j) {
+
+ float k = kcomb[j];
+
+ /* Shall we correct the position of the k-space bin
+ * to account for the different weights of the modes entering the bin? */
+ if (pow_data->shift_centre_small_k_bins && j < number_of_corrected_bins) {
+ k *= correction_shift_k_values[j];
+ }
+
fprintf(outputfile, "%15.8f %15.8e %15.8e %15.8e\n", s->e->cosmology->z,
- kcomb[j], (pcomb[j] - shot), shot);
+ k, (pcomb[j] - shot), shot);
}
fclose(outputfile);
}
@@ -1273,6 +1315,9 @@ void power_init(struct power_spectrum_data* p, struct swift_params* params,
"WARNING: fold factor is recommended not to exceed 6 for a "
"mass assignment order of 3 (TSC) or below.");
+ p->shift_centre_small_k_bins = parser_get_opt_param_int(
+ params, "PowerSpectrum:shift_centre_small_k_bins", 1);
+
/* Make sensible choices for the k-cuts */
const int kcutn = (p->windoworder >= 3) ? 90 : 70;
const int kcutleft = (int)(p->Ngrid / 256.0 * kcutn);
diff --git a/src/power_spectrum.h b/src/power_spectrum.h
index 97b7feefde6225e590a05f0331b607838a5a24bb..7b2cb73020221b67ae71a5a746fc4d422fc7145d 100644
--- a/src/power_spectrum.h
+++ b/src/power_spectrum.h
@@ -71,6 +71,9 @@ struct power_spectrum_data {
/*! The order of the mass assignment window */
int windoworder;
+ /* Shall we correct the position of the k-space bin? */
+ int shift_centre_small_k_bins;
+
/*! Array of component types to correlate on the "left" side */
enum power_type* types1;