diff --git a/TurTLE/DNS_statistics.py b/TurTLE/DNS_statistics.py
index 452816244d8fce70f83e02d541e842a7a66bd225..b0f3aeb3e27544c32f98feed64c4af09bd2016ad 100644
--- a/TurTLE/DNS_statistics.py
+++ b/TurTLE/DNS_statistics.py
@@ -301,6 +301,7 @@ def read_histogram(
         base_group = 'statistics',
         pp_file = None,
         iteration_range = None,
+        min_estimate = None,
         max_estimate = None,
         niter_stat = None,
         nbins = None,
@@ -326,6 +327,17 @@ def read_histogram(
     if (hh.shape[-1] == 4):
         hh = hh[..., :3]
         max_estimate /= 3**0.5
+    if type(min_estimate) == type(None):
+        if quantity[:5] == 'log2_':
+            if 'min_' + quantity[5:] + '_estimate' in self.parameters:
+                min_estimate = np.log2(self.parameters['min_' + quantity[5:] + '_estimate'])
+            else:
+                min_estimate = -max_estimate # for compatibility with old data sets
+        else:
+            if 'min_' + quantity + '_estimate' in self.parameters:
+                min_estimate = self.parameters['min_' + quantity + '_estimate']
+            else:
+                min_estimate = -max_estimate # for compatibility with old data sets
     ngrid = self.parameters['nx']*self.parameters['ny']*self.parameters['nz']
     if return_full_history:
         npoints_expected = ngrid
@@ -334,7 +346,7 @@ def read_histogram(
         npoints_expected = ngrid*hh.shape[0]
         hh = np.sum(hh, axis = 0)
         npoints = np.sum(hh, axis = 0, keepdims=True)
-    bb0 = np.linspace(-max_estimate, max_estimate, nbins+1)
+    bb0 = np.linspace(min_estimate, max_estimate, nbins+1)
     bb = (bb0[1:] + bb0[:-1])*0.5
     if (undersample_warning_on):
         if not (npoints == npoints_expected).all():
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 8389561470fd7d5fbdfcbd67c07adaeaa7491a90..5f384f3cb673203e79b0721739f1710357e132e2 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1082,6 +1082,7 @@ void field<rnumber, be, fc>::compute_rspace_stats(
                 const hid_t group,
                 const std::string dset_name,
                 const hsize_t toffset,
+                const std::vector<double> min_estimate,
                 const std::vector<double> max_estimate,
                 const unsigned int maximum_moment)
 {
@@ -1150,14 +1151,26 @@ void field<rnumber, be, fc>::compute_rspace_stats(
         MPI_Bcast(&nvals, 1, MPI_INT, 0, this->comm);
         MPI_Bcast(&nbins, 1, MPI_INT, 0, this->comm);
     }
+    assert(nvals == int(min_estimate.size()));
     assert(nvals == int(max_estimate.size()));
+    for(int i=0; i<nvals; i++)
+        assert(max_estimate[i] > min_estimate[i]);
 
     start_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
     shared_array<double> local_moments_threaded(
             nmoments*nvals,
             [&](double* local_moments){
                 std::fill_n(local_moments, nmoments*nvals, 0);
-                if (nvals == 4) local_moments[3] = max_estimate[3];
+                for (int i=0; i<nvals; i++)
+                {
+                    // see https://en.cppreference.com/w/cpp/types/numeric_limits
+                    local_moments[0*nvals+i] = std::numeric_limits<double>::max();
+                    local_moments[(nmoments-1)*nvals+i] = std::numeric_limits<double>::lowest();
+                }
+                if (nvals == 4) {
+                    local_moments[0*nvals+3] = std::numeric_limits<double>::max();
+                    local_moments[(nmoments-1)*nvals+3] = std::numeric_limits<double>::min();
+                }
             });
 
     shared_array<double> val_tmp_threaded(
@@ -1176,7 +1189,7 @@ void field<rnumber, be, fc>::compute_rspace_stats(
 
     double *binsize = new double[nvals];
     for (int i=0; i<nvals; i++)
-        binsize[i] = 2*max_estimate[i] / nbins;
+        binsize[i] = (max_estimate[i] - min_estimate[i]) / nbins;
 
     {
         {
@@ -1202,27 +1215,19 @@ void field<rnumber, be, fc>::compute_rspace_stats(
             if (nvals == int(4))
             {
                 val_tmp[3] = std::sqrt(val_tmp[3]);
-                if (val_tmp[3] < local_moments[0*nvals+3])
-                    local_moments[0*nvals+3] = val_tmp[3];
-                if (val_tmp[3] > local_moments[9*nvals+3])
-                    local_moments[9*nvals+3] = val_tmp[3];
-                const int bin = int(std::floor(val_tmp[3]*2/binsize[3]));
-                if (bin >= 0 && bin < nbins){
-                    local_hist[bin*nvals+3]++;
-                }
             }
-            for (unsigned int i=0; i<ncomp(fc); i++)
+            for (int i=0; i<nvals; i++)
             {
                 if (val_tmp[i] < local_moments[0*nvals+i])
                     local_moments[0*nvals+i] = val_tmp[i];
                 if (val_tmp[i] > local_moments[(nmoments-1)*nvals+i])
                     local_moments[(nmoments-1)*nvals+i] = val_tmp[i];
-                const int bin = int(std::floor((val_tmp[i] + max_estimate[i]) / binsize[i]));
+                const int bin = int(std::floor((val_tmp[i] - min_estimate[i]) / binsize[i]));
                 if (bin >= 0 && bin < nbins)
                     local_hist[bin*nvals+i]++;
             }
             for (unsigned int n=1; n < maximum_moment; n++){
-                for (auto i=0; i<nvals; i++){
+                for (int i=0; i<nvals; i++){
                     local_moments[n*nvals + i] += (pow_tmp[i] = val_tmp[i]*pow_tmp[i]);
                 }
             }
@@ -1235,16 +1240,10 @@ void field<rnumber, be, fc>::compute_rspace_stats(
         {
         TIMEZONE("FIELD_RLOOP::Merge");
         local_moments_threaded.mergeParallel([&](const int idx, const double& v1, const double& v2) -> double {
-            if(nvals == int(4) && idx == 0*nvals+3){
-                return std::min(v1, v2);
-            }
-            if(nvals == int(4) && idx == 9*nvals+3){
-                return std::max(v1, v2);
-            }
-            if(idx < int(ncomp(fc))){
+            if(idx < int(nvals)){
                 return std::min(v1, v2);
             }
-            if(int(nmoments-1)*nvals <= idx && idx < int(int(nmoments-1)*nvals+ncomp(fc))){
+            if(int(nmoments-1)*nvals <= idx){
                 return std::max(v1, v2);
             }
             return v1 + v2;
diff --git a/cpp/field.hpp b/cpp/field.hpp
index cb43d2ae0a6b6298140406f296ca1ddf95248b15..c136761e780eac73c1eada16e5ecece5c78f4e53 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -128,13 +128,36 @@ class field
                 const hsize_t toffset,
                 const std::vector<double> max_estimate,
                 field<rnumber, be, fc> *tmp_field = NULL);
+    
+        void compute_rspace_stats(
+                const hid_t group,
+                const std::string dset_name,
+                const hsize_t toffset,
+                const std::vector<double> max_estimate,
+                const unsigned int maximum_moment = 9) {
+            std::vector<double> min_estimate;
+            min_estimate.resize(max_estimate.size());
+            for (std::size_t ii = 0; ii < max_estimate.size(); ii++)
+                min_estimate[ii] = -max_estimate[ii];
+            if (max_estimate.size() == 4)
+                min_estimate[3] = 0;
+            this->compute_rspace_stats(
+                    group,
+                    dset_name,
+                    toffset,
+                    min_estimate,
+                    max_estimate,
+                    maximum_moment);
+        }
 
         void compute_rspace_stats(
                 const hid_t group,
                 const std::string dset_name,
                 const hsize_t toffset,
+                const std::vector<double> min_estimate,
                 const std::vector<double> max_estimate,
                 const unsigned int maximum_moment = 9);
+ 
 
 
         double compute_CFL_velocity();
diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 4448294eba55c4ed7ae8869b77fcc2c0e4fbdf4d..2930c39b67ccbdd2a8cac751f3d39464e70ee222 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -31,6 +31,7 @@
 template <typename rnumber>
 int NSVE_field_stats<rnumber>::read_parameters(void)
 {
+    TIMEZONE("NSVE_field_stats::read_parameters");
     this->postprocess::read_parameters();
     hid_t parameter_file = H5Fopen(
             (this->simname + std::string(".h5")).c_str(),
diff --git a/cpp/full_code/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index 9e4731c7447b6c10ea9d545ea3a5996955e48921..c1ca04b3f8bdd1cd718b79ed8be80b29c4078bb1 100644
--- a/cpp/full_code/field_single_to_double.cpp
+++ b/cpp/full_code/field_single_to_double.cpp
@@ -46,12 +46,9 @@ int field_single_to_double<rnumber>::initialize(void)
             (this->simname + std::string(".h5")).c_str(),
             H5F_ACC_RDONLY,
             H5P_DEFAULT);
-    hid_t dset = H5Dopen(parameter_file, "/parameters/niter_out", H5P_DEFAULT);
-    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_out);
-    H5Dclose(dset);
     if (H5Lexists(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT))
     {
-        dset = H5Dopen(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT);
+        hid_t dset = H5Dopen(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT);
         H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->checkpoints_per_file);
         H5Dclose(dset);
     }
diff --git a/cpp/full_code/field_single_to_double.hpp b/cpp/full_code/field_single_to_double.hpp
index 1cb0fdf67059b98c33611d0bb924f83c10a749f9..abfb70f6c90061c7b662b45f27bf8a26300171cd 100644
--- a/cpp/full_code/field_single_to_double.hpp
+++ b/cpp/full_code/field_single_to_double.hpp
@@ -33,7 +33,6 @@ class field_single_to_double: public NSVE_field_stats<rnumber>
 {
     public:
         int checkpoints_per_file;
-        int niter_out;
         kspace<FFTW, SMOOTH> *kk;
 
         field<double, FFTW, THREE> *vec_field_double;
diff --git a/cpp/full_code/get_3D_correlations.hpp b/cpp/full_code/get_3D_correlations.hpp
index 6d92890d9c2caa541461296c34ab01f9255518a4..84187186a612f3bea7e32b80d040cb318e3bc0ad 100644
--- a/cpp/full_code/get_3D_correlations.hpp
+++ b/cpp/full_code/get_3D_correlations.hpp
@@ -34,7 +34,6 @@ class get_3D_correlations: public NSVE_field_stats<rnumber>
     public:
         hid_t stat_file;
         int checkpoints_per_file;
-        int niter_out;
         bool full_snapshot_output;
         kspace<FFTW, SMOOTH> *kk;
         field<rnumber, FFTW, THREE> *vel;
diff --git a/cpp/full_code/get_velocity.hpp b/cpp/full_code/get_velocity.hpp
index 40befc377e0af2492324c04a4fee38a1f24e0e0d..ff453a9a0cc9c50e9c003f902b2cb4ab3dff0970 100644
--- a/cpp/full_code/get_velocity.hpp
+++ b/cpp/full_code/get_velocity.hpp
@@ -33,7 +33,6 @@ class get_velocity: public NSVE_field_stats<rnumber>
 {
     public:
         int checkpoints_per_file;
-        int niter_out;
         kspace<FFTW, SMOOTH> *kk;
         field<rnumber, FFTW, THREE> *vel;
 
diff --git a/cpp/full_code/joint_acc_vel_stats.hpp b/cpp/full_code/joint_acc_vel_stats.hpp
index 2abf116835266efe571136f626582ce4e3b262d2..be5f42746801fbb8c9599f06467564b38f9032d2 100644
--- a/cpp/full_code/joint_acc_vel_stats.hpp
+++ b/cpp/full_code/joint_acc_vel_stats.hpp
@@ -37,7 +37,6 @@ class joint_acc_vel_stats: public NSVE_field_stats<rnumber>
         vorticity_equation<rnumber, FFTW> *ve;
 
         int checkpoints_per_file;
-        int niter_out;
         double max_acceleration_estimate;
         double max_velocity_estimate;
 
diff --git a/cpp/full_code/pressure_stats.hpp b/cpp/full_code/pressure_stats.hpp
index 785013a7431db42511c629c07dfca6460a585600..86ae9457b7e0f9c475604c9f7cfdc859c7b07b36 100644
--- a/cpp/full_code/pressure_stats.hpp
+++ b/cpp/full_code/pressure_stats.hpp
@@ -46,7 +46,6 @@ class pressure_stats: public NSVE_field_stats<rnumber>
         field<rnumber, FFTW, ONE> *scalar;
 
         int checkpoints_per_file;
-        int niter_out;
         double max_pressure_estimate;
 
         pressure_stats(
diff --git a/cpp/full_code/resize.cpp b/cpp/full_code/resize.cpp
index 46b035f85bca1c277916df8da6caf06ca69c65a2..ffb19bdcbf5a4bdeb3d08c365a61567d9dc87369 100644
--- a/cpp/full_code/resize.cpp
+++ b/cpp/full_code/resize.cpp
@@ -40,8 +40,6 @@ int resize<rnumber>::initialize(void)
             H5F_ACC_RDONLY,
             H5P_DEFAULT);
 
-    this->niter_out = hdf5_tools::read_value<int>(
-            parameter_file, "/parameters/niter_out");
     H5Fclose(parameter_file);
     // the following ensures the file is free for rank 0 to open in read/write mode
     MPI_Barrier(this->comm);
diff --git a/cpp/full_code/resize.hpp b/cpp/full_code/resize.hpp
index 521c04618fa34433ee43a8f474320e5afec7606c..afb51cd83c8c74b479b1b819d430365be1d647c2 100644
--- a/cpp/full_code/resize.hpp
+++ b/cpp/full_code/resize.hpp
@@ -38,8 +38,6 @@ class resize: public NSVE_field_stats<rnumber>
         int new_ny;
         int new_nz;
 
-        int niter_out;
-
         field<rnumber, FFTW, THREE> *new_field;
 
         resize(
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
index 38556fb20c0f46e9e62fd8a05e6b4ad614fe9ec6..eb61ca171c972c0160b28c00943fe662e59c94ba 100644
--- a/cpp/full_code/write_rpressure.hpp
+++ b/cpp/full_code/write_rpressure.hpp
@@ -34,8 +34,6 @@ template <typename rnumber>
 class write_rpressure: public NSVE_field_stats<rnumber>
 {
     public:
-        int niter_out;
-
         vorticity_equation<rnumber, FFTW> *ve;
         field<rnumber, FFTW, ONE> *pressure;
         field_binary_IO <rnumber, REAL, ONE> *bin_IO_scalar;