diff --git a/Base/Axis/DiscreteAxis.cpp b/Base/Axis/DiscreteAxis.cpp
index 00725d36b634570558d5e3fcea976a2193b26f96..1917be99a70301e4d05c7e7fbad19f2fecbf1f9e 100644
--- a/Base/Axis/DiscreteAxis.cpp
+++ b/Base/Axis/DiscreteAxis.cpp
@@ -66,15 +66,6 @@ void DiscreteAxis::clip(double lower, double upper)
     m_bins = centers2bins(m_coordinates);
 }
 
-void DiscreteAxis::print(std::ostream& ostr) const
-{
-    ostr << "DiscreteAxis(\"" << axisName() << "\", "
-         << ", [";
-    for (size_t i = 0, fin = m_coordinates.size() - 1; i < fin; ++i)
-        ostr << std::setprecision(17) << m_coordinates[i] << ",";
-    ostr << m_coordinates.back() << "])";
-}
-
 void DiscreteAxis::checkIndex(size_t index) const
 {
     if (m_coordinates.size() > index)
diff --git a/Base/Axis/DiscreteAxis.h b/Base/Axis/DiscreteAxis.h
index 9bb8436c91e8e4682cc4f39bdb4fd33e1915633a..13352637aedee428ced19e0560b30650950bf20c 100644
--- a/Base/Axis/DiscreteAxis.h
+++ b/Base/Axis/DiscreteAxis.h
@@ -29,8 +29,6 @@ public:
     void clip(double lower, double upper) override;
 
 private:
-    void print(std::ostream& ostr) const override;
-
     void checkIndex(size_t index) const;
     void sanityCheck() const;
 
diff --git a/Base/Axis/FixedBinAxis.cpp b/Base/Axis/FixedBinAxis.cpp
index ce5ee8f1393395d8251565209fe1a8c506d09b72..63eae6dd718d6631fe6825b66de5cf9a3bfd608c 100644
--- a/Base/Axis/FixedBinAxis.cpp
+++ b/Base/Axis/FixedBinAxis.cpp
@@ -75,12 +75,6 @@ void FixedBinAxis::clip(double lower, double upper)
     m_bins = bounds2bins(m_nbins, m_start, m_end);
 }
 
-void FixedBinAxis::print(std::ostream& ostr) const
-{
-    ostr << "FixedBinAxis(\"" << axisName() << "\", " << size() << ", " << std::setprecision(17)
-         << min() << ", " << max() << ")";
-}
-
 //... In global namespace.
 
 FixedBinAxis FixedScanAxis(const std::string& name, size_t nbins, double scan_start,
diff --git a/Base/Axis/FixedBinAxis.h b/Base/Axis/FixedBinAxis.h
index f9353522c95a09e7ac22b52ef8e54d2abf1082b1..f28b148c30fd69909aecf59c15be57ea029ac2fd 100644
--- a/Base/Axis/FixedBinAxis.h
+++ b/Base/Axis/FixedBinAxis.h
@@ -34,9 +34,6 @@ public:
 
     void clip(double lower, double upper) override;
 
-protected:
-    void print(std::ostream& ostr) const override;
-
 private:
     size_t m_nbins;
     double m_start;
diff --git a/Base/Axis/IAxis.cpp b/Base/Axis/IAxis.cpp
index 42c28577ace815c14d95848ed8d7bd1c2f71a147..25cf816acfea273b76c97248eee67cb20c5a8096 100644
--- a/Base/Axis/IAxis.cpp
+++ b/Base/Axis/IAxis.cpp
@@ -13,7 +13,10 @@
 //  ************************************************************************************************
 
 #include "Base/Axis/IAxis.h"
+#include "Base/Util/Algorithms.h"
 #include "Base/Util/Assert.h"
+#include <iomanip>
+#include <optional>
 
 IAxis::~IAxis() = default;
 
@@ -98,3 +101,38 @@ void IAxis::clip(const std::pair<double, double> bounds)
 {
     return clip(bounds.first, bounds.second);
 }
+
+std::ostream& operator<<(std::ostream& ostr, const IAxis& ax)
+{
+    ASSERT(ax.size() > 0);
+    std::optional<double> width;
+    for (const Bin1D& b : ax.bins()) {
+        if (!width)
+            width = b.binSize();
+        else if (!BaseUtils::algo::almostEqual(b.binSize(), width.value(), 3 * ax.size())) {
+            width.reset();
+            break;
+        }
+    }
+    if (width) {
+        if (width.value() == 0) {
+            ostr << "DiscreteAxis(\"" << ax.axisName() << "\", "
+                 << ", [" << std::setprecision(17);
+            for (double v : ax.binCenters())
+                ostr << v << ",";
+            ostr << "])";
+        } else {
+            ASSERT(
+                BaseUtils::algo::almostEqual(width.value(), ax.span() / ax.size(), 3 * ax.size()));
+            ostr << "FixedBinAxis(\"" << ax.axisName() << "\", " << ax.size() << ", "
+                 << std::setprecision(17) << ax.min() << ", " << ax.max() << ")";
+        }
+    } else {
+        std::cerr << "DEBUG no uniform width\n" << std::setprecision(17);
+        for (const Bin1D& b : ax.bins())
+            std::cerr << "  b from " << b.lowerBound() << " to " << b.upperBound() << " center "
+                      << b.center() << " width " << b.binSize() << std::endl;
+        ASSERT(false);
+    }
+    return ostr;
+}
diff --git a/Base/Axis/IAxis.h b/Base/Axis/IAxis.h
index 632edea6522b4e185371c7ca44839467be5e6bc9..a1a42f4f3a10370adc89ffa6b583d58e42947ee3 100644
--- a/Base/Axis/IAxis.h
+++ b/Base/Axis/IAxis.h
@@ -64,7 +64,7 @@ public:
 
     const Bin1D& bin(size_t i) const;
     double binCenter(size_t i) const;
-
+    const std::vector<Bin1D>& bins() const { return m_bins; }
     std::vector<double> binCenters() const;
 
     //! find bin index which is best match for given value
@@ -79,15 +79,9 @@ public:
 
     bool operator==(const IAxis& right) const;
 
-    friend std::ostream& operator<<(std::ostream& ostr, const IAxis& m)
-    {
-        m.print(ostr);
-        return ostr;
-    }
+    friend std::ostream& operator<<(std::ostream& ostr, const IAxis& m);
 
 protected:
-    virtual void print(std::ostream& ostr) const = 0;
-
     // private:
     std::string m_name;
     std::vector<Bin1D> m_bins;
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index f12810e1dee278aec0648e2aa2a2c7ee8852b6f9..14e35e0a6cb5c99a6cb7d3367b4b71ab4300c71c 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -1898,6 +1898,10 @@ class IAxis(object):
         r"""binCenter(IAxis self, size_t i) -> double"""
         return _libBornAgainBase.IAxis_binCenter(self, i)
 
+    def bins(self):
+        r"""bins(IAxis self) -> std::vector< Bin1D,std::allocator< Bin1D > > const &"""
+        return _libBornAgainBase.IAxis_bins(self)
+
     def binCenters(self):
         r"""binCenters(IAxis self) -> vdouble1d_t"""
         return _libBornAgainBase.IAxis_binCenters(self)
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index 239022d14f23341105701eac0a2a611d94f57867..200719b045d193478b52a5888ca0790a323e8e64 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -3425,24 +3425,25 @@ namespace Swig {
 #define SWIGTYPE_p_std__lessT_std__string_t swig_types[36]
 #define SWIGTYPE_p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t swig_types[37]
 #define SWIGTYPE_p_std__pairT_double_double_t swig_types[38]
-#define SWIGTYPE_p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t swig_types[39]
-#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[40]
-#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[41]
-#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[42]
-#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[43]
-#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[44]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[45]
-#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[46]
-#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[47]
-#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[48]
-#define SWIGTYPE_p_swig__SwigPyIterator swig_types[49]
-#define SWIGTYPE_p_unsigned_char swig_types[50]
-#define SWIGTYPE_p_unsigned_int swig_types[51]
-#define SWIGTYPE_p_unsigned_long_long swig_types[52]
-#define SWIGTYPE_p_unsigned_short swig_types[53]
-#define SWIGTYPE_p_value_type swig_types[54]
-static swig_type_info *swig_types[56];
-static swig_module_info swig_module = {swig_types, 55, 0, 0, 0, 0};
+#define SWIGTYPE_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t swig_types[39]
+#define SWIGTYPE_p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t swig_types[40]
+#define SWIGTYPE_p_std__vectorT_double_std__allocatorT_double_t_t swig_types[41]
+#define SWIGTYPE_p_std__vectorT_int_std__allocatorT_int_t_t swig_types[42]
+#define SWIGTYPE_p_std__vectorT_std__complexT_double_t_std__allocatorT_std__complexT_double_t_t_t swig_types[43]
+#define SWIGTYPE_p_std__vectorT_std__pairT_double_double_t_std__allocatorT_std__pairT_double_double_t_t_t swig_types[44]
+#define SWIGTYPE_p_std__vectorT_std__string_std__allocatorT_std__string_t_t swig_types[45]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_double_std__allocatorT_double_t_t_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t_t swig_types[46]
+#define SWIGTYPE_p_std__vectorT_std__vectorT_int_std__allocatorT_int_t_t_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t_t swig_types[47]
+#define SWIGTYPE_p_std__vectorT_unsigned_int_std__allocatorT_unsigned_int_t_t swig_types[48]
+#define SWIGTYPE_p_std__vectorT_unsigned_long_std__allocatorT_unsigned_long_t_t swig_types[49]
+#define SWIGTYPE_p_swig__SwigPyIterator swig_types[50]
+#define SWIGTYPE_p_unsigned_char swig_types[51]
+#define SWIGTYPE_p_unsigned_int swig_types[52]
+#define SWIGTYPE_p_unsigned_long_long swig_types[53]
+#define SWIGTYPE_p_unsigned_short swig_types[54]
+#define SWIGTYPE_p_value_type swig_types[55]
+static swig_type_info *swig_types[57];
+static swig_module_info swig_module = {swig_types, 56, 0, 0, 0, 0};
 #define SWIG_TypeQuery(name) SWIG_TypeQueryModule(&swig_module, &swig_module, name)
 #define SWIG_MangledTypeQuery(name) SWIG_MangledTypeQueryModule(&swig_module, &swig_module, name)
 
@@ -25577,6 +25578,29 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_IAxis_bins(PyObject *self, PyObject *args) {
+  PyObject *resultobj = 0;
+  IAxis *arg1 = (IAxis *) 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  PyObject *swig_obj[1] ;
+  std::vector< Bin1D,std::allocator< Bin1D > > *result = 0 ;
+  
+  if (!args) SWIG_fail;
+  swig_obj[0] = args;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_IAxis, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "IAxis_bins" "', argument " "1"" of type '" "IAxis const *""'"); 
+  }
+  arg1 = reinterpret_cast< IAxis * >(argp1);
+  result = (std::vector< Bin1D,std::allocator< Bin1D > > *) &((IAxis const *)arg1)->bins();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t, 0 |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_IAxis_binCenters(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   IAxis *arg1 = (IAxis *) 0 ;
@@ -28727,6 +28751,7 @@ static PyMethodDef SwigMethods[] = {
 	 { "IAxis_center", _wrap_IAxis_center, METH_O, "IAxis_center(IAxis self) -> double"},
 	 { "IAxis_bin", _wrap_IAxis_bin, METH_VARARGS, "IAxis_bin(IAxis self, size_t i) -> Bin1D"},
 	 { "IAxis_binCenter", _wrap_IAxis_binCenter, METH_VARARGS, "IAxis_binCenter(IAxis self, size_t i) -> double"},
+	 { "IAxis_bins", _wrap_IAxis_bins, METH_O, "IAxis_bins(IAxis self) -> std::vector< Bin1D,std::allocator< Bin1D > > const &"},
 	 { "IAxis_binCenters", _wrap_IAxis_binCenters, METH_O, "IAxis_binCenters(IAxis self) -> vdouble1d_t"},
 	 { "IAxis_closestIndex", _wrap_IAxis_closestIndex, METH_VARARGS, "IAxis_closestIndex(IAxis self, double value) -> size_t"},
 	 { "IAxis_clip", _wrap_IAxis_clip, METH_VARARGS, "\n"
@@ -28880,6 +28905,7 @@ static swig_type_info _swigt__p_std__invalid_argument = {"_p_std__invalid_argume
 static swig_type_info _swigt__p_std__lessT_std__string_t = {"_p_std__lessT_std__string_t", "std::less< std::string > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t = {"_p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t", "std::map< std::string,double,std::less< std::string >,std::allocator< std::pair< std::string const,double > > > *|std::map< std::string,double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__pairT_double_double_t = {"_p_std__pairT_double_double_t", "std::pair< double,double > *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t = {"_p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t", "std::vector< Bin1D,std::allocator< Bin1D > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t = {"_p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t", "std::vector< IAxis const *,std::allocator< IAxis const * > > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_double_std__allocatorT_double_t_t = {"_p_std__vectorT_double_std__allocatorT_double_t_t", "std::vector< double,std::allocator< double > > *|std::vector< double > *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_std__vectorT_int_std__allocatorT_int_t_t = {"_p_std__vectorT_int_std__allocatorT_int_t_t", "std::vector< int,std::allocator< int > > *|std::vector< int > *", 0, 0, (void*)0, 0};
@@ -28937,6 +28963,7 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_std__lessT_std__string_t,
   &_swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t,
   &_swigt__p_std__pairT_double_double_t,
+  &_swigt__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t,
   &_swigt__p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t,
   &_swigt__p_std__vectorT_double_std__allocatorT_double_t_t,
   &_swigt__p_std__vectorT_int_std__allocatorT_int_t_t,
@@ -28994,6 +29021,7 @@ static swig_cast_info _swigc__p_std__invalid_argument[] = {  {&_swigt__p_std__in
 static swig_cast_info _swigc__p_std__lessT_std__string_t[] = {  {&_swigt__p_std__lessT_std__string_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t[] = {  {&_swigt__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__pairT_double_double_t[] = {  {&_swigt__p_std__pairT_double_double_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t[] = {  {&_swigt__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t[] = {  {&_swigt__p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_double_std__allocatorT_double_t_t[] = {  {&_swigt__p_std__vectorT_double_std__allocatorT_double_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__vectorT_int_std__allocatorT_int_t_t[] = {  {&_swigt__p_std__vectorT_int_std__allocatorT_int_t_t, 0, 0, 0},{0, 0, 0, 0}};
@@ -29051,6 +29079,7 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_std__lessT_std__string_t,
   _swigc__p_std__mapT_std__string_double_std__lessT_std__string_t_std__allocatorT_std__pairT_std__string_const_double_t_t_t,
   _swigc__p_std__pairT_double_double_t,
+  _swigc__p_std__vectorT_Bin1D_std__allocatorT_Bin1D_t_t,
   _swigc__p_std__vectorT_IAxis_const_p_std__allocatorT_IAxis_const_p_t_t,
   _swigc__p_std__vectorT_double_std__allocatorT_double_t_t,
   _swigc__p_std__vectorT_int_std__allocatorT_int_t_t,