diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index e5c84a69d2d2d5a9529683fdfdcfc0b82e4a2ee3..e6333f1474ecf7c41827e7693b4d5c97952ab567 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -123,12 +123,15 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
                     double *spec_vorticity = new double[fs->nshells*9];
                     fs->cospectrum(fs->cvelocity, fs->cvelocity, spec_velocity);
                     fs->cospectrum(fs->cvorticity, fs->cvorticity, spec_vorticity);
+                    max_estimates[0] = max_trS2_estimate;
+                    max_estimates[1] = max_Q_estimate;
+                    max_estimates[2] = max_R_estimate;
                     fs->compute_gradient_statistics(
                         fs->cvelocity,
                         trS2_Q_R_moments,
                         hist_trS2_Q_R,
                         hist_QR2D,
-                        max_vorticity_estimate*max_vorticity_estimate,
+                        max_estimates,
                         histogram_bins,
                         QR2D_histogram_bins);
                     fs->ift_velocity();
diff --git a/bfps/cpp/fluid_solver.cpp b/bfps/cpp/fluid_solver.cpp
index ae89158e3bf641de179aa88615315e39196c9acd..3a6dca4d84f751a5097d21bace03816b70ef4457 100644
--- a/bfps/cpp/fluid_solver.cpp
+++ b/bfps/cpp/fluid_solver.cpp
@@ -544,7 +544,7 @@ void fluid_solver<R>::compute_gradient_statistics( \
         double *moments, \
         ptrdiff_t *hist, \
         ptrdiff_t *QR2D_hist, \
-        double max_estimate, \
+        double max_estimates[], \
         int nbins, \
         int QR2D_nbins) \
 { \
@@ -573,9 +573,9 @@ void fluid_solver<R>::compute_gradient_statistics( \
     dz_u = ra + 4*this->cd->local_size; \
     double binsize[2]; \
     double tmp_max_estimate[4]; \
-    tmp_max_estimate[0] = max_estimate; \
-    tmp_max_estimate[1] = max_estimate; \
-    tmp_max_estimate[2] = max_estimate*sqrt(max_estimate); \
+    tmp_max_estimate[0] = max_estimates[0]; \
+    tmp_max_estimate[1] = max_estimates[1]; \
+    tmp_max_estimate[2] = max_estimates[2]; \
     tmp_max_estimate[3] = 1.0; /* 4th component is gonna be disregarded anyway... */ \
     binsize[0] = 2*tmp_max_estimate[2] / QR2D_nbins; \
     binsize[1] = 2*tmp_max_estimate[1] / QR2D_nbins; \
diff --git a/bfps/cpp/fluid_solver.hpp b/bfps/cpp/fluid_solver.hpp
index 883a3a9664730c970e8383b5a58aeb5709702507..0d184bcc1475c90d55c3f15e09d732f6abda1afe 100644
--- a/bfps/cpp/fluid_solver.hpp
+++ b/bfps/cpp/fluid_solver.hpp
@@ -87,7 +87,7 @@ class fluid_solver:public fluid_solver_base<rnumber>
                 double *moments,
                 ptrdiff_t *histograms_1D,
                 ptrdiff_t *histogram_QR2D,
-                double trS2_max_estimate = 1.0,
+                double max_estimates[3],
                 int nbins_1D = 256,
                 int nbins_2D = 64);
 
diff --git a/bfps/fluid_base.py b/bfps/fluid_base.py
index 9c6b9132c9b40bf1e4e1a16320a693c6eda52d78..ba1c10a39289824e530e6401c58f159d95908fae 100644
--- a/bfps/fluid_base.py
+++ b/bfps/fluid_base.py
@@ -79,6 +79,9 @@ class fluid_particle_base(bfps.code):
         self.parameters['QR2D_histogram_bins'] = 64
         self.parameters['max_velocity_estimate'] = 1.0
         self.parameters['max_vorticity_estimate'] = 1.0
+        self.parameters['max_trS2_estimate'] = 1.0
+        self.parameters['max_Q_estimate'] = 1.0
+        self.parameters['max_R_estimate'] = 1.0
         self.parameters['dealias_type'] = 1
         self.fluid_includes = '#include "fluid_solver.hpp"\n'
         self.fluid_variables = ''