diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index 53e84ec390b64bdb25927c4d60909588cb9a4310..b47438f97c1415049f3d4178f3c79cff0815cb3e 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -127,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;
@@ -139,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 2ee7cfa70cbdc88e4f12bc08af36c548b7cd1284..df0a752f4662bfb6f2781a0fba423abe5c03f3b9 100644
--- a/Device/Histo/SimulationResult.cpp
+++ b/Device/Histo/SimulationResult.cpp
@@ -49,21 +49,6 @@ 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);
diff --git a/Device/Histo/SimulationResult.h b/Device/Histo/SimulationResult.h
index 782a2d21a0b0aa22df6a50e71d1bc43b386792a9..68d8d1d10ef0034b948d92f62d9be02de6b0f035 100644
--- a/Device/Histo/SimulationResult.h
+++ b/Device/Histo/SimulationResult.h
@@ -38,11 +38,6 @@ public:
 
     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
@@ -53,10 +48,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/Tests/Py/Functional/PolarizedNoAnalyzer.py b/Tests/Py/Functional/PolarizedNoAnalyzer.py
index 535407dad12e6a61c6f0e773a980c8d3619d3775..4b2e8a2760f7ffd697e44df3bd39566e3acfbe12 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(), result.npArray()
 
 
 if __name__ == '__main__':
diff --git a/Wrap/Python/ba_fitmonitor.py b/Wrap/Python/ba_fitmonitor.py
index 7d1f9b0149d9972b603d4befbefbb7c2545ce619..fa19fc8464d29a6cc1b0bc37d1729a31f3599711 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,9 +155,9 @@ 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 }
diff --git a/Wrap/Python/ba_plot.py b/Wrap/Python/ba_plot.py
index 1e456d7dd6876ab86756dc234b66cb8543d49d76..1d405d86e49bbf3ba622dcce5c953f4708e11e90 100644
--- a/Wrap/Python/ba_plot.py
+++ b/Wrap/Python/ba_plot.py
@@ -195,7 +195,7 @@ def plot_specular_curve(result):
     :param result: SimulationResult from SpecularSimulation
     Used internally.
     """
-    intensity = result.array()
+    intensity = result.npArray()
     x_axis = result.convertedBinCenters()
 
     xlabel = plotargs.pop('xlabel', get_axes_labels(result)[0])
@@ -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,
@@ -522,7 +522,7 @@ def plot_multicurve(results, xlabel, ylabel):
     legend = []
     for result in results:
         x = result.convertedBinCenters()
-        y = result.array()
+        y = result.npArray()
         legend.append(result.title())
         plt.plot(x, y)
 
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/PolarizedSpinAsymmetry.py b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..ebd700d4f500eacf1d831938988a6164c0d52257 100755
--- a/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/fit/specular/PolarizedSpinAsymmetry.py
@@ -102,7 +102,7 @@ def qr(result):
     reflectivity from a given simulation result
     """
     q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/Examples/specular/PolarizedSpinAsymmetry.py b/auto/Examples/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..ebd700d4f500eacf1d831938988a6164c0d52257 100755
--- a/auto/Examples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/Examples/specular/PolarizedSpinAsymmetry.py
@@ -102,7 +102,7 @@ def qr(result):
     reflectivity from a given simulation result
     """
     q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    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/PolarizedSpinAsymmetry.py b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..ebd700d4f500eacf1d831938988a6164c0d52257 100755
--- a/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/fit/specular/PolarizedSpinAsymmetry.py
@@ -102,7 +102,7 @@ def qr(result):
     reflectivity from a given simulation result
     """
     q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..ebd700d4f500eacf1d831938988a6164c0d52257 100755
--- a/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
+++ b/auto/MiniExamples/specular/PolarizedSpinAsymmetry.py
@@ -102,7 +102,7 @@ def qr(result):
     reflectivity from a given simulation result
     """
     q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    r = numpy.array(result.npArray())
 
     return q, r
 
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index faefe2cbfd7f42a928a4d293266f51ab839321a9..5c2ab8d8c4fa04b3c9e94f74057501209b525d4f 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -3061,10 +3061,6 @@ class SimulationResult(Datafield):
         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
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index 6ec966e32b282fac9795a281cab2d15553fda37b..277b9305f4c4957454eb635fb1a93614c90ac335 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -37157,29 +37157,6 @@ fail:
 }
 
 
-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 ;
@@ -38348,7 +38325,6 @@ static PyMethodDef SwigMethods[] = {
 	 { "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"
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/specular/PolarizedSpinAsymmetry.py b/rawEx/specular/PolarizedSpinAsymmetry.py
index 5f4235c6d2924485ad136ddbe1e60afc881522b6..ebd700d4f500eacf1d831938988a6164c0d52257 100755
--- a/rawEx/specular/PolarizedSpinAsymmetry.py
+++ b/rawEx/specular/PolarizedSpinAsymmetry.py
@@ -102,7 +102,7 @@ def qr(result):
     reflectivity from a given simulation result
     """
     q = numpy.array(result.convertedBinCenters())
-    r = numpy.array(result.array())
+    r = numpy.array(result.npArray())
 
     return q, r