diff --git a/Device/Analyze/Fourier.cpp b/Device/Analyze/Fourier.cpp
index b3ec5f790854da5c1f9bf10153660c541876a0ee..232d774a30a8c3febeefc6971cb6c9f22dc4c884 100644
--- a/Device/Analyze/Fourier.cpp
+++ b/Device/Analyze/Fourier.cpp
@@ -24,7 +24,7 @@ Datafield Analyze::createFFT(const Datafield& data)
     const auto signal = data.values2D();
 
     FourierTransform ft;
-    std::vector<std::vector<double>> signal2 = ft.ramplitude(signal);
+    double2d_t signal2 = ft.ramplitude(signal);
     signal2 = ft.fftshift(signal2); // low frequency to center of array
 
     return {"~" + data.xAxis().axisLabel(), "~" + data.yAxis().axisLabel(), signal2};
diff --git a/Device/Analyze/Peaks.cpp b/Device/Analyze/Peaks.cpp
index fc012a565454917a43301d7604b181f641f0b36b..f0d8078dd88959cb6e7cb2ee7776eb10596ad6cd 100644
--- a/Device/Analyze/Peaks.cpp
+++ b/Device/Analyze/Peaks.cpp
@@ -21,7 +21,7 @@
 std::vector<std::pair<double, double>>
 Analyze::FindPeaks(const Datafield& data, double sigma, const std::string& option, double threshold)
 {
-    const std::vector<std::vector<double>> arr = data.values2D();
+    const double2d_t arr = data.values2D();
     tspectrum::Spectrum2D spec;
     const auto peaks = spec.find_peaks(arr, sigma, option, threshold);
 
diff --git a/Device/Data/DataUtil.cpp b/Device/Data/DataUtil.cpp
index e842c755eca50f7aaf4784b88787cdf9694e5520..06db72b4701acbc84e91b3c875bb9eada7e726ad 100644
--- a/Device/Data/DataUtil.cpp
+++ b/Device/Data/DataUtil.cpp
@@ -21,10 +21,9 @@
 #include <algorithm>
 #include <cmath>
 
-std::vector<std::vector<double>>
-DataUtil::invertAxis(int axis, const std::vector<std::vector<double>>& original)
+double2d_t DataUtil::invertAxis(int axis, const double2d_t& original)
 {
-    std::vector<std::vector<double>> inverse = original;
+    double2d_t inverse = original;
 
     const size_t orig_rows = original.size();
 
@@ -43,15 +42,14 @@ DataUtil::invertAxis(int axis, const std::vector<std::vector<double>>& original)
     return inverse;
 }
 
-std::vector<std::vector<double>>
-DataUtil::transpose(const std::vector<std::vector<double>>& original)
+double2d_t DataUtil::transpose(const double2d_t& original)
 {
     ASSERT(!original.empty());
 
     const size_t orig_rows = original.size();
     const size_t orig_cols = original.front().size();
 
-    std::vector<std::vector<double>> transposed(orig_cols, std::vector<double>(orig_rows));
+    double2d_t transposed(orig_cols, std::vector<double>(orig_rows));
 
     for (size_t i = 0; i < orig_rows; ++i)
         for (size_t j = 0; j < orig_cols; ++j)
diff --git a/Device/Data/DataUtil.h b/Device/Data/DataUtil.h
index 5b30b829259992667d99426726f2d5b3a3fa324c..fe7dee302f87c41f30c7ea3838feb41c5bc98c59 100644
--- a/Device/Data/DataUtil.h
+++ b/Device/Data/DataUtil.h
@@ -18,6 +18,7 @@
 #ifndef BORNAGAIN_DEVICE_DATA_DATAUTIL_H
 #define BORNAGAIN_DEVICE_DATA_DATAUTIL_H
 
+#include "Base/Type/Field2D.h"
 #include <string>
 #include <vector>
 
@@ -32,14 +33,13 @@ namespace DataUtil {
 Datafield rotatedDatafield(const Datafield& data, int n);
 
 //! Changes the direction of rows or columns
-std::vector<std::vector<double>> invertAxis(int axis,
-                                            const std::vector<std::vector<double>>& original);
+double2d_t invertAxis(int axis, const double2d_t& original);
 
 //! Transposes the matrix
-std::vector<std::vector<double>> transpose(const std::vector<std::vector<double>>& original);
+double2d_t transpose(const double2d_t& original);
 
 //! Creates a vector of vectors of double (2D Array) from Datafield.
-std::vector<std::vector<double>> create2DArrayfromDatafield(const Datafield& data);
+double2d_t create2DArrayfromDatafield(const Datafield& data);
 
 //! Returns Datafield representing relative difference of two histograms.
 Datafield relativeDifferenceField(const Datafield& dat, const Datafield& ref);
diff --git a/Device/Data/Datafield.cpp b/Device/Data/Datafield.cpp
index 4655dd07783ef7c5e9ef7f5e2ae449044fa70766..08763d2e91369beba4b1f404806a2ac759d40203 100644
--- a/Device/Data/Datafield.cpp
+++ b/Device/Data/Datafield.cpp
@@ -85,8 +85,7 @@ Datafield::Datafield(const std::vector<const Scale*>& axes)
 {
 }
 
-Datafield::Datafield(const std::string& xlabel, const std::string& ylabel,
-                     const std::vector<std::vector<double>>& vec)
+Datafield::Datafield(const std::string& xlabel, const std::string& ylabel, const double2d_t& vec)
 {
     const size_t nrows = vec.size();
     const size_t ncols = vec[0].size();
@@ -165,7 +164,7 @@ void Datafield::setVector(const std::vector<double>& vector)
     m_values = vector;
 }
 
-void Datafield::setVector2D(const std::vector<std::vector<double>>& in)
+void Datafield::setVector2D(const double2d_t& in)
 {
     const size_t ncols = axis(0).size();
     const size_t nrows = axis(1).size();
@@ -422,10 +421,10 @@ Datafield Datafield::normalizedToMax() const
     return {title(), frame().clone(), outval, errval};
 }
 
-std::vector<std::vector<double>> Datafield::values2D() const
+double2d_t Datafield::values2D() const
 {
     ASSERT(rank() == 2);
-    std::vector<std::vector<double>> result;
+    double2d_t result;
 
     const size_t nrows = axis(1).size();
     const size_t ncols = axis(0).size();
diff --git a/Device/Data/Datafield.h b/Device/Data/Datafield.h
index c50e40f3898a907e0e0d268172a76c7cc6efce8c..6a900a570b87c85f5f8d4d8cafcc6074a3a0053c 100644
--- a/Device/Data/Datafield.h
+++ b/Device/Data/Datafield.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 #define BORNAGAIN_DEVICE_DATA_DATAFIELD_H
 
+#include "Base/Type/Field2D.h"
 #include <memory>
 #include <string>
 #include <vector>
@@ -51,8 +52,7 @@ public:
     Datafield(const std::vector<const Scale*>& axes);
 
     //! Constructor that creates equidistant grid from 2D array.
-    Datafield(const std::string& xlabel, const std::string& ylabel,
-              const std::vector<std::vector<double>>& vec);
+    Datafield(const std::string& xlabel, const std::string& ylabel, const double2d_t& vec);
 
     Datafield(const Datafield&);
 #ifndef SWIG
@@ -155,7 +155,7 @@ public:
 
     //! Sets new values to raw data vector.
     void setVector(const std::vector<double>& data_vector);
-    void setVector2D(const std::vector<std::vector<double>>& in);
+    void setVector2D(const double2d_t& in);
 
     bool hasErrorSigmas() const;
     std::vector<double>& errorSigmas()
@@ -163,7 +163,7 @@ public:
         return m_err_sigmas;
     }
 
-    std::vector<std::vector<double>> values2D() const;
+    double2d_t values2D() const;
 
 private:
     std::string m_title;
diff --git a/Device/IO/ReadReflectometry.cpp b/Device/IO/ReadReflectometry.cpp
index 7807a794440d5d7f814b97c1734ba77bc483df95..333f9bcaeca0b195ffb41376cce909cf06e3f8cd 100644
--- a/Device/IO/ReadReflectometry.cpp
+++ b/Device/IO/ReadReflectometry.cpp
@@ -37,7 +37,7 @@ Datafield Util::RW::readReflectometryTable(std::istream& s, const ImportSettings
     std::string line;
     int lineno = 0;
 
-    std::vector<std::vector<double>> rowsVec;
+    double2d_t rowsVec;
 
     // Read numbers from file:
     while (std::getline(s, line)) {
diff --git a/Device/IO/ReadWrite2DTable.cpp b/Device/IO/ReadWrite2DTable.cpp
index 51282614b670fc2377305e87bd6d67120e20a835..3d95d18a79bbb0d934f1ce79de243983d34ecaef 100644
--- a/Device/IO/ReadWrite2DTable.cpp
+++ b/Device/IO/ReadWrite2DTable.cpp
@@ -56,7 +56,7 @@ void write2DRepresentation(const Datafield& data, std::ostream& output_stream)
 
     output_stream << "# [nrows=" << nrows << ", ncols=" << ncols << "]" << std::endl;
 
-    std::vector<std::vector<double>> dataArray = data.values2D();
+    double2d_t dataArray = data.values2D();
     output_stream.imbue(std::locale::classic());
     output_stream << std::scientific << std::setprecision(12);
 
@@ -83,10 +83,10 @@ bool getNextLine(std::istream& input_stream, std::string& line)
     return false;
 }
 
-std::vector<std::vector<double>> parseFile(std::istream& input_stream)
+double2d_t parseFile(std::istream& input_stream)
 {
     std::string line;
-    std::vector<std::vector<double>> data;
+    double2d_t data;
 
     // Read numbers from input stream:
     size_t nrows = 0;
@@ -105,7 +105,7 @@ std::vector<std::vector<double>> parseFile(std::istream& input_stream)
 
 Datafield readBareIntensity(std::istream& input_stream)
 {
-    std::vector<std::vector<double>> data = parseFile(input_stream);
+    double2d_t data = parseFile(input_stream);
     size_t nrows = data.size();
     size_t ncols = nrows ? data[0].size() : 0;
 
@@ -143,7 +143,7 @@ Datafield Util::RW::read2DTable(std::istream& input_stream, const ImportSettings
         return readBareIntensity(input_stream);
 
     // read table with axes info
-    std::vector<std::vector<double>> data = parseFile(input_stream);
+    double2d_t data = parseFile(input_stream);
     size_t nrows = data.size();
     size_t ncols = nrows ? data.front().size() : 0;
 
diff --git a/Device/Resolution/ConvolutionDetectorResolution.cpp b/Device/Resolution/ConvolutionDetectorResolution.cpp
index ea45fc45bc123a58057f763ccb5f51194659675c..016f848156aa8a485cf3dea3ba1be76b59078f71 100644
--- a/Device/Resolution/ConvolutionDetectorResolution.cpp
+++ b/Device/Resolution/ConvolutionDetectorResolution.cpp
@@ -109,7 +109,7 @@ void ConvolutionDetectorResolution::apply2dConvolution(Datafield* df) const
         return;
 
     // Construct kernel vector from resolution function
-    std::vector<std::vector<double>> kernel;
+    double2d_t kernel;
     kernel.resize(ny);
     double mid_value1 = X.binCenter(nx / 2); // because Convolve expects zero at midpoint
     double mid_value2 = Y.binCenter(ny / 2); // because Convolve expects zero at midpoint
@@ -127,7 +127,7 @@ void ConvolutionDetectorResolution::apply2dConvolution(Datafield* df) const
     }
 
     // Calculate convolution
-    std::vector<std::vector<double>> result;
+    double2d_t result;
     Convolve().fftconvolve2D(df->values2D(), kernel, result);
 
     df->setVector2D(result);
diff --git a/Fit/Adapter/MinimizerAdapter.cpp b/Fit/Adapter/MinimizerAdapter.cpp
index ea0cbfbca404604773a514ac42ebab3956ab38e3..0482d7004b3c792c45d4d0231bf302500b069e58 100644
--- a/Fit/Adapter/MinimizerAdapter.cpp
+++ b/Fit/Adapter/MinimizerAdapter.cpp
@@ -124,7 +124,7 @@ void MinimizerAdapter::propagateResults(mumufit::Parameters& parameters)
     parameters.setErrors(parErrorsAtMinimum());
     // sets correlation matrix
     if (providesError()) {
-        mumufit::Parameters::corr_matrix_t matrix;
+        double2d_t matrix;
         matrix.resize(fitRank());
 
         for (size_t i = 0; i < fitRank(); ++i) {
diff --git a/Fit/Minimizer/MinimizerResult.cpp b/Fit/Minimizer/MinimizerResult.cpp
index 388fdda5c5552edb82f7046b5401416655f69ec9..fa5399ade6e87503bafd9763d2fb440bf7d00b49 100644
--- a/Fit/Minimizer/MinimizerResult.cpp
+++ b/Fit/Minimizer/MinimizerResult.cpp
@@ -33,7 +33,7 @@ std::string reportParameters(const mumufit::Parameters& parameters)
                       % par.startValue() % par.limits().toString() % par.value() % par.error();
     }
 
-    mumufit::Parameters::corr_matrix_t matrix = parameters.correlationMatrix();
+    double2d_t matrix = parameters.correlationMatrix();
     if (!matrix.empty()) {
         result << mumufit::utils::sectionString("Correlations");
         for (size_t i = 0; i < matrix.size(); ++i) {
diff --git a/Fit/Param/Parameters.cpp b/Fit/Param/Parameters.cpp
index 12df12a211610e381f6fff9a9184366f3572af2a..3f1c20d548e4be3f478d2ad84b8806dfb3b2643c 100644
--- a/Fit/Param/Parameters.cpp
+++ b/Fit/Param/Parameters.cpp
@@ -114,12 +114,12 @@ const Parameter& Parameters::operator[](size_t index) const
     return m_parameters[check_index(index)];
 }
 
-Parameters::corr_matrix_t Parameters::correlationMatrix() const
+double2d_t Parameters::correlationMatrix() const
 {
     return m_corr_matrix;
 }
 
-void Parameters::setCorrelationMatrix(const Parameters::corr_matrix_t& matrix)
+void Parameters::setCorrelationMatrix(const double2d_t& matrix)
 {
     if (matrix.size() != size())
         throw std::runtime_error("Parameters::setCorrelationMatrix -> Error. Wrong "
diff --git a/Fit/Param/Parameters.h b/Fit/Param/Parameters.h
index f71f197361f3d25b3f4e93cfc0643156204e5698..f038fc4fae402854624038a40e8602ea1ebffb83 100644
--- a/Fit/Param/Parameters.h
+++ b/Fit/Param/Parameters.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_FIT_PARAM_PARAMETERS_H
 #define BORNAGAIN_FIT_PARAM_PARAMETERS_H
 
+#include "Base/Type/Field2D.h"
 #include "Fit/Param/Parameter.h"
 #include <vector>
 
@@ -27,7 +28,6 @@ public:
     using parameters_t = std::vector<Parameter>;
     using const_iterator = parameters_t::const_iterator;
     using iterator = parameters_t::iterator;
-    using corr_matrix_t = std::vector<std::vector<double>>;
 
     Parameters() = default;
 
@@ -50,8 +50,8 @@ public:
     const Parameter& operator[](const std::string& name) const;
     const Parameter& operator[](size_t index) const;
 
-    corr_matrix_t correlationMatrix() const;
-    void setCorrelationMatrix(const corr_matrix_t& matrix);
+    double2d_t correlationMatrix() const;
+    void setCorrelationMatrix(const double2d_t& matrix);
 
     size_t freeParameterCount() const;
 
@@ -61,7 +61,7 @@ private:
     size_t check_index(size_t index) const;
 
     parameters_t m_parameters;
-    corr_matrix_t m_corr_matrix; //!< correlation matrix
+    double2d_t m_corr_matrix; //!< correlation matrix
 };
 
 } // namespace mumufit
diff --git a/Fit/Residual/ResidualFunctionAdapter.h b/Fit/Residual/ResidualFunctionAdapter.h
index d696eec24f161f2ad225f7bf03c8e5e8f4c6eb6f..ebba0b3cb2dd5a6505abfc74d6559e9614f6f2e5 100644
--- a/Fit/Residual/ResidualFunctionAdapter.h
+++ b/Fit/Residual/ResidualFunctionAdapter.h
@@ -53,7 +53,7 @@ private:
     fcn_residual_t m_fcn; //!< user function to minimize
     Parameters m_parameters;
     std::vector<double> m_residuals;
-    std::vector<std::vector<double>> m_gradients; // [m_npars][m_ndatasize]
+    double2d_t m_gradients; // [m_npars][m_ndatasize]
     std::unique_ptr<RootResidualFunction> m_root_objective;
 };
 
diff --git a/GUI/View/Realspace/RealspaceBuilder.cpp b/GUI/View/Realspace/RealspaceBuilder.cpp
index 267770e628b2614b2b2944b0dcdedefe111f3a8f..ed2d4a3038ac9fef655f315ce825782153f2065b 100644
--- a/GUI/View/Realspace/RealspaceBuilder.cpp
+++ b/GUI/View/Realspace/RealspaceBuilder.cpp
@@ -314,7 +314,7 @@ void RealspaceBuilder::translateContainer(Model* model,
 }
 
 void RealspaceBuilder::populateParticlesInLattice(
-    const std::vector<std::vector<double>>& lattice_positions,
+    const double2d_t& lattice_positions,
     const std::vector<Particle3DContainer>& particle3DContainer_vector, Model* model,
     const SceneGeometry& sceneGeometry, unsigned& numParticles) const
 {
diff --git a/GUI/View/Realspace/RealspaceBuilder.h b/GUI/View/Realspace/RealspaceBuilder.h
index 45ae620b6ae4bb1482f4067cdce45b767c87ac99..858af20c990339d1a497f515e6bfb200173b30b8 100644
--- a/GUI/View/Realspace/RealspaceBuilder.h
+++ b/GUI/View/Realspace/RealspaceBuilder.h
@@ -69,7 +69,7 @@ private:
                             unsigned& numParticles, const Img3D::F3& lattice_position = {}) const;
 
     void populateParticlesInLattice(
-        const std::vector<std::vector<double>>& lattice_positions,
+        const double2d_t& lattice_positions,
         const std::vector<Img3D::Particle3DContainer>& particle3DContainer_vector,
         Img3D::Model* model, const SceneGeometry& sceneGeometry, unsigned& numParticles) const;
 };
diff --git a/Img3D/Build/ParacrystalLatticePositions.cpp b/Img3D/Build/ParacrystalLatticePositions.cpp
index 8c4553ffe4631e48238c02175f7bb34430694e75..e4730cbfb6dfc6f891fbb50b89d742b15074ff7c 100644
--- a/Img3D/Build/ParacrystalLatticePositions.cpp
+++ b/Img3D/Build/ParacrystalLatticePositions.cpp
@@ -18,8 +18,7 @@
 
 namespace {
 
-void ResizeLatticePositions(std::vector<std::vector<double>>& lattice_positions, double l1,
-                            double l2, double layer_size)
+void ResizeLatticePositions(double2d_t& lattice_positions, double l1, double l2, double layer_size)
 {
     // Estimate the limit n1 and n2 of the integer multiple j and i of the lattice vectors l1 and l2
     // required for populating particles correctly within the 3D model's boundaries
@@ -58,9 +57,10 @@ void FindLatticePositionsIndex(size_t& index, size_t& index_prev, int i, int j,
     }
 }
 
-std::pair<double, double> ComputePositionAlongPositiveLatticeVector(
-    const size_t index_prev, std::vector<std::vector<double>>& lattice_positions,
-    const IProfile2D* pdf, double l, double l_xi, double l_alpha)
+std::pair<double, double> ComputePositionAlongPositiveLatticeVector(const size_t index_prev,
+                                                                    double2d_t& lattice_positions,
+                                                                    const IProfile2D* pdf, double l,
+                                                                    double l_xi, double l_alpha)
 {
     double gamma_pdf = pdf->gamma();
     std::pair<double, double> sampleXYpdf = pdf->createSampler()->randomSample();
@@ -78,9 +78,10 @@ std::pair<double, double> ComputePositionAlongPositiveLatticeVector(
     return std::make_pair(x, y);
 }
 
-std::pair<double, double> ComputePositionAlongNegativeLatticeVector(
-    const size_t index_prev, std::vector<std::vector<double>>& lattice_positions,
-    const IProfile2D* pdf, double l, double l_xi, double l_alpha)
+std::pair<double, double> ComputePositionAlongNegativeLatticeVector(const size_t index_prev,
+                                                                    double2d_t& lattice_positions,
+                                                                    const IProfile2D* pdf, double l,
+                                                                    double l_xi, double l_alpha)
 {
     double gamma_pdf = pdf->gamma();
     std::pair<double, double> sampleXYpdf = pdf->createSampler()->randomSample();
@@ -98,10 +99,10 @@ std::pair<double, double> ComputePositionAlongNegativeLatticeVector(
     return std::make_pair(x, y);
 }
 
-std::pair<double, double>
-ComputeLatticePosition(const size_t index_prev, int i, int j,
-                       std::vector<std::vector<double>>& lattice_positions, const IProfile2D* pdf,
-                       double l, double l_xi, double l_alpha)
+std::pair<double, double> ComputeLatticePosition(const size_t index_prev, int i, int j,
+                                                 double2d_t& lattice_positions,
+                                                 const IProfile2D* pdf, double l, double l_xi,
+                                                 double l_alpha)
 {
     if (std::sin(l_alpha) == 0) {
         if (!(j % 2 == 0)) // along +l1
@@ -119,9 +120,8 @@ ComputeLatticePosition(const size_t index_prev, int i, int j,
                                                      l_alpha);
 }
 
-void ComputePositionsAlongLatticeVectorAxes(std::vector<std::vector<double>>& lattice_positions,
-                                            const IProfile2D* pdf, double l, double l_xi,
-                                            double l_alpha)
+void ComputePositionsAlongLatticeVectorAxes(double2d_t& lattice_positions, const IProfile2D* pdf,
+                                            double l, double l_xi, double l_alpha)
 {
     int n = static_cast<int>((std::sqrt(lattice_positions.size()) - 1) / 2);
 
@@ -156,9 +156,9 @@ void ComputePositionsAlongLatticeVectorAxes(std::vector<std::vector<double>>& la
     }
 }
 
-void ComputePositionsInsideLatticeQuadrants(std::vector<std::vector<double>>& lattice_positions,
-                                            const IProfile2D* pdf1, const IProfile2D* pdf2,
-                                            double l1, double l2, double l_xi, double l_alpha)
+void ComputePositionsInsideLatticeQuadrants(double2d_t& lattice_positions, const IProfile2D* pdf1,
+                                            const IProfile2D* pdf2, double l1, double l2,
+                                            double l_xi, double l_alpha)
 {
     int n = static_cast<int>((std::sqrt(lattice_positions.size()) - 1) / 2);
 
@@ -189,8 +189,7 @@ void ComputePositionsInsideLatticeQuadrants(std::vector<std::vector<double>>& la
 //  implement namespace Img3D::Paracrystal2D
 //  ************************************************************************************************
 
-std::vector<std::vector<double>>
-Paracrystal::latticePositions(const Interference2DParacrystal* p_iff, double layer_size)
+double2d_t Paracrystal::latticePositions(const Interference2DParacrystal* p_iff, double layer_size)
 {
     const auto& lattice = p_iff->lattice();
     double l1 = lattice.length1();
@@ -198,7 +197,7 @@ Paracrystal::latticePositions(const Interference2DParacrystal* p_iff, double lay
     double alpha = lattice.latticeAngle();
     double xi = lattice.rotationAngle();
 
-    std::vector<std::vector<double>> lattice_positions;
+    double2d_t lattice_positions;
     ResizeLatticePositions(lattice_positions, l1, l2, layer_size);
 
     ComputePositionsAlongLatticeVectorAxes(lattice_positions, p_iff->pdf1(), l1, xi, 0);
diff --git a/Img3D/Build/ParacrystalLatticePositions.h b/Img3D/Build/ParacrystalLatticePositions.h
index ba89b2b9ed376ff08e12b5686f95284e6cecb863..1195cd8f1d12adafb1f2f2d9243fd163e1970ec9 100644
--- a/Img3D/Build/ParacrystalLatticePositions.h
+++ b/Img3D/Build/ParacrystalLatticePositions.h
@@ -15,14 +15,14 @@
 #ifndef BORNAGAIN_IMG3D_BUILD_PARACRYSTALLATTICEPOSITIONS_H
 #define BORNAGAIN_IMG3D_BUILD_PARACRYSTALLATTICEPOSITIONS_H
 
+#include "Base/Type/Field2D.h"
 #include <vector>
 
 class Interference2DParacrystal;
 
 namespace Paracrystal {
 
-std::vector<std::vector<double>> latticePositions(const Interference2DParacrystal*,
-                                                  double layer_size);
+double2d_t latticePositions(const Interference2DParacrystal*, double layer_size);
 }
 
 #endif // BORNAGAIN_IMG3D_BUILD_PARACRYSTALLATTICEPOSITIONS_H
diff --git a/Img3D/Build/PositionBuilders.cpp b/Img3D/Build/PositionBuilders.cpp
index dabbb1ca767723b211d991e42897b958043bf115..c2be21c9ac0321df7bb0b0b225bbdcbdb60928e1 100644
--- a/Img3D/Build/PositionBuilders.cpp
+++ b/Img3D/Build/PositionBuilders.cpp
@@ -19,10 +19,10 @@
 
 namespace {
 
-std::vector<std::vector<double>> Generate2DLatticePoints(double l1, double l2, double alpha,
-                                                         double xi, unsigned n1, unsigned n2)
+double2d_t Generate2DLatticePoints(double l1, double l2, double alpha, double xi, unsigned n1,
+                                   unsigned n2)
 {
-    std::vector<std::vector<double>> lattice_positions;
+    double2d_t lattice_positions;
     std::vector<double> position;
 
     const unsigned nn1 = std::max(1u, n1);
@@ -57,10 +57,9 @@ std::vector<std::vector<double>> Generate2DLatticePoints(double l1, double l2, d
 
 IPositionBuilder::~IPositionBuilder() = default;
 
-std::vector<std::vector<double>> IPositionBuilder::generatePositions(double layer_size,
-                                                                     double density) const
+double2d_t IPositionBuilder::generatePositions(double layer_size, double density) const
 {
-    std::vector<std::vector<double>> positions = generatePositionsImpl(layer_size, density);
+    double2d_t positions = generatePositionsImpl(layer_size, density);
     const double pos_var = positionVariance();
     if (pos_var > 0.0) {
         // random generator and distribution
@@ -84,10 +83,9 @@ RandomPositionBuilder::RandomPositionBuilder() = default;
 
 RandomPositionBuilder::~RandomPositionBuilder() = default;
 
-std::vector<std::vector<double>> RandomPositionBuilder::generatePositionsImpl(double layer_size,
-                                                                              double density) const
+double2d_t RandomPositionBuilder::generatePositionsImpl(double layer_size, double density) const
 {
-    std::vector<std::vector<double>> lattice_positions;
+    double2d_t lattice_positions;
     std::vector<double> position;
 
     // to compute total number of particles we use the total particle density
@@ -127,8 +125,7 @@ Lattice1DPositionBuilder::Lattice1DPositionBuilder(const Interference1DLattice*
 
 Lattice1DPositionBuilder::~Lattice1DPositionBuilder() = default;
 
-std::vector<std::vector<double>> Lattice1DPositionBuilder::generatePositionsImpl(double layer_size,
-                                                                                 double) const
+double2d_t Lattice1DPositionBuilder::generatePositionsImpl(double layer_size, double) const
 {
     const double length = m_iff->length();
     const double xi = m_iff->xi();
@@ -158,8 +155,7 @@ Lattice2DPositionBuilder::Lattice2DPositionBuilder(const Interference2DLattice*
 
 Lattice2DPositionBuilder::~Lattice2DPositionBuilder() = default;
 
-std::vector<std::vector<double>> Lattice2DPositionBuilder::generatePositionsImpl(double layer_size,
-                                                                                 double) const
+double2d_t Lattice2DPositionBuilder::generatePositionsImpl(double layer_size, double) const
 {
     const Lattice2D& lattice = m_iff->lattice();
     const double l1 = lattice.length1();
@@ -198,8 +194,7 @@ Paracrystal2DPositionBuilder::Paracrystal2DPositionBuilder(const Interference2DP
 
 Paracrystal2DPositionBuilder::~Paracrystal2DPositionBuilder() = default;
 
-std::vector<std::vector<double>>
-Paracrystal2DPositionBuilder::generatePositionsImpl(double layer_size, double) const
+double2d_t Paracrystal2DPositionBuilder::generatePositionsImpl(double layer_size, double) const
 {
     return Paracrystal::latticePositions(m_iff.get(), layer_size);
 }
@@ -222,8 +217,7 @@ Finite2DLatticePositionBuilder::Finite2DLatticePositionBuilder(
 
 Finite2DLatticePositionBuilder::~Finite2DLatticePositionBuilder() = default;
 
-std::vector<std::vector<double>>
-Finite2DLatticePositionBuilder::generatePositionsImpl(double layer_size, double) const
+double2d_t Finite2DLatticePositionBuilder::generatePositionsImpl(double layer_size, double) const
 {
     const auto& lattice = m_iff->lattice();
     const double l1 = lattice.length1();
@@ -264,10 +258,9 @@ RadialParacrystalPositionBuilder::RadialParacrystalPositionBuilder(
 
 RadialParacrystalPositionBuilder::~RadialParacrystalPositionBuilder() = default;
 
-std::vector<std::vector<double>>
-RadialParacrystalPositionBuilder::generatePositionsImpl(double layer_size, double) const
+double2d_t RadialParacrystalPositionBuilder::generatePositionsImpl(double layer_size, double) const
 {
-    std::vector<std::vector<double>> lattice_positions;
+    double2d_t lattice_positions;
 
     const double distance = m_iff->peakDistance();
 
diff --git a/Img3D/Build/PositionBuilders.h b/Img3D/Build/PositionBuilders.h
index 20d56bcaf848d0fed9a8de71fc996ba38a3cd931..6ae39c2ee5ec831e8843cc48e2164468b26385ab 100644
--- a/Img3D/Build/PositionBuilders.h
+++ b/Img3D/Build/PositionBuilders.h
@@ -15,6 +15,7 @@
 #ifndef BORNAGAIN_IMG3D_BUILD_POSITIONBUILDERS_H
 #define BORNAGAIN_IMG3D_BUILD_POSITIONBUILDERS_H
 
+#include "Base/Type/Field2D.h"
 #include <memory>
 #include <vector>
 
@@ -28,11 +29,10 @@ class IPositionBuilder {
 public:
     virtual ~IPositionBuilder();
 
-    std::vector<std::vector<double>> generatePositions(double layer_size, double density) const;
+    double2d_t generatePositions(double layer_size, double density) const;
 
 private:
-    virtual std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                                   double density) const = 0;
+    virtual double2d_t generatePositionsImpl(double layer_size, double density) const = 0;
     virtual double positionVariance() const = 0;
 };
 
@@ -43,8 +43,7 @@ public:
     ~RandomPositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
 };
 
@@ -55,8 +54,7 @@ public:
     ~Lattice1DPositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
     std::unique_ptr<Interference1DLattice> m_iff;
 };
@@ -68,8 +66,7 @@ public:
     ~Lattice2DPositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
     std::unique_ptr<Interference2DLattice> m_iff;
 };
@@ -81,8 +78,7 @@ public:
     ~Paracrystal2DPositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
     std::unique_ptr<Interference2DParacrystal> m_iff;
 };
@@ -94,8 +90,7 @@ public:
     ~Finite2DLatticePositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
     std::unique_ptr<InterferenceFinite2DLattice> m_iff;
 };
@@ -107,8 +102,7 @@ public:
     ~RadialParacrystalPositionBuilder() override;
 
 private:
-    std::vector<std::vector<double>> generatePositionsImpl(double layer_size,
-                                                           double density) const override;
+    double2d_t generatePositionsImpl(double layer_size, double density) const override;
     double positionVariance() const override;
     std::unique_ptr<InterferenceRadialParacrystal> m_iff;
 };
diff --git a/Tests/Unit/Device/DatafieldTest.cpp b/Tests/Unit/Device/DatafieldTest.cpp
index 86c7e112a799e4101c477dd8d103215d581f1523..87c1c0e02845886ea0c94779b8dd5924c50dfcf8 100644
--- a/Tests/Unit/Device/DatafieldTest.cpp
+++ b/Tests/Unit/Device/DatafieldTest.cpp
@@ -70,14 +70,14 @@ TEST(Datafield, create2DArrayfromDatafield)
 
     auto v3 = a.values2D();
 
-    std::vector<std::vector<double>> v3expected{{0, 1, 2}, {3, 4, 5}};
+    double2d_t v3expected{{0, 1, 2}, {3, 4, 5}};
     EXPECT_EQ(v3, v3expected);
 }
 
 TEST(Datafield, ctor2D)
 {
     std::vector<double> arr_in{1, 2, 3, 4, 5, 6};
-    std::vector<std::vector<double>> array_2d{{1., 2., 3.}, {4., 5., 6.}};
+    double2d_t array_2d{{1., 2., 3.}, {4., 5., 6.}};
     const Datafield data("x ()", "y ()", array_2d);
     EXPECT_EQ(arr_in, data.flatVector());
 }
@@ -90,6 +90,6 @@ TEST(Datafield, DatafieldToVector2D)
     data.setVector(values);
 
     auto vec = data.values2D();
-    const std::vector<std::vector<double>> expected = {{0., 1., 2.}, {10., 11., 12.}};
+    const double2d_t expected = {{0., 1., 2.}, {10., 11., 12.}};
     EXPECT_EQ(vec, expected);
 }
diff --git a/Tests/Unit/Device/SpectrumTest.cpp b/Tests/Unit/Device/SpectrumTest.cpp
index 7324620c290beca414b8dd974058886e9a455f28..a916324366c7d620ecf641a61721b0973eadd3db 100644
--- a/Tests/Unit/Device/SpectrumTest.cpp
+++ b/Tests/Unit/Device/SpectrumTest.cpp
@@ -6,12 +6,12 @@
 
 TEST(Spectrum, arrayPeaks)
 {
-    std::vector<std::vector<double>> data = {{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
-                                             {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
-                                             {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
-                                             {1.0, 1.0, 1.0, 1.0, 10.0, 1.0, 1.0, 1.0, 1.0, 1.0},
-                                             {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
-                                             {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}};
+    double2d_t data = {{1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+                       {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+                       {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+                       {1.0, 1.0, 1.0, 1.0, 10.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+                       {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
+                       {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}};
 
     tspectrum::Spectrum2D finder;
     auto peaks = finder.find_peaks(data, 3, "nomarkov", 0.1);
diff --git a/Tests/Unit/Numeric/FourierTransformTest.cpp b/Tests/Unit/Numeric/FourierTransformTest.cpp
index f12544f52266b5528f0a6cb7dfc092d34e611463..abfecc54adaa653aebb19d7cefcfcf9849b8cb8d 100644
--- a/Tests/Unit/Numeric/FourierTransformTest.cpp
+++ b/Tests/Unit/Numeric/FourierTransformTest.cpp
@@ -87,12 +87,11 @@ TEST(FourierTransform, fft1D)
 // 3x5 input of all zeros
 TEST(FourierTransform, fft2DTest1)
 {
-    std::vector<std::vector<double>> signal({{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
-    std::vector<std::vector<double>> expected_result(
-        {{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
+    double2d_t signal({{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
+    double2d_t expected_result({{0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
 
     FourierTransform ft;
-    std::vector<std::vector<double>> result = ft.ramplitude(signal);
+    double2d_t result = ft.ramplitude(signal);
     result = ft.fftshift(result);
 
     for (size_t i = 0; i < signal.size(); ++i)
@@ -103,13 +102,12 @@ TEST(FourierTransform, fft2DTest1)
 // 4x5 input of all zeros except for 1 element
 TEST(FourierTransform, fft2DTest2)
 {
-    std::vector<std::vector<double>> signal(
-        {{0, 0, 0, 0, 0}, {0, 0, 2, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
-    std::vector<std::vector<double>> expected_result(
+    double2d_t signal({{0, 0, 0, 0, 0}, {0, 0, 2, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}});
+    double2d_t expected_result(
         {{2, 2, 2, 2, 2}, {2, 2, 2, 2, 2}, {2, 2, 2, 2, 2}, {2, 2, 2, 2, 2}});
 
     FourierTransform ft;
-    std::vector<std::vector<double>> result = ft.ramplitude(signal);
+    double2d_t result = ft.ramplitude(signal);
     result = ft.fftshift(result);
 
     for (size_t i = 0; i < signal.size(); ++i)
@@ -120,21 +118,21 @@ TEST(FourierTransform, fft2DTest2)
 // 6x6 input of all ones except for 1 element
 TEST(FourierTransform, fft2DTest3)
 {
-    std::vector<std::vector<double>> signal({{1, 1, 1, 1, 1, 1},
-                                             {1, 1, 1, 1, 1, 1},
-                                             {1, 1, 1, 1, 1, 1},
-                                             {1, 1, 1, 1, 1, 1},
-                                             {1, 1, 1, 1, 1, 32},
-                                             {1, 1, 1, 1, 1, 1}});
-    std::vector<std::vector<double>> expected_result({{31, 31, 31, 31, 31, 31},
-                                                      {31, 31, 31, 31, 31, 31},
-                                                      {31, 31, 31, 31, 31, 31},
-                                                      {31, 31, 31, 67, 31, 31},
-                                                      {31, 31, 31, 31, 31, 31},
-                                                      {31, 31, 31, 31, 31, 31}});
+    double2d_t signal({{1, 1, 1, 1, 1, 1},
+                       {1, 1, 1, 1, 1, 1},
+                       {1, 1, 1, 1, 1, 1},
+                       {1, 1, 1, 1, 1, 1},
+                       {1, 1, 1, 1, 1, 32},
+                       {1, 1, 1, 1, 1, 1}});
+    double2d_t expected_result({{31, 31, 31, 31, 31, 31},
+                                {31, 31, 31, 31, 31, 31},
+                                {31, 31, 31, 31, 31, 31},
+                                {31, 31, 31, 67, 31, 31},
+                                {31, 31, 31, 31, 31, 31},
+                                {31, 31, 31, 31, 31, 31}});
 
     FourierTransform ft;
-    std::vector<std::vector<double>> result = ft.ramplitude(signal);
+    double2d_t result = ft.ramplitude(signal);
     result = ft.fftshift(result);
 
     for (size_t i = 0; i < signal.size(); ++i)
@@ -145,24 +143,24 @@ TEST(FourierTransform, fft2DTest3)
 // 3x5 input with 1 row of all zeros
 TEST(FourierTransform, fft2DTest4)
 {
-    std::vector<std::vector<double>> signal({{1, 88, 0, 1, 0}, {0, 1, 1, 1, 0}, {0, 0, 0, 0, 0}});
+    double2d_t signal({{1, 88, 0, 1, 0}, {0, 1, 1, 1, 0}, {0, 0, 0, 0, 0}});
 
     FourierTransform ft;
 
     // shift
-    std::vector<std::vector<double>> signal_shifted = ft.fftshift(signal);
+    double2d_t signal_shifted = ft.fftshift(signal);
     signal_shifted = ft.ifftshift(signal_shifted);
     for (size_t i = 0; i < signal.size(); ++i)
         for (size_t j = 0; j < signal[0].size(); ++j)
             EXPECT_NEAR(signal_shifted[i][j], signal[i][j], 1e-13);
 
     // amplitude, shifted
-    std::vector<std::vector<double>> expected_result_amp(
+    double2d_t expected_result_amp(
         {{87.56947917, 85.92017286, 88.53812738, 88.59651195, 86.95382846},
          {88.02055461, 88.00785173, 93., 88.00785173, 88.02055461},
          {86.95382846, 88.59651195, 88.53812738, 85.92017286, 87.56947917}});
 
-    std::vector<std::vector<double>> result_amp = ft.ramplitude(signal);
+    double2d_t result_amp = ft.ramplitude(signal);
     result_amp = ft.fftshift(result_amp);
 
     for (size_t i = 0; i < signal.size(); ++i)
@@ -183,7 +181,7 @@ TEST(FourierTransform, fft2DTest4)
             EXPECT_NEAR(abs(result[i][j] - expected_result[i][j]), 0, abs(result[i][j]) * 1e-8);
 
     // inverse
-    std::vector<std::vector<double>> result_inv = ft.irfft(result, signal[0].size());
+    double2d_t result_inv = ft.irfft(result, signal[0].size());
 
     EXPECT_EQ(signal[0].size(), result_inv[0].size());
     for (size_t i = 0; i < signal.size(); ++i)
@@ -194,26 +192,24 @@ TEST(FourierTransform, fft2DTest4)
 // 4x4 input
 TEST(FourierTransform, fft2DTest5)
 {
-    std::vector<std::vector<double>> signal(
-        {{1, 0, 0, 5}, {1, 0, 0, 0}, {0, 1, 1, 1}, {0, 1, 1, 1}});
+    double2d_t signal({{1, 0, 0, 5}, {1, 0, 0, 0}, {0, 1, 1, 1}, {0, 1, 1, 1}});
 
     FourierTransform ft;
 
     // shift
-    std::vector<std::vector<double>> signal_shifted = ft.fftshift(signal);
+    double2d_t signal_shifted = ft.fftshift(signal);
     signal_shifted = ft.ifftshift(signal_shifted);
     for (size_t i = 0; i < signal.size(); ++i)
         for (size_t j = 0; j < signal[0].size(); ++j)
             EXPECT_NEAR(signal_shifted[i][j], signal[i][j], 1e-13);
 
     // amplitude, shifted
-    std::vector<std::vector<double>> expected_result_amp(
-        {{5., 5., 5., 5.},
-         {3.60555128, 3.60555128, 3.60555128, 7.28010989},
-         {5., 5., 13., 5.},
-         {3.60555128, 7.28010989, 3.60555128, 3.60555128}});
+    double2d_t expected_result_amp({{5., 5., 5., 5.},
+                                    {3.60555128, 3.60555128, 3.60555128, 7.28010989},
+                                    {5., 5., 13., 5.},
+                                    {3.60555128, 7.28010989, 3.60555128, 3.60555128}});
 
-    std::vector<std::vector<double>> result_amp = ft.ramplitude(signal);
+    double2d_t result_amp = ft.ramplitude(signal);
     result_amp = ft.fftshift(result_amp);
 
     for (size_t i = 0; i < signal.size(); ++i)
@@ -234,7 +230,7 @@ TEST(FourierTransform, fft2DTest5)
             EXPECT_NEAR(abs(result[i][j] - expected_result[i][j]), 0, abs(result[i][j]) * 1e-8);
 
     // inverse
-    std::vector<std::vector<double>> result_inv = ft.irfft(result, signal[0].size());
+    double2d_t result_inv = ft.irfft(result, signal[0].size());
 
     EXPECT_EQ(signal[0].size(), result_inv[0].size());
     for (size_t i = 0; i < signal.size(); ++i)
@@ -245,12 +241,12 @@ TEST(FourierTransform, fft2DTest5)
 // 7x7 input
 TEST(FourierTransform, fft2DTest6)
 {
-    std::vector<std::vector<double>> signal{
-        {1., 0., 0., 0., 0., 0., 0.}, {1., 0., 0., 0., 0., 0., 0.}, {1., 0., 5., 0., 0., 0., 0.},
-        {1., 0., 0., 0., 0., 0., 0.}, {1., 0., 0., 0., 0., 5., 0.}, {1., 0., 0., 0., 0., 0., 0.},
-        {1., 0., 0., 0., 0., 0., 0.}};
+    double2d_t signal{{1., 0., 0., 0., 0., 0., 0.}, {1., 0., 0., 0., 0., 0., 0.},
+                      {1., 0., 5., 0., 0., 0., 0.}, {1., 0., 0., 0., 0., 0., 0.},
+                      {1., 0., 0., 0., 0., 5., 0.}, {1., 0., 0., 0., 0., 0., 0.},
+                      {1., 0., 0., 0., 0., 0., 0.}};
 
-    std::vector<std::vector<double>> expected_result(
+    double2d_t expected_result(
         {{9.00968868, 6.23489802, 6.23489802, 9.00968868, 2.22520934, 10., 2.22520934},
          {9.00968868, 2.22520934, 10., 2.22520934, 9.00968868, 6.23489802, 6.23489802},
          {2.22520934, 9.00968868, 6.23489802, 6.23489802, 9.00968868, 2.22520934, 10.},
@@ -260,7 +256,7 @@ TEST(FourierTransform, fft2DTest6)
          {2.22520934, 10., 2.22520934, 9.00968868, 6.23489802, 6.23489802, 9.00968868}});
 
     FourierTransform ft;
-    std::vector<std::vector<double>> result = ft.ramplitude(signal);
+    double2d_t result = ft.ramplitude(signal);
     result = ft.fftshift(result);
 
     for (size_t i = 0; i < signal.size(); ++i)
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index e2a577d678854ad82a192cb2745216b1386fb4b4..cbc0df311880ae9cfcb20f637fffd2b3be76ca44 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -2080,7 +2080,7 @@ class Datafield(object):
         __init__(Datafield self, Frame frame) -> Datafield
         __init__(Datafield self, std::vector< Scale const *,std::allocator< Scale const * > > const & axes, vdouble1d_t values, vdouble1d_t errSigmas={}) -> Datafield
         __init__(Datafield self, std::vector< Scale const *,std::allocator< Scale const * > > const & axes) -> Datafield
-        __init__(Datafield self, std::string const & xlabel, std::string const & ylabel, vdouble2d_t vec) -> Datafield
+        __init__(Datafield self, std::string const & xlabel, std::string const & ylabel, double2d_t const & vec) -> Datafield
         __init__(Datafield self, Datafield arg2) -> Datafield
         """
         _libBornAgainDevice.Datafield_swiginit(self, _libBornAgainDevice.new_Datafield(*args))
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index a80054b35971f6da1696e60615b60f70db73cc3c..814b5c902668147fc61ee34180ef330a1a88cf2b 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -3682,8 +3682,8 @@ namespace Swig {
 #define SWIGTYPE_p_allocator_type swig_types[34]
 #define SWIGTYPE_p_char swig_types[35]
 #define SWIGTYPE_p_const_iterator swig_types[36]
-#define SWIGTYPE_p_corr_matrix_t swig_types[37]
-#define SWIGTYPE_p_difference_type swig_types[38]
+#define SWIGTYPE_p_difference_type swig_types[37]
+#define SWIGTYPE_p_double2d_t swig_types[38]
 #define SWIGTYPE_p_first_type swig_types[39]
 #define SWIGTYPE_p_int swig_types[40]
 #define SWIGTYPE_p_iterator swig_types[41]
@@ -29796,10 +29796,11 @@ SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_9(PyObject *self, Py_ssize_t nobj
   PyObject *resultobj = 0;
   std::string *arg1 = 0 ;
   std::string *arg2 = 0 ;
-  std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *arg3 = 0 ;
+  double2d_t *arg3 = 0 ;
   int res1 = SWIG_OLDOBJ ;
   int res2 = SWIG_OLDOBJ ;
-  int res3 = SWIG_OLDOBJ ;
+  void *argp3 = 0 ;
+  int res3 = 0 ;
   Datafield *result = 0 ;
   
   (void)self;
@@ -29826,20 +29827,17 @@ SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_9(PyObject *self, Py_ssize_t nobj
     }
     arg2 = ptr;
   }
-  {
-    std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *ptr = (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *)0;
-    res3 = swig::asptr(swig_obj[2], &ptr);
-    if (!SWIG_IsOK(res3)) {
-      SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "new_Datafield" "', argument " "3"" of type '" "std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_Datafield" "', argument " "3"" of type '" "std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &""'"); 
-    }
-    arg3 = ptr;
+  res3 = SWIG_ConvertPtr(swig_obj[2], &argp3, SWIGTYPE_p_double2d_t,  0  | 0);
+  if (!SWIG_IsOK(res3)) {
+    SWIG_exception_fail(SWIG_ArgError(res3), "in method '" "new_Datafield" "', argument " "3"" of type '" "double2d_t const &""'"); 
   }
+  if (!argp3) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "new_Datafield" "', argument " "3"" of type '" "double2d_t const &""'"); 
+  }
+  arg3 = reinterpret_cast< double2d_t * >(argp3);
   {
     try {
-      result = (Datafield *)new Datafield((std::string const &)*arg1,(std::string const &)*arg2,(std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &)*arg3);
+      result = (Datafield *)new Datafield((std::string const &)*arg1,(std::string const &)*arg2,(double2d_t const &)*arg3);
     } catch (const std::exception& ex) {
       // message shown in the Python interpreter
       const std::string msg {
@@ -29851,12 +29849,10 @@ SWIGINTERN PyObject *_wrap_new_Datafield__SWIG_9(PyObject *self, Py_ssize_t nobj
   resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Datafield, SWIG_POINTER_NEW |  0 );
   if (SWIG_IsNewObj(res1)) delete arg1;
   if (SWIG_IsNewObj(res2)) delete arg2;
-  if (SWIG_IsNewObj(res3)) delete arg3;
   return resultobj;
 fail:
   if (SWIG_IsNewObj(res1)) delete arg1;
   if (SWIG_IsNewObj(res2)) delete arg2;
-  if (SWIG_IsNewObj(res3)) delete arg3;
   return NULL;
 }
 
@@ -30025,7 +30021,7 @@ SWIGINTERN PyObject *_wrap_new_Datafield(PyObject *self, PyObject *args) {
       int res = SWIG_AsPtr_std_string(argv[1], (std::string**)(0));
       _v = SWIG_CheckState(res);
       if (_v) {
-        int res = swig::asptr(argv[2], (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > >**)(0));
+        int res = SWIG_ConvertPtr(argv[2], 0, SWIGTYPE_p_double2d_t, SWIG_POINTER_NO_NULL | 0);
         _v = SWIG_CheckState(res);
         if (_v) {
           return _wrap_new_Datafield__SWIG_9(self, argc, argv);
@@ -30067,7 +30063,7 @@ fail:
     "    Datafield::Datafield(std::vector< Scale const *,std::allocator< Scale const * > > const &,std::vector< double,std::allocator< double > > const &,std::vector< double,std::allocator< double > > const &)\n"
     "    Datafield::Datafield(std::vector< Scale const *,std::allocator< Scale const * > > const &,std::vector< double,std::allocator< double > > const &)\n"
     "    Datafield::Datafield(std::vector< Scale const *,std::allocator< Scale const * > > const &)\n"
-    "    Datafield::Datafield(std::string const &,std::string const &,std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > const &)\n"
+    "    Datafield::Datafield(std::string const &,std::string const &,double2d_t const &)\n"
     "    Datafield::Datafield(Datafield const &)\n");
   return 0;
 }
@@ -43807,7 +43803,7 @@ static PyMethodDef SwigMethods[] = {
 		"Datafield(Frame frame)\n"
 		"Datafield(std::vector< Scale const *,std::allocator< Scale const * > > const & axes, vdouble1d_t values, vdouble1d_t errSigmas={})\n"
 		"Datafield(std::vector< Scale const *,std::allocator< Scale const * > > const & axes)\n"
-		"Datafield(std::string const & xlabel, std::string const & ylabel, vdouble2d_t vec)\n"
+		"Datafield(std::string const & xlabel, std::string const & ylabel, double2d_t const & vec)\n"
 		"new_Datafield(Datafield arg1) -> Datafield\n"
 		""},
 	 { "delete_Datafield", _wrap_delete_Datafield, METH_O, "delete_Datafield(Datafield self)"},
@@ -44240,8 +44236,8 @@ static swig_type_info _swigt__p_VerticalLine = {"_p_VerticalLine", "VerticalLine
 static swig_type_info _swigt__p_allocator_type = {"_p_allocator_type", "allocator_type *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_char = {"_p_char", "char *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_const_iterator = {"_p_const_iterator", "const_iterator *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_corr_matrix_t = {"_p_corr_matrix_t", "corr_matrix_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_difference_type = {"_p_difference_type", "difference_type *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_double2d_t = {"_p_double2d_t", "double2d_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_first_type = {"_p_first_type", "first_type *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_int = {"_p_int", "int32_t *|int_fast16_t *|int_fast32_t *|int_least32_t *|intptr_t *|int *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_iterator = {"_p_iterator", "iterator *", 0, 0, (void*)0, 0};
@@ -44327,8 +44323,8 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_allocator_type,
   &_swigt__p_char,
   &_swigt__p_const_iterator,
-  &_swigt__p_corr_matrix_t,
   &_swigt__p_difference_type,
+  &_swigt__p_double2d_t,
   &_swigt__p_first_type,
   &_swigt__p_int,
   &_swigt__p_iterator,
@@ -44414,8 +44410,8 @@ static swig_cast_info _swigc__p_VerticalLine[] = {  {&_swigt__p_VerticalLine, 0,
 static swig_cast_info _swigc__p_allocator_type[] = {  {&_swigt__p_allocator_type, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_char[] = {  {&_swigt__p_char, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_const_iterator[] = {  {&_swigt__p_const_iterator, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_corr_matrix_t[] = {  {&_swigt__p_corr_matrix_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_difference_type[] = {  {&_swigt__p_difference_type, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_double2d_t[] = {  {&_swigt__p_double2d_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_first_type[] = {  {&_swigt__p_first_type, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_int[] = {  {&_swigt__p_int, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_iterator[] = {  {&_swigt__p_iterator, 0, 0, 0},{0, 0, 0, 0}};
@@ -44501,8 +44497,8 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_allocator_type,
   _swigc__p_char,
   _swigc__p_const_iterator,
-  _swigc__p_corr_matrix_t,
   _swigc__p_difference_type,
+  _swigc__p_double2d_t,
   _swigc__p_first_type,
   _swigc__p_int,
   _swigc__p_iterator,
diff --git a/auto/Wrap/libBornAgainFit.py b/auto/Wrap/libBornAgainFit.py
index a7d33684b2f5f2b0e605417c0731aa0b02908d7e..536c4ff40b05dfb076c6f45544a4294b4ead6009 100644
--- a/auto/Wrap/libBornAgainFit.py
+++ b/auto/Wrap/libBornAgainFit.py
@@ -1956,11 +1956,11 @@ class Parameters(object):
         return _libBornAgainFit.Parameters_setErrors(self, errors)
 
     def correlationMatrix(self):
-        r"""correlationMatrix(Parameters self) -> vdouble2d_t"""
+        r"""correlationMatrix(Parameters self) -> double2d_t"""
         return _libBornAgainFit.Parameters_correlationMatrix(self)
 
     def setCorrelationMatrix(self, matrix):
-        r"""setCorrelationMatrix(Parameters self, vdouble2d_t matrix)"""
+        r"""setCorrelationMatrix(Parameters self, double2d_t const & matrix)"""
         return _libBornAgainFit.Parameters_setCorrelationMatrix(self, matrix)
 
     def freeParameterCount(self):
diff --git a/auto/Wrap/libBornAgainFit_wrap.cpp b/auto/Wrap/libBornAgainFit_wrap.cpp
index 4fdf8885d54c5619ed6bbd53ccadfb43af127436..075b0b65dfaf85e45eb2a482f0bb14b84841f05e 100644
--- a/auto/Wrap/libBornAgainFit_wrap.cpp
+++ b/auto/Wrap/libBornAgainFit_wrap.cpp
@@ -3653,8 +3653,8 @@ namespace Swig {
 #define SWIGTYPE_p_allocator_type swig_types[5]
 #define SWIGTYPE_p_char swig_types[6]
 #define SWIGTYPE_p_const_iterator swig_types[7]
-#define SWIGTYPE_p_corr_matrix_t swig_types[8]
-#define SWIGTYPE_p_difference_type swig_types[9]
+#define SWIGTYPE_p_difference_type swig_types[8]
+#define SWIGTYPE_p_double2d_t swig_types[9]
 #define SWIGTYPE_p_fcn_residual_t swig_types[10]
 #define SWIGTYPE_p_fcn_scalar_t swig_types[11]
 #define SWIGTYPE_p_first_type swig_types[12]
@@ -27146,7 +27146,7 @@ SWIGINTERN PyObject *_wrap_Parameters_correlationMatrix(PyObject *self, PyObject
   void *argp1 = 0 ;
   int res1 = 0 ;
   PyObject *swig_obj[1] ;
-  mumufit::Parameters::corr_matrix_t result;
+  double2d_t result;
   
   (void)self;
   if (!args) SWIG_fail;
@@ -27167,7 +27167,7 @@ SWIGINTERN PyObject *_wrap_Parameters_correlationMatrix(PyObject *self, PyObject
       SWIG_exception(SWIG_RuntimeError, msg.c_str());
     }
   }
-  resultobj = swig::from(static_cast< std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > >(result));
+  resultobj = SWIG_NewPointerObj((new double2d_t(result)), SWIGTYPE_p_double2d_t, SWIG_POINTER_OWN |  0 );
   return resultobj;
 fail:
   return NULL;
@@ -27177,10 +27177,11 @@ fail:
 SWIGINTERN PyObject *_wrap_Parameters_setCorrelationMatrix(PyObject *self, PyObject *args) {
   PyObject *resultobj = 0;
   mumufit::Parameters *arg1 = (mumufit::Parameters *) 0 ;
-  mumufit::Parameters::corr_matrix_t *arg2 = 0 ;
+  double2d_t *arg2 = 0 ;
   void *argp1 = 0 ;
   int res1 = 0 ;
-  int res2 = SWIG_OLDOBJ ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
   PyObject *swig_obj[2] ;
   
   (void)self;
@@ -27190,20 +27191,17 @@ SWIGINTERN PyObject *_wrap_Parameters_setCorrelationMatrix(PyObject *self, PyObj
     SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "Parameters_setCorrelationMatrix" "', argument " "1"" of type '" "mumufit::Parameters *""'"); 
   }
   arg1 = reinterpret_cast< mumufit::Parameters * >(argp1);
-  {
-    std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *ptr = (std::vector< std::vector< double,std::allocator< double > >,std::allocator< std::vector< double,std::allocator< double > > > > *)0;
-    res2 = swig::asptr(swig_obj[1], &ptr);
-    if (!SWIG_IsOK(res2)) {
-      SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "Parameters_setCorrelationMatrix" "', argument " "2"" of type '" "mumufit::Parameters::corr_matrix_t const &""'"); 
-    }
-    if (!ptr) {
-      SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "Parameters_setCorrelationMatrix" "', argument " "2"" of type '" "mumufit::Parameters::corr_matrix_t const &""'"); 
-    }
-    arg2 = ptr;
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_double2d_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "Parameters_setCorrelationMatrix" "', argument " "2"" of type '" "double2d_t const &""'"); 
   }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "Parameters_setCorrelationMatrix" "', argument " "2"" of type '" "double2d_t const &""'"); 
+  }
+  arg2 = reinterpret_cast< double2d_t * >(argp2);
   {
     try {
-      (arg1)->setCorrelationMatrix((mumufit::Parameters::corr_matrix_t const &)*arg2);
+      (arg1)->setCorrelationMatrix((double2d_t const &)*arg2);
     } catch (const std::exception& ex) {
       // message shown in the Python interpreter
       const std::string msg {
@@ -27213,10 +27211,8 @@ SWIGINTERN PyObject *_wrap_Parameters_setCorrelationMatrix(PyObject *self, PyObj
     }
   }
   resultobj = SWIG_Py_Void();
-  if (SWIG_IsNewObj(res2)) delete arg2;
   return resultobj;
 fail:
-  if (SWIG_IsNewObj(res2)) delete arg2;
   return NULL;
 }
 
@@ -30054,8 +30050,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "Parameters_setValues", _wrap_Parameters_setValues, METH_VARARGS, "Parameters_setValues(Parameters self, vdouble1d_t values)"},
 	 { "Parameters_errors", _wrap_Parameters_errors, METH_O, "Parameters_errors(Parameters self) -> vdouble1d_t"},
 	 { "Parameters_setErrors", _wrap_Parameters_setErrors, METH_VARARGS, "Parameters_setErrors(Parameters self, vdouble1d_t errors)"},
-	 { "Parameters_correlationMatrix", _wrap_Parameters_correlationMatrix, METH_O, "Parameters_correlationMatrix(Parameters self) -> vdouble2d_t"},
-	 { "Parameters_setCorrelationMatrix", _wrap_Parameters_setCorrelationMatrix, METH_VARARGS, "Parameters_setCorrelationMatrix(Parameters self, vdouble2d_t matrix)"},
+	 { "Parameters_correlationMatrix", _wrap_Parameters_correlationMatrix, METH_O, "Parameters_correlationMatrix(Parameters self) -> double2d_t"},
+	 { "Parameters_setCorrelationMatrix", _wrap_Parameters_setCorrelationMatrix, METH_VARARGS, "Parameters_setCorrelationMatrix(Parameters self, double2d_t const & matrix)"},
 	 { "Parameters_freeParameterCount", _wrap_Parameters_freeParameterCount, METH_O, "Parameters_freeParameterCount(Parameters self) -> size_t"},
 	 { "Parameters___getitem__", _wrap_Parameters___getitem__, METH_VARARGS, "\n"
 		"Parameters___getitem__(Parameters self, std::string name) -> Parameter\n"
@@ -30127,8 +30123,8 @@ static swig_type_info _swigt__p_RealLimits = {"_p_RealLimits", "RealLimits *", 0
 static swig_type_info _swigt__p_allocator_type = {"_p_allocator_type", "allocator_type *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_char = {"_p_char", "char *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_const_iterator = {"_p_const_iterator", "const_iterator *", 0, 0, (void*)0, 0};
-static swig_type_info _swigt__p_corr_matrix_t = {"_p_corr_matrix_t", "corr_matrix_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_difference_type = {"_p_difference_type", "difference_type *", 0, 0, (void*)0, 0};
+static swig_type_info _swigt__p_double2d_t = {"_p_double2d_t", "double2d_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_fcn_residual_t = {"_p_fcn_residual_t", "fcn_residual_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_fcn_scalar_t = {"_p_fcn_scalar_t", "fcn_scalar_t *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_first_type = {"_p_first_type", "first_type *", 0, 0, (void*)0, 0};
@@ -30186,8 +30182,8 @@ static swig_type_info *swig_type_initial[] = {
   &_swigt__p_allocator_type,
   &_swigt__p_char,
   &_swigt__p_const_iterator,
-  &_swigt__p_corr_matrix_t,
   &_swigt__p_difference_type,
+  &_swigt__p_double2d_t,
   &_swigt__p_fcn_residual_t,
   &_swigt__p_fcn_scalar_t,
   &_swigt__p_first_type,
@@ -30245,8 +30241,8 @@ static swig_cast_info _swigc__p_RealLimits[] = {  {&_swigt__p_RealLimits, 0, 0,
 static swig_cast_info _swigc__p_allocator_type[] = {  {&_swigt__p_allocator_type, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_char[] = {  {&_swigt__p_char, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_const_iterator[] = {  {&_swigt__p_const_iterator, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_corr_matrix_t[] = {  {&_swigt__p_corr_matrix_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_difference_type[] = {  {&_swigt__p_difference_type, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_double2d_t[] = {  {&_swigt__p_double2d_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_fcn_residual_t[] = {  {&_swigt__p_fcn_residual_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_fcn_scalar_t[] = {  {&_swigt__p_fcn_scalar_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_first_type[] = {  {&_swigt__p_first_type, 0, 0, 0},{0, 0, 0, 0}};
@@ -30304,8 +30300,8 @@ static swig_cast_info *swig_cast_initial[] = {
   _swigc__p_allocator_type,
   _swigc__p_char,
   _swigc__p_const_iterator,
-  _swigc__p_corr_matrix_t,
   _swigc__p_difference_type,
+  _swigc__p_double2d_t,
   _swigc__p_fcn_residual_t,
   _swigc__p_fcn_scalar_t,
   _swigc__p_first_type,