diff --git a/Base/Axis/Frame.cpp b/Base/Axis/Frame.cpp
index 7dfe30e3b0bc7c4b7c47ad5de98256524d1d74c9..13d6a0ea95a3a607e3c1c05d559995d95bdd2ae4 100644
--- a/Base/Axis/Frame.cpp
+++ b/Base/Axis/Frame.cpp
@@ -54,6 +54,7 @@ size_t Frame::size() const
 
 const Scale& Frame::axis(size_t k_axis) const
 {
+    ASSERT(k_axis < rank());
     return *m_axes.at(k_axis);
 }
 const Scale& Frame::xAxis() const
@@ -62,11 +63,13 @@ const Scale& Frame::xAxis() const
 }
 const Scale& Frame::yAxis() const
 {
+    ASSERT(1 < rank());
     return *m_axes.at(1);
 }
 
 size_t Frame::projectedSize(size_t k_axis) const
 {
+    ASSERT(k_axis < rank());
     return m_axes[k_axis]->size();
 }
 
diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index 856112f09ee49b6a1e9dca897e4016ce76eb69c8..2c8b5fcb3332a26021e08a3cbdb9b2d409b91561 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -18,8 +18,6 @@
 #include "Base/Util/Assert.h"
 #include <algorithm>
 
-Datafield::Datafield() = default;
-
 Datafield::Datafield(const Frame* frame, const std::vector<double>& values,
                      const std::vector<double>& errSigmas)
     : m_frame(frame)
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index db4c4fd312d2b5b7f39970eee07655e16ce71484..b47438f97c1415049f3d4178f3c79cff0815cb3e 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -15,23 +15,20 @@
 #ifndef BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 #define BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 
-#include <memory>
-#include <vector>
-
 #ifdef BORNAGAIN_PYTHON
 #include "PyCore/Embed/PyObjectDecl.h"
 #endif
-
+#include <memory>
+#include <vector>
 using std::size_t;
 
 class Scale;
 class Frame;
 
-//! Stores radiation power per bin.
+//! Stores n-dimensional data with n associated scales.
 
 class Datafield {
 public:
-    Datafield(); // required by SimulationResult() required by Swig
     //! Constructor that takes ownership of supplied frame and initializes values and errorbars
     Datafield(const Frame* frame, const std::vector<double>& values = {},
               const std::vector<double>& errSigmas = {});
@@ -130,10 +127,11 @@ public:
     //! Sets content of output data to specific value
     void setAllTo(const double& value);
 
-protected:
+private:
     std::unique_ptr<const Frame> m_frame;
+    std::vector<double> m_values;
+    std::vector<double> m_errSigmas;
 
-private:
     //! Creates projection along X. The projections is made by collecting the data in the range
     //! between [ybinlow, ybinup].
     Datafield* create_xProjection(int ybinlow, int ybinup) const;
@@ -142,9 +140,6 @@ private:
     //! between [xbinlow, xbinup].
     Datafield* create_yProjection(int xbinlow, int xbinup) const;
 
-    std::vector<double> m_values;
-    std::vector<double> m_errSigmas;
-
 #endif // SWIG
 };
 
diff --git a/Device/Histo/SimulationResult.cpp b/Device/Histo/SimulationResult.cpp
index 99d01e45dc243e1ecaffd0c96c106712424847ac..693c49b2be776b960551067865ba7b176c39aae7 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -17,13 +17,9 @@
 #include "Base/Axis/Scale.h"
 #include "Base/Util/Assert.h"
 
-SimulationResult::SimulationResult() = default;
-
-SimulationResult::SimulationResult(const Datafield& data, const Frame* coords)
-    : Datafield(coords, data.flatVector(), data.errorSigmas())
+SimulationResult::SimulationResult(const Datafield& data)
+    : Datafield(data)
 {
-    ASSERT(coords);
-    ASSERT(data.rank() == coords->rank());
 }
 
 SimulationResult::SimulationResult(const SimulationResult& other)
@@ -36,8 +32,6 @@ SimulationResult::SimulationResult(SimulationResult&& other) noexcept = default;
 
 SimulationResult::~SimulationResult() = default;
 
-SimulationResult& SimulationResult::operator=(SimulationResult&& other) noexcept = default;
-
 Datafield SimulationResult::extracted_field() const
 {
     return Datafield(frame().clonedAxes(), flatVector());
@@ -48,34 +42,8 @@ std::pair<double, double> SimulationResult::axisMinMax(size_t i) const
     return {frame().axis(i).min(), frame().axis(i).max()};
 }
 
-std::string SimulationResult::name_of_axis(size_t i) const
-{
-    return frame().axis(i).axisName();
-}
-
-const Frame& SimulationResult::frame() const
-{
-    const auto* coordsys = dynamic_cast<const Frame*>(m_frame.get());
-    ASSERT(coordsys);
-    return *coordsys;
-}
-
-#ifdef BORNAGAIN_PYTHON
-PyObject* SimulationResult::array() const
-{
-    Datafield data(frame().clonedAxes(), flatVector());
-    return data.npArray();
-}
-#endif
-
-std::vector<double> SimulationResult::convertedBinCenters() const
-{
-    return convertedBinCenters(0);
-}
-
 std::vector<double> SimulationResult::convertedBinCenters(size_t i_axis) const
 {
-    ASSERT(i_axis < frame().rank());
     return frame().axis(i_axis).binCenters();
 }
 
diff --git a/Device/Histo/SimulationResult.h b/Device/Histo/SimulationResult.h
index d72a650e12be703e8cf5d6973852823758fd2cfa..733328a31c48a74c3b7192235b8911143f9ffeef 100644
--- a/Device/Histo/SimulationResult.h
+++ b/Device/Histo/SimulationResult.h
@@ -24,28 +24,16 @@ class Frame;
 
 class SimulationResult : public Datafield {
 public:
-    SimulationResult(); // required by Swig
-    SimulationResult(const Datafield& data, const Frame* coords);
+    SimulationResult(const Datafield& data);
 
     SimulationResult(const SimulationResult& other);
     SimulationResult(SimulationResult&& other) noexcept;
     ~SimulationResult() override;
 
-    SimulationResult& operator=(SimulationResult&& other) noexcept;
-
     Datafield extracted_field() const;
 
     std::pair<double, double> axisMinMax(size_t i) const;
 
-    std::string name_of_axis(size_t i) const;
-
-    //! Returns intensity data as Python numpy array
-#ifdef BORNAGAIN_PYTHON
-    PyObject* array() const;
-#endif
-
-    std::vector<double> convertedBinCenters() const;
-
     //! Returns axis coordinates as a numpy array. With no parameters given
     //! Returns coordinates of x-axis in default units.
     std::vector<double> convertedBinCenters(size_t i_axis) const;
@@ -54,10 +42,6 @@ public:
     std::string title();
 
 #ifndef SWIG
-
-    //! Returns underlying unit frame
-    const Frame& frame() const;
-
 private:
     std::string m_title;
 #endif // SWIG
diff --git a/Sim/Fitting/SimDataPair.cpp b/Sim/Fitting/SimDataPair.cpp
index a7f941b6f01e4f9c87c9507c102591efd6897bcd..781f82a10f8b81f7c00f554f179583b850d37d76 100644
--- a/Sim/Fitting/SimDataPair.cpp
+++ b/Sim/Fitting/SimDataPair.cpp
@@ -69,7 +69,7 @@ SimulationResult convertData(const ScatteringSimulation& simulation, const Dataf
         throw std::runtime_error(
             "FitObject::init_dataset: Detector and experimental data have different shape");
 
-    return SimulationResult(*roi_data, coordSystem);
+    return SimulationResult(*roi_data);
 }
 
 } // namespace
@@ -122,9 +122,8 @@ void SimDataPair::execSimulation(const mumufit::Parameters& params)
         m_user_weights =
             std::make_unique<SimulationResult>(convertData(*sim2d, *m_raw_user_weights));
     } else {
-        const Frame& frame = m_sim_data->frame();
-        m_exp_data = std::make_unique<SimulationResult>(*m_raw_data, frame.clone());
-        m_user_weights = std::make_unique<SimulationResult>(*m_raw_user_weights, frame.clone());
+        m_exp_data = std::make_unique<SimulationResult>(*m_raw_data);
+        m_user_weights = std::make_unique<SimulationResult>(*m_raw_user_weights);
     }
 
     if (sim2d && containsUncertainties())
@@ -133,7 +132,7 @@ void SimDataPair::execSimulation(const mumufit::Parameters& params)
     else {
         const Frame& frame = m_sim_data->frame();
         auto dummy_array = std::make_unique<Datafield>(frame.clonedAxes());
-        m_uncertainties = std::make_unique<SimulationResult>(*dummy_array, frame.clone());
+        m_uncertainties = std::make_unique<SimulationResult>(*dummy_array);
     }
 }
 
@@ -205,9 +204,7 @@ SimulationResult SimDataPair::relativeDifference() const
     for (size_t i = 0; i < N; ++i)
         data[i] = Numeric::relativeDifference((*m_sim_data)[i], (*m_exp_data)[i]);
 
-    const Frame* f = m_sim_data->frame().clone();
-    Datafield df(f, data);
-    return {df, f};
+    return {Datafield(m_sim_data->frame().clone(), data)};
 }
 
 SimulationResult SimDataPair::absoluteDifference() const
@@ -222,9 +219,7 @@ SimulationResult SimDataPair::absoluteDifference() const
     for (size_t i = 0; i < N; ++i)
         data[i] = std::abs((*m_sim_data)[i] - (*m_exp_data)[i]);
 
-    const Frame* f = m_sim_data->frame().clone();
-    Datafield df(f, data);
-    return {df, f};
+    return {Datafield(m_sim_data->frame().clone(), data)};
 }
 
 void SimDataPair::validate() const
diff --git a/Sim/Simulation/DepthprobeSimulation.cpp b/Sim/Simulation/DepthprobeSimulation.cpp
index c392812e34696a9ed36262c2d6a34f1e31ee9943..3cc8bcc5479cf22ec2df09565143427ff35c1471 100644
--- a/Sim/Simulation/DepthprobeSimulation.cpp
+++ b/Sim/Simulation/DepthprobeSimulation.cpp
@@ -178,5 +178,5 @@ SimulationResult DepthprobeSimulation::packResult()
     std::vector<const Scale*> axes{m_scan->coordinateAxis()->clone(), m_z_axis->clone()};
     auto data = std::make_unique<Datafield>(std::move(axes), m_cache);
 
-    return {*data, simCoordSystem()};
+    return {*data};
 }
diff --git a/Sim/Simulation/OffspecSimulation.cpp b/Sim/Simulation/OffspecSimulation.cpp
index db5e5e2830edd93e6197ca6f3a23958ee923df11..781d25c2e5f93baaa12d855e0858c47721110b34 100644
--- a/Sim/Simulation/OffspecSimulation.cpp
+++ b/Sim/Simulation/OffspecSimulation.cpp
@@ -149,5 +149,5 @@ SimulationResult OffspecSimulation::packResult()
             intensity_map[j * ny + i % ny] += detector_image[i];
     }
 
-    return {intensity_map, simCoordSystem()};
+    return {intensity_map};
 }
diff --git a/Sim/Simulation/ScatteringSimulation.cpp b/Sim/Simulation/ScatteringSimulation.cpp
index 200ea1fd1b5177fe8f804e26a1bdc048ec72a816..28c5d931c32c78cabceecb702edec2aa7d575dd1 100644
--- a/Sim/Simulation/ScatteringSimulation.cpp
+++ b/Sim/Simulation/ScatteringSimulation.cpp
@@ -129,5 +129,5 @@ SimulationResult ScatteringSimulation::packResult()
         [&](const auto it) { detectorMap[it.roiIndex()] = m_cache[elementIndex++]; });
     m_detector->applyDetectorResolution(&detectorMap);
 
-    return {detectorMap, simCoordSystem()};
+    return {detectorMap};
 }
diff --git a/Sim/Simulation/SpecularSimulation.cpp b/Sim/Simulation/SpecularSimulation.cpp
index 253803885efe8aec30718d3323cb05538971a19d..dbab73e1ad520fc70f0f18336758f0d06cd4476d 100644
--- a/Sim/Simulation/SpecularSimulation.cpp
+++ b/Sim/Simulation/SpecularSimulation.cpp
@@ -106,5 +106,5 @@ SimulationResult SpecularSimulation::packResult()
     }
 
     Datafield data({m_scan->coordinateAxis()->clone()}, vec);
-    return {data, simCoordSystem()};
+    return {data};
 }
diff --git a/Tests/Py/Functional/PolarizedNoAnalyzer.py b/Tests/Py/Functional/PolarizedNoAnalyzer.py
index 535407dad12e6a61c6f0e773a980c8d3619d3775..ac2e0a718bda2a1b6021b652b8be73f0c32f1f35 100755
--- a/Tests/Py/Functional/PolarizedNoAnalyzer.py
+++ b/Tests/Py/Functional/PolarizedNoAnalyzer.py
@@ -46,7 +46,7 @@ def run_simulation(polarizer_dir, analyzer_dir=None):
 
     result = simulation.simulate()
 
-    return result.convertedBinCenters(), result.array()
+    return result.convertedBinCenters(0), result.npArray()
 
 
 if __name__ == '__main__':
diff --git a/Wrap/Python/ba_fitmonitor.py b/Wrap/Python/ba_fitmonitor.py
index 7d1f9b0149d9972b603d4befbefbb7c2545ce619..ee101cfbce6dcb8b1d8d7e41e6e564bdf4423d5c 100644
--- a/Wrap/Python/ba_fitmonitor.py
+++ b/Wrap/Python/ba_fitmonitor.py
@@ -80,7 +80,7 @@ class PlotterGISAS(Plotter):
         self.make_subplot(1)
 
         # same limits for both plots
-        arr = real_data.array()
+        arr = real_data.npArray()
         zmax = np.amax(arr) if self._zmax is None else self._zmax
         zmin = zmax*1e-6 if self._zmin is None else self._zmin
 
@@ -155,26 +155,26 @@ class PlotterSpecular:
         unc_data = fit_objective.uncertaintyData()
 
         # data values
-        sim_values = sim_data.array()
-        exp_values = exp_data.array()
-        unc_values = None if unc_data is None else unc_data.array()
+        sim_values = sim_data.npArray()
+        exp_values = exp_data.npArray()
+        unc_values = None if unc_data is None else unc_data.npArray()
 
         # default font properties dictionary to use
         font = { 'size': 16 }
 
         plt.yscale('log')
         plt.ylim((0.5*np.min(exp_values), 5*np.max(exp_values)))
-        plt.plot(exp_data.convertedBinCenters(), exp_values, 'k--')
+        plt.plot(exp_data.convertedBinCenters(0), exp_values, 'k--')
         if unc_values is not None:
-            plt.plot(exp_data.convertedBinCenters(),
+            plt.plot(exp_data.convertedBinCenters(0),
                      exp_values - unc_values,
                      'xkcd:grey',
                      alpha=0.6)
-            plt.plot(exp_data.convertedBinCenters(),
+            plt.plot(exp_data.convertedBinCenters(0),
                      exp_values + unc_values,
                      'xkcd:grey',
                      alpha=0.6)
-        plt.plot(sim_data.convertedBinCenters(), sim_values, 'b')
+        plt.plot(sim_data.convertedBinCenters(0), sim_values, 'b')
 
         xlabel = bp.get_axes_labels(exp_data)[0]
         legend = ['Experiment', 'BornAgain']
diff --git a/Wrap/Python/ba_plot.py b/Wrap/Python/ba_plot.py
index 1e456d7dd6876ab86756dc234b66cb8543d49d76..27c5756979b051c44077938b72c178aa57cbe865 100644
--- a/Wrap/Python/ba_plot.py
+++ b/Wrap/Python/ba_plot.py
@@ -164,7 +164,7 @@ def get_axes_labels(result):
     """
     labels = []
     for i in range(result.rank()):
-        labels.append(translate_axis_label(result.name_of_axis(i)))
+        labels.append(translate_axis_label(result.axis(i).axisName()))
 
     return labels
 
@@ -195,8 +195,8 @@ def plot_specular_curve(result):
     :param result: SimulationResult from SpecularSimulation
     Used internally.
     """
-    intensity = result.array()
-    x_axis = result.convertedBinCenters()
+    intensity = result.npArray()
+    x_axis = result.convertedBinCenters(0)
 
     xlabel = plotargs.pop('xlabel', get_axes_labels(result)[0])
     ylabel = plotargs.pop('ylabel', "Intensity")
@@ -409,7 +409,7 @@ def plot_simres(result, **kwargs):
     if 'ylabel' not in kwargs:
         kwargs['ylabel'] = axes_labels[1]
 
-    return plot_array(result.array(), axes_limits=axes_limits, **kwargs)
+    return plot_array(result.npArray(), axes_limits=axes_limits, **kwargs)
 
 #  **************************************************************************  #
 #  standard user calls
@@ -453,7 +453,7 @@ def plot_simulation_result(result):
     if datfile:
         check_or_save(result, datfile)
 
-    if len(result.array().shape) == 1:
+    if len(result.npArray().shape) == 1:
         # 1D data => assume specular simulation
         plot_specular_curve(result, **plotargs)
     else:
@@ -489,7 +489,7 @@ def make_plot(results, ncol):
         ax = multiPlot.axes[i]
         axes_limits = get_axes_limits(result)
 
-        im = ax.imshow(result.array(),
+        im = ax.imshow(result.npArray(),
                        cmap=cmap,
                        norm=norm,
                        extent=axes_limits,
@@ -521,8 +521,8 @@ def plot_multicurve(results, xlabel, ylabel):
 
     legend = []
     for result in results:
-        x = result.convertedBinCenters()
-        y = result.array()
+        x = result.convertedBinCenters(0)
+        y = result.npArray()
         legend.append(result.title())
         plt.plot(x, y)
 
diff --git a/Wrap/Swig/libBornAgainSim.i b/Wrap/Swig/libBornAgainSim.i
index 68def7b4104cc8e37c16accc468c25216f504837..18bf3b538e4a6d677d96aff0808cf1ec907b5688 100644
--- a/Wrap/Swig/libBornAgainSim.i
+++ b/Wrap/Swig/libBornAgainSim.i
@@ -54,6 +54,7 @@
 #include "Base/Axis/Scale.h"
 #include "Base/Axis/Frame.h"
 
+#include "Device/Data/Datafield.h"
 #include "Device/Histo/SimulationResult.h"
 
 #include "Fit/Minimizer/MinimizerResult.h"
@@ -86,6 +87,8 @@
 %include "fromParam.i"
 
 %import(module="libBornAgainSample") "Sample/Scattering/ISampleNode.h"
+%import(module="libBornAgainDevice") "Device/Data/Datafield.h"
+%import(module="libBornAgainDevice") "Device/Histo/SimulationResult.h"
 
 %template(swig_dummy_type_const_inode_vector) std::vector<const INode*>;
 
diff --git a/auto/Examples/fit/scatter2d/fit2d.py b/auto/Examples/fit/scatter2d/fit2d.py
index e8c237c22db2c317b97bdf15591c049838cb609c..624f8ce0e2f392e669c026b408123bd881d01afa 100755
--- a/auto/Examples/fit/scatter2d/fit2d.py
+++ b/auto/Examples/fit/scatter2d/fit2d.py
@@ -34,7 +34,7 @@ def real_data():
     """
     simulation = get_simulation({'radius': 5 * nm})
     result = simulation.simulate()
-    return result.array()
+    return result.npArray()
 
 
 if __name__ == '__main__':
diff --git a/auto/Examples/fit/scatter2d/model2_hexlattice.py b/auto/Examples/fit/scatter2d/model2_hexlattice.py
index 7cd0469f216afcc6699780cc6fb0fe97dfc1a164..a08df601e1bae78be55847c7f7e26d3c652812f6 100755
--- a/auto/Examples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/Examples/fit/scatter2d/model2_hexlattice.py
@@ -57,7 +57,7 @@ def create_real_data():
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
-    real_data = result.array()
+    real_data = result.npArray()
 
     # spoiling simulated data with noise to produce "real" data
     np.random.seed(0)
diff --git a/auto/Examples/fit/specular/Honeycomb_fit.py b/auto/Examples/fit/specular/Honeycomb_fit.py
index e9dc5fa803305e2f3892401c3e4d4c717664c62a..f8f9a96073799af461e0c7037003be3cffb0da1a 100755
--- a/auto/Examples/fit/specular/Honeycomb_fit.py
+++ b/auto/Examples/fit/specular/Honeycomb_fit.py
@@ -109,7 +109,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = numpy.array(result.convertedBinCenters())
+    q = numpy.array(result.convertedBinCenters(0))
     r = numpy.array(result.array())
 
     return q, r
diff --git a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..01e333315874c14b27efde71b556f2ec51a28da3 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
@@ -101,8 +101,8 @@ def qr(result):
     Returns two arrays that hold the q-values as well as the
     reflectivity from a given simulation result
     """
-    q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    q = numpy.array(result.convertedBinCenters(0))
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/Examples/fit/specular/Pt_layer_fit.py b/auto/Examples/fit/specular/Pt_layer_fit.py
index de39eb13832791629717cc2040e15e6a41b1357b..170527753bf1d6f747602d8bd058eb69b21bbe20 100755
--- a/auto/Examples/fit/specular/Pt_layer_fit.py
+++ b/auto/Examples/fit/specular/Pt_layer_fit.py
@@ -88,7 +88,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = np.array(result.convertedBinCenters())
+    q = np.array(result.convertedBinCenters(0))
     r = np.array(result.array())
 
     return q, r
diff --git a/auto/Examples/specular/PolarizedSpinAsymmetry.py b/auto/Examples/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..01e333315874c14b27efde71b556f2ec51a28da3 100755
--- a/auto/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/specular/PolarizedSpinAsymmetry.py
@@ -101,8 +101,8 @@ def qr(result):
     Returns two arrays that hold the q-values as well as the
     reflectivity from a given simulation result
     """
-    q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    q = numpy.array(result.convertedBinCenters(0))
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/MiniExamples/fit/scatter2d/fit2d.py b/auto/MiniExamples/fit/scatter2d/fit2d.py
index 7ef1a5f8e928e78fc0095eac91ab19ef25726eae..edf243eee571ee3903426fb28a663e25ee2e4fd4 100755
--- a/auto/MiniExamples/fit/scatter2d/fit2d.py
+++ b/auto/MiniExamples/fit/scatter2d/fit2d.py
@@ -34,7 +34,7 @@ def real_data():
     """
     simulation = get_simulation({'radius': 5 * nm})
     result = simulation.simulate()
-    return result.array()
+    return result.npArray()
 
 
 if __name__ == '__main__':
diff --git a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
index c6d690e3e868b3331cb67419bca8fb62d3b020c5..cad7b9101b15f28e95558ce14f0238f4f69f03e0 100755
--- a/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
+++ b/auto/MiniExamples/fit/scatter2d/model2_hexlattice.py
@@ -57,7 +57,7 @@ def create_real_data():
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
-    real_data = result.array()
+    real_data = result.npArray()
 
     # spoiling simulated data with noise to produce "real" data
     np.random.seed(0)
diff --git a/auto/MiniExamples/fit/specular/Honeycomb_fit.py b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
index e9dc5fa803305e2f3892401c3e4d4c717664c62a..f8f9a96073799af461e0c7037003be3cffb0da1a 100755
--- a/auto/MiniExamples/fit/specular/Honeycomb_fit.py
+++ b/auto/MiniExamples/fit/specular/Honeycomb_fit.py
@@ -109,7 +109,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = numpy.array(result.convertedBinCenters())
+    q = numpy.array(result.convertedBinCenters(0))
     r = numpy.array(result.array())
 
     return q, r
diff --git a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..01e333315874c14b27efde71b556f2ec51a28da3 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
@@ -101,8 +101,8 @@ def qr(result):
     Returns two arrays that hold the q-values as well as the
     reflectivity from a given simulation result
     """
-    q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    q = numpy.array(result.convertedBinCenters(0))
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/MiniExamples/fit/specular/Pt_layer_fit.py b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
index de39eb13832791629717cc2040e15e6a41b1357b..170527753bf1d6f747602d8bd058eb69b21bbe20 100755
--- a/auto/MiniExamples/fit/specular/Pt_layer_fit.py
+++ b/auto/MiniExamples/fit/specular/Pt_layer_fit.py
@@ -88,7 +88,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = np.array(result.convertedBinCenters())
+    q = np.array(result.convertedBinCenters(0))
     r = np.array(result.array())
 
     return q, r
diff --git a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..01e333315874c14b27efde71b556f2ec51a28da3 100755
--- a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
@@ -101,8 +101,8 @@ def qr(result):
     Returns two arrays that hold the q-values as well as the
     reflectivity from a given simulation result
     """
-    q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    q = numpy.array(result.convertedBinCenters(0))
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index d4a72c5a3c83a2d6f38877b81e3260f5afdf10e0..e6a3a0c7d0043602739e62f60019258109dd3027 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2041,7 +2041,6 @@ class Datafield(object):
 
     def __init__(self, *args):
         r"""
-        __init__(Datafield self) -> Datafield
         __init__(Datafield self, Frame frame, vdouble1d_t values={}, vdouble1d_t errSigmas={}) -> Datafield
         __init__(Datafield self, std::vector< Scale const *,std::allocator< Scale const * > > && axes, vdouble1d_t values={}, vdouble1d_t errSigmas={}) -> Datafield
         __init__(Datafield self, Datafield arg2) -> Datafield
@@ -3043,8 +3042,7 @@ class SimulationResult(Datafield):
 
     def __init__(self, *args):
         r"""
-        __init__(SimulationResult self) -> SimulationResult
-        __init__(SimulationResult self, Datafield data, Frame coords) -> SimulationResult
+        __init__(SimulationResult self, Datafield data) -> SimulationResult
         __init__(SimulationResult self, SimulationResult other) -> SimulationResult
         __init__(SimulationResult self, SimulationResult other) -> SimulationResult
         """
@@ -3059,20 +3057,9 @@ class SimulationResult(Datafield):
         r"""axisMinMax(SimulationResult self, size_t i) -> pvacuum_double_t"""
         return _libBornAgainDevice.SimulationResult_axisMinMax(self, i)
 
-    def name_of_axis(self, i):
-        r"""name_of_axis(SimulationResult self, size_t i) -> std::string"""
-        return _libBornAgainDevice.SimulationResult_name_of_axis(self, i)
-
-    def array(self):
-        r"""array(SimulationResult self) -> PyObject *"""
-        return _libBornAgainDevice.SimulationResult_array(self)
-
-    def convertedBinCenters(self, *args):
-        r"""
-        convertedBinCenters(SimulationResult self) -> vdouble1d_t
-        convertedBinCenters(SimulationResult self, size_t i_axis) -> vdouble1d_t
-        """
-        return _libBornAgainDevice.SimulationResult_convertedBinCenters(self, *args)
+    def convertedBinCenters(self, i_axis):
+        r"""convertedBinCenters(SimulationResult self, size_t i_axis) -> vdouble1d_t"""
+        return _libBornAgainDevice.SimulationResult_convertedBinCenters(self, i_axis)
 
     def setTitle(self, title):
         r"""setTitle(SimulationResult self, std::string const & title)"""
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 8843dd031ffd967d1399cf36a8418fb86d8d33fe..e7bf89bfdcaec0e9e7aa3baf7285ebfcd72a8aef 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -27360,20 +27360,7 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Datafield *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Datafield *)new Datafield();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Datafield, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Frame *arg1 = (Frame *) 0 ;
   std::vector< double,std::allocator< double > > *arg2 = 0 ;
@@ -27424,7 +27411,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Frame *arg1 = (Frame *) 0 ;
   std::vector< double,std::allocator< double > > *arg2 = 0 ;
@@ -27460,7 +27447,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_3(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Frame *arg1 = (Frame *) 0 ;
   void *argp1 = 0 ;
@@ -27481,7 +27468,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_4(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_3(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::vector< Scale const *,std::allocator< Scale const * > > *arg1 = 0 ;
   std::vector< double,std::allocator< double > > *arg2 = 0 ;
@@ -27541,7 +27528,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_5(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_4(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::vector< Scale const *,std::allocator< Scale const * > > *arg1 = 0 ;
   std::vector< double,std::allocator< double > > *arg2 = 0 ;
@@ -27586,7 +27573,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_6(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_5(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::vector< Scale const *,std::allocator< Scale const * > > *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -27616,7 +27603,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_7(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_6(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Datafield *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -27640,7 +27627,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_8(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_7(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Datafield *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -27678,16 +27665,13 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
   
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_Datafield", 0, 3, argv))) SWIG_fail;
   --argc;
-  if (argc == 0) {
-    return _wrap_new_Datafield__SWIG_0(self, argc, argv);
-  }
   if (argc == 1) {
     int _v = 0;
     void *vptr = 0;
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_Frame, 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_Datafield__SWIG_3(self, argc, argv);
+      return _wrap_new_Datafield__SWIG_2(self, argc, argv);
     }
   }
   if (argc == 1) {
@@ -27696,7 +27680,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_std__vectorT_Scale_const_p_std__allocatorT_Scale_const_p_t_t, SWIG_POINTER_NO_NULL);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_Datafield__SWIG_6(self, argc, argv);
+      return _wrap_new_Datafield__SWIG_5(self, argc, argv);
     }
   }
   if (argc == 1) {
@@ -27704,7 +27688,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_Datafield__SWIG_7(self, argc, argv);
+      return _wrap_new_Datafield__SWIG_6(self, argc, argv);
     }
   }
   if (argc == 1) {
@@ -27713,7 +27697,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_Datafield__SWIG_8(self, argc, argv);
+      return _wrap_new_Datafield__SWIG_7(self, argc, argv);
     }
   }
   if (argc == 2) {
@@ -27725,7 +27709,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
       int res = swig::asptr(argv[1], (std::vector< double,std::allocator< double > >**)(0));
       _v = SWIG_CheckState(res);
       if (_v) {
-        return _wrap_new_Datafield__SWIG_5(self, argc, argv);
+        return _wrap_new_Datafield__SWIG_4(self, argc, argv);
       }
     }
   }
@@ -27738,7 +27722,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
       int res = swig::asptr(argv[1], (std::vector< double,std::allocator< double > >**)(0));
       _v = SWIG_CheckState(res);
       if (_v) {
-        return _wrap_new_Datafield__SWIG_2(self, argc, argv);
+        return _wrap_new_Datafield__SWIG_1(self, argc, argv);
       }
     }
   }
@@ -27754,7 +27738,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
         int res = swig::asptr(argv[2], (std::vector< double,std::allocator< double > >**)(0));
         _v = SWIG_CheckState(res);
         if (_v) {
-          return _wrap_new_Datafield__SWIG_1(self, argc, argv);
+          return _wrap_new_Datafield__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -27771,7 +27755,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
         int res = swig::asptr(argv[2], (std::vector< double,std::allocator< double > >**)(0));
         _v = SWIG_CheckState(res);
         if (_v) {
-          return _wrap_new_Datafield__SWIG_4(self, argc, argv);
+          return _wrap_new_Datafield__SWIG_3(self, argc, argv);
         }
       }
     }
@@ -27780,7 +27764,6 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_Datafield'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Datafield::Datafield()\n"
     "    Datafield::Datafield(Frame const *,std::vector< double,std::allocator< double > > const &,std::vector< double,std::allocator< double > > const &)\n"
     "    Datafield::Datafield(Frame const *,std::vector< double,std::allocator< double > > const &)\n"
     "    Datafield::Datafield(Frame const *)\n"
@@ -36947,30 +36930,14 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  SimulationResult *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (SimulationResult *)new SimulationResult();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_SimulationResult, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   Datafield *arg1 = 0 ;
-  Frame *arg2 = (Frame *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  void *argp2 = 0 ;
-  int res2 = 0 ;
   SimulationResult *result = 0 ;
   
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1, SWIGTYPE_p_Datafield,  0  | 0);
   if (!SWIG_IsOK(res1)) {
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "new_SimulationResult" "', argument " "1"" of type '" "Datafield const &""'"); 
@@ -36979,12 +36946,7 @@ SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_1(PyObject *self, Py_ssize
     SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_SimulationResult" "', argument " "1"" of type '" "Datafield const &""'"); 
   }
   arg1 = reinterpret_cast< Datafield * >(argp1);
-  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2,SWIGTYPE_p_Frame, 0 |  0 );
-  if (!SWIG_IsOK(res2)) {
-    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "new_SimulationResult" "', argument " "2"" of type '" "Frame const *""'"); 
-  }
-  arg2 = reinterpret_cast< Frame * >(argp2);
-  result = (SimulationResult *)new SimulationResult((Datafield const &)*arg1,(Frame const *)arg2);
+  result = (SimulationResult *)new SimulationResult((Datafield const &)*arg1);
   resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_SimulationResult, SWIG_POINTER_NEW |  0 );
   return resultobj;
 fail:
@@ -36992,7 +36954,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -37016,7 +36978,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_3(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_SimulationResult__SWIG_2(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = 0 ;
   void *argp1 = 0 ;
@@ -37048,21 +37010,18 @@ fail:
 
 SWIGINTERN PyObject *_wrap_new_SimulationResult(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
-  PyObject *argv[3] = {
+  PyObject *argv[2] = {
     0
   };
   
-  if (!(argc = SWIG_Python_UnpackTuple(args, "new_SimulationResult", 0, 2, argv))) SWIG_fail;
+  if (!(argc = SWIG_Python_UnpackTuple(args, "new_SimulationResult", 0, 1, argv))) SWIG_fail;
   --argc;
-  if (argc == 0) {
-    return _wrap_new_SimulationResult__SWIG_0(self, argc, argv);
-  }
   if (argc == 1) {
     int _v = 0;
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_SimulationResult, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_SimulationResult__SWIG_2(self, argc, argv);
+      return _wrap_new_SimulationResult__SWIG_1(self, argc, argv);
     }
   }
   if (argc == 1) {
@@ -37071,28 +37030,22 @@ SWIGINTERN PyObject *_wrap_new_SimulationResult(PyObject *self, PyObject *args)
     int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_SimulationResult, SWIG_POINTER_NO_NULL);
     _v = SWIG_CheckState(res);
     if (_v) {
-      return _wrap_new_SimulationResult__SWIG_3(self, argc, argv);
+      return _wrap_new_SimulationResult__SWIG_2(self, argc, argv);
     }
   }
-  if (argc == 2) {
+  if (argc == 1) {
     int _v = 0;
     int res = SWIG_ConvertPtr(argv[0], 0, SWIGTYPE_p_Datafield, SWIG_POINTER_NO_NULL | 0);
     _v = SWIG_CheckState(res);
     if (_v) {
-      void *vptr = 0;
-      int res = SWIG_ConvertPtr(argv[1], &vptr, SWIGTYPE_p_Frame, 0);
-      _v = SWIG_CheckState(res);
-      if (_v) {
-        return _wrap_new_SimulationResult__SWIG_1(self, argc, argv);
-      }
+      return _wrap_new_SimulationResult__SWIG_0(self, argc, argv);
     }
   }
   
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_SimulationResult'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    SimulationResult::SimulationResult()\n"
-    "    SimulationResult::SimulationResult(Datafield const &,Frame const *)\n"
+    "    SimulationResult::SimulationResult(Datafield const &)\n"
     "    SimulationResult::SimulationResult(SimulationResult const &)\n"
     "    SimulationResult::SimulationResult(SimulationResult &&)\n");
   return 0;
@@ -37127,7 +37080,7 @@ SWIGINTERN PyObject *_wrap_SimulationResult_extracted_field(PyObject *self, PyOb
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
-  Datafield result;
+  SwigValueWrapper< Datafield > result;
   
   if (!args) SWIG_fail;
   swig_obj[0] = args;
@@ -37174,7 +37127,7 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_SimulationResult_name_of_axis(PyObject *self, PyObject *args) {
+SWIGINTERN PyObject *_wrap_SimulationResult_convertedBinCenters(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = (SimulationResult *) 0 ;
   size_t arg2 ;
@@ -37183,82 +37136,9 @@ SWIGINTERN PyObject *_wrap_SimulationResult_name_of_axis(PyObject *self, PyObjec
   size_t val2 ;
   int ecode2 = 0 ;
   PyObject *swig_obj[2] ;
-  std::string result;
-  
-  if (!SWIG_Python_UnpackTuple(args, "SimulationResult_name_of_axis", 2, 2, swig_obj)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_name_of_axis" "', argument " "1"" of type '" "SimulationResult const *""'"); 
-  }
-  arg1 = reinterpret_cast< SimulationResult * >(argp1);
-  ecode2 = SWIG_AsVal_size_t(swig_obj[1], &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "SimulationResult_name_of_axis" "', argument " "2"" of type '" "size_t""'");
-  } 
-  arg2 = static_cast< size_t >(val2);
-  result = ((SimulationResult const *)arg1)->name_of_axis(arg2);
-  resultobj = SWIG_From_std_string(static_cast< std::string >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_SimulationResult_array(PyObject *self, PyObject *args) {
-  PyObject *resultobj = 0;
-  SimulationResult *arg1 = (SimulationResult *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  PyObject *swig_obj[1] ;
-  PyObject *result = 0 ;
-  
-  if (!args) SWIG_fail;
-  swig_obj[0] = args;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_array" "', argument " "1"" of type '" "SimulationResult const *""'"); 
-  }
-  arg1 = reinterpret_cast< SimulationResult * >(argp1);
-  result = (PyObject *)((SimulationResult const *)arg1)->array();
-  resultobj = result;
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_SimulationResult_convertedBinCenters__SWIG_0(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  SimulationResult *arg1 = (SimulationResult *) 0 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
   std::vector< double,std::allocator< double > > result;
   
-  if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
-  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
-  if (!SWIG_IsOK(res1)) {
-    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_convertedBinCenters" "', argument " "1"" of type '" "SimulationResult const *""'"); 
-  }
-  arg1 = reinterpret_cast< SimulationResult * >(argp1);
-  result = ((SimulationResult const *)arg1)->convertedBinCenters();
-  resultobj = swig::from(static_cast< std::vector< double,std::allocator< double > > >(result));
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_SimulationResult_convertedBinCenters__SWIG_1(PyObject *self, Py_ssize_t nobjs, PyObject **swig_obj) {
-  PyObject *resultobj = 0;
-  SimulationResult *arg1 = (SimulationResult *) 0 ;
-  size_t arg2 ;
-  void *argp1 = 0 ;
-  int res1 = 0 ;
-  size_t val2 ;
-  int ecode2 = 0 ;
-  std::vector< double,std::allocator< double > > result;
-  
-  if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
+  if (!SWIG_Python_UnpackTuple(args, "SimulationResult_convertedBinCenters", 2, 2, swig_obj)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_SimulationResult, 0 |  0 );
   if (!SWIG_IsOK(res1)) {
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "SimulationResult_convertedBinCenters" "', argument " "1"" of type '" "SimulationResult const *""'"); 
@@ -37277,48 +37157,6 @@ fail:
 }
 
 
-SWIGINTERN PyObject *_wrap_SimulationResult_convertedBinCenters(PyObject *self, PyObject *args) {
-  Py_ssize_t argc;
-  PyObject *argv[3] = {
-    0
-  };
-  
-  if (!(argc = SWIG_Python_UnpackTuple(args, "SimulationResult_convertedBinCenters", 0, 2, argv))) SWIG_fail;
-  --argc;
-  if (argc == 1) {
-    int _v = 0;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_SimulationResult, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      return _wrap_SimulationResult_convertedBinCenters__SWIG_0(self, argc, argv);
-    }
-  }
-  if (argc == 2) {
-    int _v = 0;
-    void *vptr = 0;
-    int res = SWIG_ConvertPtr(argv[0], &vptr, SWIGTYPE_p_SimulationResult, 0);
-    _v = SWIG_CheckState(res);
-    if (_v) {
-      {
-        int res = SWIG_AsVal_size_t(argv[1], NULL);
-        _v = SWIG_CheckState(res);
-      }
-      if (_v) {
-        return _wrap_SimulationResult_convertedBinCenters__SWIG_1(self, argc, argv);
-      }
-    }
-  }
-  
-fail:
-  SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'SimulationResult_convertedBinCenters'.\n"
-    "  Possible C/C++ prototypes are:\n"
-    "    SimulationResult::convertedBinCenters() const\n"
-    "    SimulationResult::convertedBinCenters(size_t) const\n");
-  return 0;
-}
-
-
 SWIGINTERN PyObject *_wrap_SimulationResult_setTitle(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   SimulationResult *arg1 = (SimulationResult *) 0 ;
@@ -38082,7 +37920,6 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
 	 { "new_Datafield", _wrap_new_Datafield, METH_VARARGS, "\n"
-		"Datafield()\n"
 		"Datafield(Frame frame, vdouble1d_t values={}, vdouble1d_t errSigmas={})\n"
 		"Datafield(std::vector< Scale const *,std::allocator< Scale const * > > && axes, vdouble1d_t values={}, vdouble1d_t errSigmas={})\n"
 		"Datafield(Datafield arg1)\n"
@@ -38388,20 +38225,14 @@ static PyMethodDef SwigMethods[] = {
 	 { "writeDatafield", _wrap_writeDatafield, METH_VARARGS, "writeDatafield(Datafield data, std::string const & file_name)"},
 	 { "dataMatchesFile", _wrap_dataMatchesFile, METH_VARARGS, "dataMatchesFile(Datafield data, std::string const & refFileName, double tol) -> bool"},
 	 { "new_SimulationResult", _wrap_new_SimulationResult, METH_VARARGS, "\n"
-		"SimulationResult()\n"
-		"SimulationResult(Datafield data, Frame coords)\n"
+		"SimulationResult(Datafield data)\n"
 		"SimulationResult(SimulationResult other)\n"
 		"new_SimulationResult(SimulationResult other) -> SimulationResult\n"
 		""},
 	 { "delete_SimulationResult", _wrap_delete_SimulationResult, METH_O, "delete_SimulationResult(SimulationResult self)"},
 	 { "SimulationResult_extracted_field", _wrap_SimulationResult_extracted_field, METH_O, "SimulationResult_extracted_field(SimulationResult self) -> Datafield"},
 	 { "SimulationResult_axisMinMax", _wrap_SimulationResult_axisMinMax, METH_VARARGS, "SimulationResult_axisMinMax(SimulationResult self, size_t i) -> pvacuum_double_t"},
-	 { "SimulationResult_name_of_axis", _wrap_SimulationResult_name_of_axis, METH_VARARGS, "SimulationResult_name_of_axis(SimulationResult self, size_t i) -> std::string"},
-	 { "SimulationResult_array", _wrap_SimulationResult_array, METH_O, "SimulationResult_array(SimulationResult self) -> PyObject *"},
-	 { "SimulationResult_convertedBinCenters", _wrap_SimulationResult_convertedBinCenters, METH_VARARGS, "\n"
-		"SimulationResult_convertedBinCenters(SimulationResult self) -> vdouble1d_t\n"
-		"SimulationResult_convertedBinCenters(SimulationResult self, size_t i_axis) -> vdouble1d_t\n"
-		""},
+	 { "SimulationResult_convertedBinCenters", _wrap_SimulationResult_convertedBinCenters, METH_VARARGS, "SimulationResult_convertedBinCenters(SimulationResult self, size_t i_axis) -> vdouble1d_t"},
 	 { "SimulationResult_setTitle", _wrap_SimulationResult_setTitle, METH_VARARGS, "SimulationResult_setTitle(SimulationResult self, std::string const & title)"},
 	 { "SimulationResult_title", _wrap_SimulationResult_title, METH_O, "SimulationResult_title(SimulationResult self) -> std::string"},
 	 { "SimulationResult_swigregister", SimulationResult_swigregister, METH_O, NULL},
diff --git a/auto/Wrap/libBornAgainSim.py b/auto/Wrap/libBornAgainSim.py
index 0949d4e806ddc8c5c5ce498b9e32780de389c898..f26346adaa4c4a4c3cc263af011e14c3e26499e8 100644
--- a/auto/Wrap/libBornAgainSim.py
+++ b/auto/Wrap/libBornAgainSim.py
@@ -2050,6 +2050,7 @@ class vector_R3(object):
 _libBornAgainSim.vector_R3_swigregister(vector_R3)
 import libBornAgainParam
 import libBornAgainSample
+import libBornAgainDevice
 class swig_dummy_type_const_inode_vector(object):
     r"""Proxy of C++ std::vector< INode const * > class."""
 
diff --git a/auto/Wrap/libBornAgainSim_wrap.cpp b/auto/Wrap/libBornAgainSim_wrap.cpp
index 35781db7296f223f010cf1dfc056691742897a18..15f8ea99be2125a7c08078a1b0b7c8a6d0be3b4f 100644
--- a/auto/Wrap/libBornAgainSim_wrap.cpp
+++ b/auto/Wrap/libBornAgainSim_wrap.cpp
@@ -7021,6 +7021,7 @@ SWIGINTERN void std_vector_Sl_std_pair_Sl_double_Sc_double_Sg__Sg__insert__SWIG_
 #include "Base/Axis/Scale.h"
 #include "Base/Axis/Frame.h"
 
+#include "Device/Data/Datafield.h"
 #include "Device/Histo/SimulationResult.h"
 
 #include "Fit/Minimizer/MinimizerResult.h"
@@ -30894,7 +30895,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_simulationResult__SWIG_0(PyObject *self,
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -30920,7 +30921,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_simulationResult__SWIG_1(PyObject *self,
   FitObjective *arg1 = (FitObjective *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -30986,7 +30987,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_experimentalData__SWIG_0(PyObject *self,
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31012,7 +31013,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_experimentalData__SWIG_1(PyObject *self,
   FitObjective *arg1 = (FitObjective *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31078,7 +31079,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_uncertaintyData_cpp__SWIG_0(PyObject *se
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31104,7 +31105,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_uncertaintyData_cpp__SWIG_1(PyObject *se
   FitObjective *arg1 = (FitObjective *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31170,7 +31171,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_relativeDifference__SWIG_0(PyObject *sel
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31196,7 +31197,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_relativeDifference__SWIG_1(PyObject *sel
   FitObjective *arg1 = (FitObjective *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31262,7 +31263,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_absoluteDifference__SWIG_0(PyObject *sel
   int res1 = 0 ;
   size_t val2 ;
   int ecode2 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 2) || (nobjs > 2)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -31288,7 +31289,7 @@ SWIGINTERN PyObject *_wrap_FitObjective_absoluteDifference__SWIG_1(PyObject *sel
   FitObjective *arg1 = (FitObjective *) 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if ((nobjs < 1) || (nobjs > 1)) SWIG_fail;
   res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_FitObjective, 0 |  0 );
@@ -33494,7 +33495,7 @@ SWIGINTERN PyObject *_wrap_ISimulation_simulate(PyObject *self, PyObject *args)
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
-  SimulationResult result;
+  SwigValueWrapper< SimulationResult > result;
   
   if (!args) SWIG_fail;
   swig_obj[0] = args;
diff --git a/rawEx/fit/scatter2d/fit2d.py b/rawEx/fit/scatter2d/fit2d.py
index b0c1383f431474434381494fc7f85b7985197014..fa876429bd968a0787c10f6efe7886472e5e9044 100755
--- a/rawEx/fit/scatter2d/fit2d.py
+++ b/rawEx/fit/scatter2d/fit2d.py
@@ -34,7 +34,7 @@ def real_data():
     """
     simulation = get_simulation({'radius': 5 * nm})
     result = simulation.simulate()
-    return result.array()
+    return result.npArray()
 
 
 if __name__ == '__main__':
diff --git a/rawEx/fit/scatter2d/model2_hexlattice.py b/rawEx/fit/scatter2d/model2_hexlattice.py
index e88dc4cf5818f26966596153d295e9198a353ccb..1ce38e56f04cfde2379b0291e7d574a25c861109 100644
--- a/rawEx/fit/scatter2d/model2_hexlattice.py
+++ b/rawEx/fit/scatter2d/model2_hexlattice.py
@@ -57,7 +57,7 @@ def create_real_data():
     result = simulation.simulate()
 
     # retrieving simulated data in the form of numpy array
-    real_data = result.array()
+    real_data = result.npArray()
 
     # spoiling simulated data with noise to produce "real" data
     np.random.seed(0)
diff --git a/rawEx/fit/specular/Honeycomb_fit.py b/rawEx/fit/specular/Honeycomb_fit.py
index e9dc5fa803305e2f3892401c3e4d4c717664c62a..f8f9a96073799af461e0c7037003be3cffb0da1a 100755
--- a/rawEx/fit/specular/Honeycomb_fit.py
+++ b/rawEx/fit/specular/Honeycomb_fit.py
@@ -109,7 +109,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = numpy.array(result.convertedBinCenters())
+    q = numpy.array(result.convertedBinCenters(0))
     r = numpy.array(result.array())
 
     return q, r
diff --git a/rawEx/fit/specular/Pt_layer_fit.py b/rawEx/fit/specular/Pt_layer_fit.py
index de39eb13832791629717cc2040e15e6a41b1357b..170527753bf1d6f747602d8bd058eb69b21bbe20 100755
--- a/rawEx/fit/specular/Pt_layer_fit.py
+++ b/rawEx/fit/specular/Pt_layer_fit.py
@@ -88,7 +88,7 @@ def qr(result):
     """
     Return q and reflectivity arrays from simulation result.
     """
-    q = np.array(result.convertedBinCenters())
+    q = np.array(result.convertedBinCenters(0))
     r = np.array(result.array())
 
     return q, r
diff --git a/rawEx/specular/PolarizedSpinAsymmetry.py b/rawEx/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..01e333315874c14b27efde71b556f2ec51a28da3 100755
--- a/rawEx/specular/PolarizedSpinAsymmetry.py
+++ b/rawEx/specular/PolarizedSpinAsymmetry.py
@@ -101,8 +101,8 @@ def qr(result):
     Returns two arrays that hold the q-values as well as the
     reflectivity from a given simulation result
     """
-    q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    q = numpy.array(result.convertedBinCenters(0))
+    r = numpy.array(result.npArray())
 
     return q, r