diff --git a/Base/CMakeLists.txt b/Base/CMakeLists.txt
index 793b0599bc0f09b1bfa4b2c3eaddf9d70bbb8aa0..9cded5346edf0e0a09906729ce061d0b2fc925f7 100644
--- a/Base/CMakeLists.txt
+++ b/Base/CMakeLists.txt
@@ -22,6 +22,7 @@ target_link_libraries(${lib}
     PUBLIC
     ${GSL_LIBRARIES}
     ${FFTW3_LIBRARIES}
+    formfactor
     )
 target_include_directories(${lib}
     PUBLIC
diff --git a/Sample/CMakeLists.txt b/Sample/CMakeLists.txt
index 525d2b2b4ce325f645da3502ddf557bda0c4296b..d7cc1b00604095226fb77ead4f573fe116d1c4f5 100644
--- a/Sample/CMakeLists.txt
+++ b/Sample/CMakeLists.txt
@@ -21,8 +21,10 @@ set(${lib}_LIBRARY ${lib} PARENT_SCOPE)
 target_link_libraries(${lib}
     PUBLIC
     ${BornAgainParam_LIBRARY}
+    ${formfactor_LIBRARY}
     )
 target_include_directories(${lib}
     PUBLIC
+    ${formfactor_INCLUDE_DIR}
     ${CMAKE_SOURCE_DIR}
     )
diff --git a/Sample/HardParticle/FormFactorAnisoPyramid.cpp b/Sample/HardParticle/FormFactorAnisoPyramid.cpp
index 7bd3f6076645a533cfa2f587bd5ab734c3a1691b..0457d088fd79681e051bbc7f306e0f047d4c9e07 100644
--- a/Sample/HardParticle/FormFactorAnisoPyramid.cpp
+++ b/Sample/HardParticle/FormFactorAnisoPyramid.cpp
@@ -16,7 +16,7 @@
 #include "Base/Math/Constants.h"
 #include "Base/Math/Functions.h"
 
-const PolyhedralTopology FormFactorAnisoPyramid::topology = {{{{3, 2, 1, 0}, true},
+const ff::PolyhedralTopology FormFactorAnisoPyramid::topology = {{{{3, 2, 1, 0}, true},
                                                               {{0, 1, 5, 4}, false},
                                                               {{1, 2, 6, 5}, false},
                                                               {{2, 3, 7, 6}, false},
diff --git a/Sample/HardParticle/FormFactorAnisoPyramid.h b/Sample/HardParticle/FormFactorAnisoPyramid.h
index 97cce35eafd81b33bc36eca46391722b3c354449..554989cd8d70ad90697e65f69589a3a57f087b74 100644
--- a/Sample/HardParticle/FormFactorAnisoPyramid.h
+++ b/Sample/HardParticle/FormFactorAnisoPyramid.h
@@ -46,7 +46,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
 
     const double& m_length;
     const double& m_width;
diff --git a/Sample/HardParticle/FormFactorCantellatedCube.cpp b/Sample/HardParticle/FormFactorCantellatedCube.cpp
index 1bb00a8ddb81191da74b9385f0bf7e24b3286ed4..2de3e3bcd75a64db3c04ad626df3ff2e780cb1fd 100644
--- a/Sample/HardParticle/FormFactorCantellatedCube.cpp
+++ b/Sample/HardParticle/FormFactorCantellatedCube.cpp
@@ -14,7 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorCantellatedCube.h"
 
-const PolyhedralTopology FormFactorCantellatedCube::topology = {
+const ff::PolyhedralTopology FormFactorCantellatedCube::topology = {
     {
         /*  0 */ {{0, 1, 2, 3}, true},
         /*  1 */ {{0, 8, 5, 1}, true},
diff --git a/Sample/HardParticle/FormFactorCantellatedCube.h b/Sample/HardParticle/FormFactorCantellatedCube.h
index f36ba0ab430f01352a293fd0e7a72fd58d880be6..2c414125209ffd1e6d77fdf73468d8b4db3bafd5 100644
--- a/Sample/HardParticle/FormFactorCantellatedCube.h
+++ b/Sample/HardParticle/FormFactorCantellatedCube.h
@@ -41,7 +41,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_length;
     const double& m_removed_length;
 };
diff --git a/Sample/HardParticle/FormFactorCone6.cpp b/Sample/HardParticle/FormFactorCone6.cpp
index 88897f2eb39c768f507826c7faedac0fd144935d..ad382e2daee8de7c7e64cbd64ee4093ce23ae2f1 100644
--- a/Sample/HardParticle/FormFactorCone6.cpp
+++ b/Sample/HardParticle/FormFactorCone6.cpp
@@ -16,7 +16,7 @@
 #include "Base/Math/Constants.h"
 #include "Base/Math/Functions.h"
 
-const PolyhedralTopology FormFactorCone6::topology = {{{{5, 4, 3, 2, 1, 0}, true},
+const ff::PolyhedralTopology FormFactorCone6::topology = {{{{5, 4, 3, 2, 1, 0}, true},
                                                        {{0, 1, 7, 6}, false},
                                                        {{1, 2, 8, 7}, false},
                                                        {{2, 3, 9, 8}, false},
diff --git a/Sample/HardParticle/FormFactorCone6.h b/Sample/HardParticle/FormFactorCone6.h
index d07db56fe35a1fcb31752328454a2ff816eb2e77..e48c041198cc232e57dc810202315cff180d7514 100644
--- a/Sample/HardParticle/FormFactorCone6.h
+++ b/Sample/HardParticle/FormFactorCone6.h
@@ -45,7 +45,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_base_edge;
     const double& m_height;
     const double& m_alpha;
diff --git a/Sample/HardParticle/FormFactorCuboctahedron.cpp b/Sample/HardParticle/FormFactorCuboctahedron.cpp
index 690422731bc37a93ae71439ac1aae3ded2de45b4..930aa488a09340f62e4b80070bcbacbdf588df68 100644
--- a/Sample/HardParticle/FormFactorCuboctahedron.cpp
+++ b/Sample/HardParticle/FormFactorCuboctahedron.cpp
@@ -17,7 +17,7 @@
 #include "Base/Math/Functions.h"
 #include "Sample/HardParticle/FormFactorPyramid.h"
 
-const PolyhedralTopology FormFactorCuboctahedron::topology = {{{{3, 2, 1, 0}, true},
+const ff::PolyhedralTopology FormFactorCuboctahedron::topology = {{{{3, 2, 1, 0}, true},
                                                                {{0, 1, 5, 4}, false},
                                                                {{1, 2, 6, 5}, false},
                                                                {{2, 3, 7, 6}, false},
diff --git a/Sample/HardParticle/FormFactorCuboctahedron.h b/Sample/HardParticle/FormFactorCuboctahedron.h
index 3dc3e7e0be92410e2d4c90d24ceffbe2f563854c..c48308f92340d0f8aaa674de87e719183a2b428e 100644
--- a/Sample/HardParticle/FormFactorCuboctahedron.h
+++ b/Sample/HardParticle/FormFactorCuboctahedron.h
@@ -46,7 +46,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
 
     const double& m_length;
     const double& m_height;
diff --git a/Sample/HardParticle/FormFactorDodecahedron.cpp b/Sample/HardParticle/FormFactorDodecahedron.cpp
index 82a4e91d2368244a24b30ac0e2964681e7d2d43d..27140ae0791e670fe5d4c414e7ba6da2c7a04263 100644
--- a/Sample/HardParticle/FormFactorDodecahedron.cpp
+++ b/Sample/HardParticle/FormFactorDodecahedron.cpp
@@ -14,7 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorDodecahedron.h"
 
-const PolyhedralTopology FormFactorDodecahedron::topology = {{// bottom:
+const ff::PolyhedralTopology FormFactorDodecahedron::topology = {{// bottom:
                                                               {{0, 4, 3, 2, 1}, false},
                                                               // lower ring:
                                                               {{0, 5, 12, 9, 4}, false},
diff --git a/Sample/HardParticle/FormFactorDodecahedron.h b/Sample/HardParticle/FormFactorDodecahedron.h
index 8dc29c39157bf93b0808cb1595bf0c1a59122348..4b6e975dd375d87d7b1a383f36efcc3e3dbd253b 100644
--- a/Sample/HardParticle/FormFactorDodecahedron.h
+++ b/Sample/HardParticle/FormFactorDodecahedron.h
@@ -37,7 +37,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_edge;
 };
 
diff --git a/Sample/HardParticle/FormFactorIcosahedron.cpp b/Sample/HardParticle/FormFactorIcosahedron.cpp
index 1f0e426f8a3c44c3c7161b158eddd98d403b705f..5f840a9e80036ec2126470f8e759528ae3473bf5 100644
--- a/Sample/HardParticle/FormFactorIcosahedron.cpp
+++ b/Sample/HardParticle/FormFactorIcosahedron.cpp
@@ -14,7 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorIcosahedron.h"
 
-const PolyhedralTopology FormFactorIcosahedron::topology = {{// bottom:
+const ff::PolyhedralTopology FormFactorIcosahedron::topology = {{// bottom:
                                                              {{0, 2, 1}, false},
                                                              // 1st row:
                                                              {{0, 5, 2}, false},
diff --git a/Sample/HardParticle/FormFactorIcosahedron.h b/Sample/HardParticle/FormFactorIcosahedron.h
index 78757be0dc7b9713b815c0c65371cedc01eac104..c8663f4f2643c898b973b36859a4be65e953bbad 100644
--- a/Sample/HardParticle/FormFactorIcosahedron.h
+++ b/Sample/HardParticle/FormFactorIcosahedron.h
@@ -37,7 +37,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_edge;
 };
 
diff --git a/Sample/HardParticle/FormFactorPyramid.cpp b/Sample/HardParticle/FormFactorPyramid.cpp
index e29053643b3f75c2b621ba8afa5f6cf9e80511d7..32f07ad07136803f5aa11e1e947005be0194b12c 100644
--- a/Sample/HardParticle/FormFactorPyramid.cpp
+++ b/Sample/HardParticle/FormFactorPyramid.cpp
@@ -16,7 +16,7 @@
 #include "Base/Math/Constants.h"
 #include "Base/Math/Functions.h"
 
-const PolyhedralTopology FormFactorPyramid::topology = {{
+const ff::PolyhedralTopology FormFactorPyramid::topology = {{
                                                             {{3, 2, 1, 0}, true}, // TODO -> true
                                                             {{0, 1, 5, 4}, false},
                                                             {{1, 2, 6, 5}, false},
diff --git a/Sample/HardParticle/FormFactorPyramid.h b/Sample/HardParticle/FormFactorPyramid.h
index 0528093335d1e8fdaac4b989c14131d5ef256219..06c01e828db877fc8ffa5acb3e7d1b0a917d31ff 100644
--- a/Sample/HardParticle/FormFactorPyramid.h
+++ b/Sample/HardParticle/FormFactorPyramid.h
@@ -45,7 +45,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
 
     const double& m_base_edge;
     const double& m_height;
diff --git a/Sample/HardParticle/FormFactorTetrahedron.cpp b/Sample/HardParticle/FormFactorTetrahedron.cpp
index f56c380ab00947cd5efa8db95dc283c4364c452c..3728fde9cb7b9e6290b9cde35475f2a16e6696d9 100644
--- a/Sample/HardParticle/FormFactorTetrahedron.cpp
+++ b/Sample/HardParticle/FormFactorTetrahedron.cpp
@@ -16,7 +16,7 @@
 #include "Base/Math/Constants.h"
 #include "Base/Math/Functions.h"
 
-const PolyhedralTopology FormFactorTetrahedron::topology = {{{{2, 1, 0}, false},
+const ff::PolyhedralTopology FormFactorTetrahedron::topology = {{{{2, 1, 0}, false},
                                                              {{0, 1, 4, 3}, false},
                                                              {{1, 2, 5, 4}, false},
                                                              {{2, 0, 3, 5}, false},
diff --git a/Sample/HardParticle/FormFactorTetrahedron.h b/Sample/HardParticle/FormFactorTetrahedron.h
index 0535a8b3545f06c02623e2e982bee3359b6923d6..546cedc137991e3f7382634da3fd470c6e6867a9 100644
--- a/Sample/HardParticle/FormFactorTetrahedron.h
+++ b/Sample/HardParticle/FormFactorTetrahedron.h
@@ -45,7 +45,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_base_edge;
     const double& m_height;
     const double& m_alpha;
diff --git a/Sample/HardParticle/FormFactorTruncatedCube.cpp b/Sample/HardParticle/FormFactorTruncatedCube.cpp
index c50d198214b66d7d1bff7306f30b2dfc1ccd9cc3..6d16fcb9e82457b9436bc9271357e18fccddec56 100644
--- a/Sample/HardParticle/FormFactorTruncatedCube.cpp
+++ b/Sample/HardParticle/FormFactorTruncatedCube.cpp
@@ -14,7 +14,7 @@
 
 #include "Sample/HardParticle/FormFactorTruncatedCube.h"
 
-const PolyhedralTopology FormFactorTruncatedCube::topology = {
+const ff::PolyhedralTopology FormFactorTruncatedCube::topology = {
     {{{0, 1, 7, 6, 9, 10, 4, 3}, true},
      {{0, 2, 1}, false},
      {{3, 4, 5}, false},
diff --git a/Sample/HardParticle/FormFactorTruncatedCube.h b/Sample/HardParticle/FormFactorTruncatedCube.h
index 63ea81f73849eeafbd468f778b87d6b4b033e7ff..874e18d4d026ffbe2838e28afffc34de609b9a93 100644
--- a/Sample/HardParticle/FormFactorTruncatedCube.h
+++ b/Sample/HardParticle/FormFactorTruncatedCube.h
@@ -41,7 +41,7 @@ protected:
     void onChange() override;
 
 private:
-    static const PolyhedralTopology topology;
+    static const ff::PolyhedralTopology topology;
     const double& m_length;
     const double& m_removed_length;
 };
diff --git a/Sample/HardParticle/IFormFactorPolyhedron.cpp b/Sample/HardParticle/IFormFactorPolyhedron.cpp
index 4570920c9e0fc122e97065d74cabfa6f44a9ede2..d99a6663cd5d060ef94d93a3d7e43bc6b97c4a69 100644
--- a/Sample/HardParticle/IFormFactorPolyhedron.cpp
+++ b/Sample/HardParticle/IFormFactorPolyhedron.cpp
@@ -17,7 +17,8 @@
 //! "Form factor (Fourier shape transform) of polygon and polyhedron."
 
 #include "Sample/HardParticle/IFormFactorPolyhedron.h"
-#include "Sample/LibFF/Polyhedron.h"
+#include <ff/Polyhedron.h>
+#include <ff/PolyhedralComponents.h>
 
 // #ifdef ALGORITHM_DIAGNOSTIC // TODO restore
 // void IFormFactorPolyhedron::setLimits(double _q, int _n)
@@ -37,10 +38,10 @@ IFormFactorPolyhedron::~IFormFactorPolyhedron() = default;
 
 //! Called by child classes to set faces and other internal variables.
 
-void IFormFactorPolyhedron::setPolyhedron(const PolyhedralTopology& topology, double z_bottom,
+void IFormFactorPolyhedron::setPolyhedron(const ff::PolyhedralTopology& topology, double z_bottom,
                                           const std::vector<R3>& vertices)
 {
-    pimpl = std::make_unique<Polyhedron>(topology, z_bottom, vertices);
+    pimpl = std::make_unique<ff::Polyhedron>(topology, z_bottom, vertices);
 }
 
 double IFormFactorPolyhedron::bottomZ(const IRotation& rotation) const
diff --git a/Sample/HardParticle/IFormFactorPolyhedron.h b/Sample/HardParticle/IFormFactorPolyhedron.h
index d9d13814274122d3fbba194aafd92927f2d08777..41465b9a79f2abbb0516624a9e6ba561dd0f02fc 100644
--- a/Sample/HardParticle/IFormFactorPolyhedron.h
+++ b/Sample/HardParticle/IFormFactorPolyhedron.h
@@ -16,11 +16,13 @@
 #ifndef BORNAGAIN_SAMPLE_HARDPARTICLE_IFORMFACTORPOLYHEDRON_H
 #define BORNAGAIN_SAMPLE_HARDPARTICLE_IFORMFACTORPOLYHEDRON_H
 
-#include "Sample/LibFF/PolyhedralTopology.h"
+#include <ff/PolyhedralTopology.h>
 #include "Sample/Scattering/IBornFF.h"
 #include <memory>
 
+namespace ff {
 class Polyhedron;
+}
 
 //! A polyhedron, for form factor computation.
 
@@ -44,11 +46,11 @@ public:
     void assert_platonic() const;
 
 protected:
-    void setPolyhedron(const PolyhedralTopology& topology, double z_bottom,
+    void setPolyhedron(const ff::PolyhedralTopology& topology, double z_bottom,
                        const std::vector<R3>& vertices);
 
 private:
-    std::unique_ptr<Polyhedron> pimpl;
+    std::unique_ptr<ff::Polyhedron> pimpl;
 };
 
 #endif // BORNAGAIN_SAMPLE_HARDPARTICLE_IFORMFACTORPOLYHEDRON_H
diff --git a/Sample/HardParticle/Prism.cpp b/Sample/HardParticle/Prism.cpp
index 25f3354cbb241a06d5f3672a19cfa0e986a89d14..de27c9ce0bca5e88f5687e93633ad4bb00bcd411 100644
--- a/Sample/HardParticle/Prism.cpp
+++ b/Sample/HardParticle/Prism.cpp
@@ -30,7 +30,7 @@ Prism::Prism(bool symmetry_Ci, double height, const std::vector<R3>& vertices)
     }
 
     try {
-        m_base = std::make_unique<PolyhedralFace>(vertices, symmetry_Ci);
+        m_base = std::make_unique<ff::PolyhedralFace>(vertices, symmetry_Ci);
     } catch (std::invalid_argument& e) {
         throw std::invalid_argument(std::string("Invalid parameterization of Prism: ") + e.what());
     } catch (std::logic_error& e) {
diff --git a/Sample/HardParticle/Prism.h b/Sample/HardParticle/Prism.h
index 4d1ac331dd281f3a56f7b58d576f9232e74850bb..263e7ca272f6ae8b7257677c9d793e9f4408d0a2 100644
--- a/Sample/HardParticle/Prism.h
+++ b/Sample/HardParticle/Prism.h
@@ -20,8 +20,8 @@
 #ifndef BORNAGAIN_SAMPLE_HARDPARTICLE_PRISM_H
 #define BORNAGAIN_SAMPLE_HARDPARTICLE_PRISM_H
 
-#include "Sample/LibFF/PolyhedralComponents.h"
-#include "Sample/LibFF/PolyhedralTopology.h"
+#include <ff/PolyhedralComponents.h>
+#include <ff/PolyhedralTopology.h>
 #include <memory>
 
 class Prism {
@@ -34,7 +34,7 @@ public:
     complex_t evaluate_for_q(const C3& q) const;
 
 private:
-    std::unique_ptr<PolyhedralFace> m_base;
+    std::unique_ptr<ff::PolyhedralFace> m_base;
     double m_height;
     std::vector<R3> m_vertices; //! for topZ, bottomZ computation only
 };
diff --git a/Sample/LibFF/PolyhedralComponents.cpp b/Sample/LibFF/PolyhedralComponents.cpp
deleted file mode 100644
index f23c7addb610f4fd53b1eaf0ae025661672ba2d1..0000000000000000000000000000000000000000
--- a/Sample/LibFF/PolyhedralComponents.cpp
+++ /dev/null
@@ -1,377 +0,0 @@
-//  ************************************************************************************************
-//
-//  BornAgain: simulate and fit reflection and scattering
-//
-//! @file      Sample/LibFF/PolyhedralComponents.cpp
-//! @brief     Implements classes PolyhedralEdge, PolyhedralFace
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-//  ************************************************************************************************
-
-#include "Sample/LibFF/PolyhedralComponents.h"
-#include "Base/Math/Functions.h"
-#include "Base/Math/Precomputed.h"
-#include <iomanip>
-#include <stdexcept> // need overlooked by g++ 5.4
-
-namespace {
-const double eps = 2e-16;
-constexpr auto ReciprocalFactorialArray = Math::generateReciprocalFactorialArray<171>();
-} // namespace
-
-#ifdef ALGORITHM_DIAGNOSTIC
-void PolyhedralDiagnosis::reset()
-{
-    order = 0;
-    algo = 0;
-    msg.clear();
-};
-std::string PolyhedralDiagnosis::message() const
-{
-    std::string ret = "algo=" + std::to_string(algo) + ", order=" + std::to_string(order);
-    if (!msg.empty())
-        ret += ", msg:\n" + msg;
-    return ret;
-}
-bool PolyhedralDiagnosis::operator==(const PolyhedralDiagnosis& other) const
-{
-    return order == other.order && algo == other.algo;
-}
-bool PolyhedralDiagnosis::operator!=(const PolyhedralDiagnosis& other) const
-{
-    return !(*this == other);
-}
-#endif
-
-//  ************************************************************************************************
-//  PolyhedralEdge implementation
-//  ************************************************************************************************
-
-PolyhedralEdge::PolyhedralEdge(R3 Vlow, R3 Vhig) : m_E((Vhig - Vlow) / 2), m_R((Vhig + Vlow) / 2)
-{
-    if (m_E.mag2() == 0)
-        throw std::invalid_argument("At least one edge has zero length");
-};
-
-//! Returns sum_l=0^M/2 u^2l v^(M-2l) / (2l+1)!(M-2l)! - vperp^M/M!
-
-complex_t PolyhedralEdge::contrib(int M, C3 qpa, complex_t qrperp) const
-{
-    complex_t u = qE(qpa);
-    complex_t v2 = m_R.dot(qpa);
-    complex_t v1 = qrperp;
-    complex_t v = v2 + v1;
-    // std::cout << std::scientific << std::showpos << std::setprecision(16) << "contrib: u=" << u
-    //              << " v1=" << v1 << " v2=" << v2 << "\n";
-    if (v == 0.) { // only 2l=M contributes
-        if (M & 1) // M is odd
-            return 0.;
-        return ReciprocalFactorialArray[M] * (pow(u, M) / (M + 1.) - pow(v1, M));
-    }
-    complex_t ret = 0;
-    // the l=0 term, minus (qperp.R)^M, which cancels under the sum over E*contrib()
-    if (v1 == 0.) {
-        ret = ReciprocalFactorialArray[M] * pow(v2, M);
-    } else if (v2 == 0.) {
-        ; // leave ret=0
-    } else {
-        // binomial expansion
-        for (int mm = 1; mm <= M; ++mm) {
-            complex_t term = ReciprocalFactorialArray[mm] * ReciprocalFactorialArray[M - mm]
-                             * pow(v2, mm) * pow(v1, M - mm);
-            ret += term;
-            // std::cout << "contrib mm=" << mm << " t=" << term << " s=" << ret << "\n";
-        }
-    }
-    if (u == 0.)
-        return ret;
-    for (int l = 1; l <= M / 2; ++l) {
-        complex_t term = ReciprocalFactorialArray[M - 2 * l] * ReciprocalFactorialArray[2 * l + 1]
-                         * pow(u, 2 * l) * pow(v, M - 2 * l);
-        ret += term;
-        // std::cout << "contrib l=" << l << " t=" << term << " s=" << ret << "\n";
-    }
-    return ret;
-}
-
-//  ************************************************************************************************
-//  PolyhedralFace implementation
-//  ************************************************************************************************
-
-double PolyhedralFace::qpa_limit_series = 1e-2;
-int PolyhedralFace::n_limit_series = 20;
-
-//! Static method, returns diameter of circle that contains all vertices.
-
-double PolyhedralFace::diameter(const std::vector<R3>& V)
-{
-    double diameterFace = 0;
-    for (size_t j = 0; j < V.size(); ++j)
-        for (size_t jj = j + 1; jj < V.size(); ++jj)
-            diameterFace = std::max(diameterFace, (V[j] - V[jj]).mag());
-    return diameterFace;
-}
-
-//! Sets internal variables for given vertex chain.
-
-//! @param V oriented vertex list
-//! @param _sym_S2 true if face has a perpedicular two-fold symmetry axis
-
-PolyhedralFace::PolyhedralFace(const std::vector<R3>& V, bool _sym_S2) : sym_S2(_sym_S2)
-{
-    size_t NV = V.size();
-    if (!NV)
-        throw std::logic_error("Face with no edges");
-    if (NV < 3)
-        throw std::logic_error("Face with less than three edges");
-
-    // compute radius in 2d and 3d
-    m_radius_2d = diameter(V) / 2;
-    m_radius_3d = 0;
-    for (const R3& v : V)
-        m_radius_3d = std::max(m_radius_3d, v.mag());
-
-    // Initialize list of 'edges'.
-    // Do not create an edge if two vertices are too close to each other.
-    // TODO This is implemented in a somewhat sloppy way: we just skip an edge if it would
-    //      be too short. This leaves tiny open edges. In a clean implementation, we
-    //      rather should merge adjacent vertices before generating edges.
-    for (size_t j = 0; j < NV; ++j) {
-        size_t jj = (j + 1) % NV;
-        if ((V[j] - V[jj]).mag() < 1e-14 * m_radius_2d)
-            continue; // distance too short -> skip this edge
-        edges.emplace_back(V[j], V[jj]);
-    }
-    size_t NE = edges.size();
-    if (NE < 3)
-        throw std::invalid_argument("Face has less than three non-vanishing edges");
-
-    // compute n_k, rperp
-    m_normal = R3();
-    for (size_t j = 0; j < NE; ++j) {
-        size_t jj = (j + 1) % NE;
-        R3 ee = edges[j].E().cross(edges[jj].E());
-        if (ee.mag2() == 0) {
-            throw std::logic_error("Two adjacent edges are parallel");
-        }
-        m_normal += ee.unit();
-    }
-    m_normal /= NE;
-    m_rperp = 0;
-    for (size_t j = 0; j < NV; ++j)
-        m_rperp += V[j].dot(m_normal);
-    m_rperp /= NV;
-    // assert that the vertices lay in a plane
-    for (size_t j = 1; j < NV; ++j)
-        if (std::abs(V[j].dot(m_normal) - m_rperp) > 1e-14 * m_radius_3d)
-            throw std::logic_error("Face is not planar");
-    // compute m_area
-    m_area = 0;
-    for (size_t j = 0; j < NV; ++j) {
-        size_t jj = (j + 1) % NV;
-        m_area += m_normal.dot(V[j].cross(V[jj])) / 2;
-    }
-    // only now deal with inversion symmetry
-    if (sym_S2) {
-        if (NE & 1)
-            throw std::logic_error("Odd #edges violates symmetry S2");
-        NE /= 2;
-        for (size_t j = 0; j < NE; ++j) {
-            if (((edges[j].R() - m_rperp * m_normal) + (edges[j + NE].R() - m_rperp * m_normal))
-                    .mag()
-                > 1e-12 * m_radius_2d)
-                throw std::logic_error("Edge centers violate symmetry S2");
-            if ((edges[j].E() + edges[j + NE].E()).mag() > 1e-12 * m_radius_2d)
-                throw std::logic_error("Edge vectors violate symmetry S2");
-        }
-        // keep only half of the egdes
-        edges.erase(edges.begin() + NE, edges.end());
-    }
-}
-
-//! Sets qperp and qpa according to argument q and to this polygon's normal.
-
-void PolyhedralFace::decompose_q(C3 q, complex_t& qperp, C3& qpa) const
-{
-    qperp = m_normal.dot(q);
-    qpa = q - qperp * m_normal;
-    // improve numeric accuracy:
-    qpa -= m_normal.dot(qpa) * m_normal;
-    if (qpa.mag() < eps * std::abs(qperp))
-        qpa = C3(0., 0., 0.);
-}
-
-//! Returns core contribution to f_n
-
-complex_t PolyhedralFace::ff_n_core(int m, C3 qpa, complex_t qperp) const
-{
-    const C3 prevec = 2. * m_normal.cross(qpa); // complex conjugation not here but in .dot
-    complex_t ret = 0;
-    const complex_t qrperp = qperp * m_rperp;
-    for (size_t i = 0; i < edges.size(); ++i) {
-        const PolyhedralEdge& e = edges[i];
-        const complex_t vfac = prevec.dot(e.E());
-        const complex_t tmp = e.contrib(m + 1, qpa, qrperp);
-        ret += vfac * tmp;
-        //     std::cout << std::scientific << std::showpos << std::setprecision(16)
-        //               << "DBX ff_n_core " << m << " " << vfac << " " << tmp
-        //               << " term=" << vfac * tmp << " sum=" << ret << "\n";
-    }
-    return ret;
-}
-
-//! Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor.
-
-complex_t PolyhedralFace::ff_n(int n, C3 q) const
-{
-    complex_t qn = q.dot(m_normal); // conj(q)*normal (dot is antilinear in 'this' argument)
-    if (std::abs(qn) < eps * q.mag())
-        return 0.;
-    complex_t qperp;
-    C3 qpa;
-    decompose_q(q, qperp, qpa);
-    double qpa_mag2 = qpa.mag2();
-    if (qpa_mag2 == 0.)
-        return qn * pow(qperp * m_rperp, n) * m_area * ReciprocalFactorialArray[n];
-    if (sym_S2)
-        return qn * (ff_n_core(n, qpa, qperp) + ff_n_core(n, -qpa, qperp)) / qpa_mag2;
-    complex_t tmp = ff_n_core(n, qpa, qperp);
-    // std::cout << "DBX ff_n " << n << " " << qn << " " << tmp << " " << qpa_mag2 << "\n";
-    return qn * tmp / qpa_mag2;
-}
-
-//! Returns sum of n>=1 terms of qpa expansion of 2d form factor
-
-complex_t PolyhedralFace::expansion(complex_t fac_even, complex_t fac_odd, C3 qpa,
-                                    double abslevel) const
-{
-#ifdef ALGORITHM_DIAGNOSTIC
-    polyhedralDiagnosis.algo += 1;
-#endif
-    complex_t sum = 0;
-    complex_t n_fac = I;
-    int count_return_condition = 0;
-    for (int n = 1; n < n_limit_series; ++n) {
-#ifdef ALGORITHM_DIAGNOSTIC
-        polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
-#endif
-        complex_t term = n_fac * (n & 1 ? fac_odd : fac_even) * ff_n_core(n, qpa, 0) / qpa.mag2();
-        // std::cout << std::setprecision(16) << "    sum=" << sum << " +term=" << term << "\n";
-        sum += term;
-        if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * abslevel)
-            ++count_return_condition;
-        else
-            count_return_condition = 0;
-        if (count_return_condition > 2)
-            return sum; // regular exit
-        n_fac = mul_I(n_fac);
-    }
-    throw std::runtime_error("Series f(q_pa) not converged");
-}
-
-//! Returns core contribution to analytic 2d form factor.
-
-complex_t PolyhedralFace::edge_sum_ff(C3 q, C3 qpa, bool sym_Ci) const
-{
-    C3 prevec = m_normal.cross(qpa); // complex conjugation will take place in .dot
-    complex_t sum = 0;
-    complex_t vfacsum = 0;
-    for (size_t i = 0; i < edges.size(); ++i) {
-        const PolyhedralEdge& e = edges[i];
-        complex_t qE = e.qE(qpa);
-        complex_t qR = e.qR(qpa);
-        complex_t Rfac = sym_S2 ? sin(qR) : (sym_Ci ? cos(e.qR(q)) : exp_I(qR));
-        complex_t vfac;
-        if (sym_S2 || i < edges.size() - 1) {
-            vfac = prevec.dot(e.E());
-            vfacsum += vfac;
-        } else {
-            vfac = -vfacsum; // to improve numeric accuracy: qcE_J = - sum_{j=0}^{J-1} qcE_j
-        }
-        complex_t term = vfac * Math::sinc(qE) * Rfac;
-        sum += term;
-        //    std::cout << std::scientific << std::showpos << std::setprecision(16)
-        //              << "    sum=" << sum << " term=" << term << " vf=" << vfac << " qE=" << qE
-        //              << " qR=" << qR << " sinc=" << Math::sinc(qE) << " Rfac=" << Rfac << "\n";
-    }
-    return sum;
-}
-
-//! Returns the contribution ff(q) of this face to the polyhedral form factor.
-
-complex_t PolyhedralFace::ff(C3 q, bool sym_Ci) const
-{
-    complex_t qperp;
-    C3 qpa;
-    decompose_q(q, qperp, qpa);
-    double qpa_red = m_radius_2d * qpa.mag();
-    complex_t qr_perp = qperp * m_rperp;
-    complex_t ff0 = (sym_Ci ? 2. * I * sin(qr_perp) : exp_I(qr_perp)) * m_area;
-    if (qpa_red == 0)
-        return ff0;
-    if (qpa_red < qpa_limit_series && !sym_S2) {
-        // summation of power series
-        complex_t fac_even;
-        complex_t fac_odd;
-        if (sym_Ci) {
-            fac_even = 2. * mul_I(sin(qr_perp));
-            fac_odd = 2. * cos(qr_perp);
-        } else {
-            fac_even = exp_I(qr_perp);
-            fac_odd = fac_even;
-        }
-        return ff0 + expansion(fac_even, fac_odd, qpa, std::abs(ff0));
-    }
-    // direct evaluation of analytic formula
-    complex_t prefac;
-    if (sym_S2)
-        prefac = sym_Ci ? -8. * sin(qr_perp) : 4. * mul_I(exp_I(qr_perp));
-    else
-        prefac = sym_Ci ? 4. : 2. * exp_I(qr_perp);
-    // std::cout << "       qrperp=" << qr_perp << " => prefac=" << prefac << "\n";
-    return prefac * edge_sum_ff(q, qpa, sym_Ci) / mul_I(qpa.mag2());
-}
-
-//! Two-dimensional form factor, for use in prism, from power series.
-
-complex_t PolyhedralFace::ff_2D_expanded(C3 qpa) const
-{
-    return m_area + expansion(1., 1., qpa, std::abs(m_area));
-}
-
-//! Two-dimensional form factor, for use in prism, from sum over edge form factors.
-
-complex_t PolyhedralFace::ff_2D_direct(C3 qpa) const
-{
-    return (sym_S2 ? 4. : 2. / I) * edge_sum_ff(qpa, qpa, false) / qpa.mag2();
-}
-
-//! Returns the two-dimensional form factor of this face, for use in a prism.
-
-complex_t PolyhedralFace::ff_2D(C3 qpa) const
-{
-    if (std::abs(qpa.dot(m_normal)) > eps * qpa.mag())
-        throw std::logic_error("ff_2D called with perpendicular q component");
-    double qpa_red = m_radius_2d * qpa.mag();
-    if (qpa_red == 0)
-        return m_area;
-    if (qpa_red < qpa_limit_series && !sym_S2)
-        return ff_2D_expanded(qpa);
-    return ff_2D_direct(qpa);
-}
-
-//! Throws if deviation from inversion symmetry is detected. Does not check vertices.
-
-void PolyhedralFace::assert_Ci(const PolyhedralFace& other) const
-{
-    if (std::abs(m_rperp - other.m_rperp) > 1e-15 * (m_rperp + other.m_rperp))
-        throw std::logic_error("Faces with different distance from origin violate symmetry Ci");
-    if (std::abs(m_area - other.m_area) > 1e-15 * (m_area + other.m_area))
-        throw std::logic_error("Faces with different areas violate symmetry Ci");
-    if ((m_normal + other.m_normal).mag() > 1e-14)
-        throw std::logic_error("Faces do not have opposite orientation, violating symmetry Ci");
-}
diff --git a/Sample/LibFF/PolyhedralComponents.h b/Sample/LibFF/PolyhedralComponents.h
deleted file mode 100644
index d645b2103fcfed58962f1813934863960f4c700b..0000000000000000000000000000000000000000
--- a/Sample/LibFF/PolyhedralComponents.h
+++ /dev/null
@@ -1,98 +0,0 @@
-//  ************************************************************************************************
-//
-//  BornAgain: simulate and fit reflection and scattering
-//
-//! @file      Sample/LibFF/PolyhedralComponents.h
-//! @brief     Defines classes PolyhedralEdge, PolyhedralFace
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-//  ************************************************************************************************
-
-#ifdef SWIG
-#error no need to expose this header to Swig
-#endif
-
-#ifndef USER_API
-#ifndef BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALCOMPONENTS_H
-#define BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALCOMPONENTS_H
-
-#include <heinz/Complex.h>
-#include <heinz/Vectors3D.h>
-#include <vector>
-
-#ifdef ALGORITHM_DIAGNOSTIC
-#include <string>
-struct PolyhedralDiagnosis {
-    int algo;
-    int order;
-    std::string msg;
-    void reset();
-    std::string message() const;
-    bool operator==(const PolyhedralDiagnosis&) const;
-    bool operator!=(const PolyhedralDiagnosis&) const;
-};
-inline PolyhedralDiagnosis polyhedralDiagnosis;
-#endif
-
-//! One edge of a polygon, for form factor computation.
-
-class PolyhedralEdge {
-public:
-    PolyhedralEdge(R3 Vlow, R3 Vhig);
-
-    R3 E() const { return m_E; }
-    R3 R() const { return m_R; }
-    complex_t qE(C3 q) const { return m_E.dot(q); }
-    complex_t qR(C3 q) const { return m_R.dot(q); }
-
-    complex_t contrib(int m, C3 qpa, complex_t qrperp) const;
-
-private:
-    R3 m_E; //!< vector pointing from mid of edge to upper vertex
-    R3 m_R; //!< position vector of edge midpoint
-};
-
-//! A polygon, for form factor computation.
-
-class PolyhedralFace {
-public:
-    static double diameter(const std::vector<R3>& V);
-
-    PolyhedralFace(const std::vector<R3>& _V = std::vector<R3>(), bool _sym_S2 = false);
-
-    double area() const { return m_area; }
-    double pyramidalVolume() const { return m_rperp * m_area / 3; }
-    double radius3d() const { return m_radius_3d; }
-    //! Returns conj(q)*normal [Vec3::dot is antilinear in 'this' argument]
-    complex_t normalProjectionConj(C3 q) const { return q.dot(m_normal); }
-    complex_t ff_n(int n, C3 q) const;
-    complex_t ff(C3 q, bool sym_Ci) const;
-    complex_t ff_2D(C3 qpa) const;
-    complex_t ff_2D_direct(C3 qpa) const;   // for TestTriangle
-    complex_t ff_2D_expanded(C3 qpa) const; // for TestTriangle
-    void assert_Ci(const PolyhedralFace& other) const;
-
-private:
-    static double qpa_limit_series; //!< determines when use power series
-    static int n_limit_series;
-
-    bool sym_S2; //!< if true, then edges obtainable by inversion are not provided
-    std::vector<PolyhedralEdge> edges;
-    double m_area;
-    R3 m_normal;        //!< normal vector of this polygon's plane
-    double m_rperp;     //!< distance of this polygon's plane from the origin, along 'm_normal'
-    double m_radius_2d; //!< radius of enclosing cylinder
-    double m_radius_3d; //!< radius of enclosing sphere
-
-    void decompose_q(C3 q, complex_t& qperp, C3& qpa) const;
-    complex_t ff_n_core(int m, C3 qpa, complex_t qperp) const;
-    complex_t edge_sum_ff(C3 q, C3 qpa, bool sym_Ci) const;
-    complex_t expansion(complex_t fac_even, complex_t fac_odd, C3 qpa, double abslevel) const;
-};
-
-#endif // BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALCOMPONENTS_H
-#endif // USER_API
diff --git a/Sample/LibFF/PolyhedralTopology.h b/Sample/LibFF/PolyhedralTopology.h
deleted file mode 100644
index 933cae137e6773d7642c9d1035453084c186719c..0000000000000000000000000000000000000000
--- a/Sample/LibFF/PolyhedralTopology.h
+++ /dev/null
@@ -1,40 +0,0 @@
-//  ************************************************************************************************
-//
-//  BornAgain: simulate and fit reflection and scattering
-//
-//! @file      Sample/LibFF/PolyhedralTopology.h
-//! @brief     Defines classes PolygonalTopology, PolyhedralTopology
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-//  ************************************************************************************************
-
-#ifdef SWIG
-#error no need to expose this header to Swig
-#endif
-
-#ifndef USER_API
-#ifndef BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALTOPOLOGY_H
-#define BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALTOPOLOGY_H
-
-#include <vector>
-
-//! For internal use in PolyhedralFace.
-class PolygonalTopology {
-public:
-    std::vector<int> vertexIndices;
-    bool symmetry_S2;
-};
-
-//! For internal use in IFormFactorPolyhedron.
-class PolyhedralTopology {
-public:
-    std::vector<PolygonalTopology> faces;
-    bool symmetry_Ci;
-};
-
-#endif // BORNAGAIN_SAMPLE_LIBFF_POLYHEDRALTOPOLOGY_H
-#endif // USER_API
diff --git a/Sample/LibFF/Polyhedron.cpp b/Sample/LibFF/Polyhedron.cpp
deleted file mode 100644
index 15ebbd824e8935f425c7cb3e791bfb9d818a3cac..0000000000000000000000000000000000000000
--- a/Sample/LibFF/Polyhedron.cpp
+++ /dev/null
@@ -1,209 +0,0 @@
-//  ************************************************************************************************
-//
-//  BornAgain: simulate and fit reflection and scattering
-//
-//! @file      Sample/LibFF/Polyhedron.cpp
-//! @brief     Implements class Polyhedron.
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-//  ************************************************************************************************
-
-//! The mathematics implemented here is described in full detail in a paper
-//! by Joachim Wuttke, entitled
-//! "Numerically stable form factor of any polygon and polyhedron"
-
-#include "Sample/LibFF/Polyhedron.h"
-#include "Base/Math/Functions.h"
-#include "Sample/LibFF/PolyhedralTopology.h"
-#include <boost/format.hpp>
-#include <iomanip>
-#include <iostream>
-#include <stdexcept> // need overlooked by g++ 5.4
-
-namespace {
-const double eps = 2e-16;
-const double q_limit_series = 1e-2;
-const int n_limit_series = 20;
-} // namespace
-
-Polyhedron::Polyhedron(const PolyhedralTopology& topology, double z_bottom,
-                       const std::vector<R3>& vertices)
-{
-
-    m_vertices.clear();
-    for (const R3& vertex : vertices)
-        m_vertices.push_back(vertex - R3{0, 0, z_bottom});
-
-    try {
-        m_z_bottom = z_bottom;
-        m_sym_Ci = topology.symmetry_Ci;
-
-        double diameter = 0;
-        for (size_t j = 0; j < vertices.size(); ++j)
-            for (size_t jj = j + 1; jj < vertices.size(); ++jj)
-                diameter = std::max(diameter, (vertices[j] - vertices[jj]).mag());
-
-        m_faces.clear();
-        for (const PolygonalTopology& tf : topology.faces) {
-            std::vector<R3> corners; // of one face
-            for (int i : tf.vertexIndices)
-                corners.push_back(vertices[i]);
-            if (PolyhedralFace::diameter(corners) <= 1e-14 * diameter)
-                continue; // skip ridiculously small face
-            m_faces.emplace_back(corners, tf.symmetry_S2);
-        }
-        if (m_faces.size() < 4)
-            throw std::logic_error("Less than four non-vanishing faces");
-
-        m_radius = 0;
-        m_volume = 0;
-        for (const PolyhedralFace& Gk : m_faces) {
-            m_radius = std::max(m_radius, Gk.radius3d());
-            m_volume += Gk.pyramidalVolume();
-        }
-        if (m_sym_Ci) {
-            if (m_faces.size() & 1)
-                throw std::logic_error("Odd #faces violates symmetry Ci");
-            size_t N = m_faces.size() / 2;
-            // for this tests, m_faces must be in a specific order
-            for (size_t k = 0; k < N; ++k)
-                m_faces[k].assert_Ci(m_faces[2 * N - 1 - k]);
-            // keep only half of the faces
-            m_faces.erase(m_faces.begin() + N, m_faces.end());
-        }
-    } catch (std::invalid_argument& e) {
-        throw std::invalid_argument(std::string("Invalid parameterization of Polyhedron: ")
-                                    + e.what());
-    } catch (std::logic_error& e) {
-        throw std::logic_error(std::string("Bug in Polyhedron: ") + e.what()
-                               + " [please report to the maintainers]");
-    } catch (std::exception& e) {
-        throw std::runtime_error(std::string("Unexpected exception in Polyhedron: ") + e.what()
-                                 + " [please report to the maintainers]");
-    }
-}
-
-Polyhedron::~Polyhedron() = default;
-
-void Polyhedron::assert_platonic() const
-{
-    // just one test; one could do much more ...
-    double pyramidal_volume = 0;
-    for (const auto& Gk : m_faces)
-        pyramidal_volume += Gk.pyramidalVolume();
-    pyramidal_volume /= m_faces.size();
-    for (const auto& Gk : m_faces)
-        if (std::abs(Gk.pyramidalVolume() - pyramidal_volume) > 160 * eps * pyramidal_volume) {
-            std::cerr << std::setprecision(16)
-                      << "Bug: pyr_volume(this face)=" << Gk.pyramidalVolume()
-                      << " vs pyr_volume(avge)=" << pyramidal_volume << "\n";
-            throw std::runtime_error("Deviant pyramidal volume in Platonic Polyhedron");
-        }
-}
-
-double Polyhedron::volume() const
-{
-    return m_volume;
-}
-double Polyhedron::radius() const
-{
-    return m_radius;
-}
-
-const std::vector<R3>& Polyhedron::vertices() const
-{
-    return m_vertices;
-}
-
-//! Returns the form factor F(q) of this polyhedron, respecting the offset z_bottom.
-
-complex_t Polyhedron::evaluate_for_q(const C3& q) const
-{
-    try {
-        return exp_I(-m_z_bottom * q.z()) * evaluate_centered(q);
-    } catch (std::logic_error& e) {
-        throw std::logic_error(std::string("Bug in Polyhedron: ") + e.what()
-                               + " [please report to the maintainers]");
-    } catch (std::runtime_error& e) {
-        throw std::runtime_error(std::string("Numeric computation failed in Polyhedron: ")
-                                 + e.what() + " [please report to the maintainers]");
-    } catch (std::exception& e) {
-        throw std::runtime_error(std::string("Unexpected exception in Polyhedron: ") + e.what()
-                                 + " [please report to the maintainers]");
-    }
-}
-
-//! Returns the form factor F(q) of this polyhedron, with origin at z=0.
-
-complex_t Polyhedron::evaluate_centered(const C3& q) const
-{
-    double q_red = m_radius * q.mag();
-#ifdef ALGORITHM_DIAGNOSTIC
-    polyhedralDiagnosis.reset();
-#endif
-    if (q_red == 0)
-        return m_volume;
-    if (q_red < q_limit_series) {
-        // summation of power series
-#ifdef ALGORITHM_DIAGNOSTIC
-        polyhedralDiagnosis.algo = 100;
-#endif
-        complex_t sum = 0;
-        complex_t n_fac = (m_sym_Ci ? -2 : -1) / q.mag2();
-        int count_return_condition = 0;
-        for (int n = 2; n < n_limit_series; ++n) {
-            if (m_sym_Ci && n & 1)
-                continue;
-#ifdef ALGORITHM_DIAGNOSTIC
-            polyhedralDiagnosis.order = std::max(polyhedralDiagnosis.order, n);
-#endif
-            complex_t term = 0;
-            for (const PolyhedralFace& Gk : m_faces) {
-                complex_t tmp = Gk.ff_n(n + 1, q);
-                term += tmp;
-            }
-            term *= n_fac;
-#ifdef ALGORITHM_DIAGNOSTIC_LEVEL2
-            polyhedralDiagnosis.msg +=
-                boost::str(boost::format("  + term(n=%2i) = %23.17e+i*%23.17e\n") % n % term.real()
-                           % term.imag());
-#endif
-            sum += term;
-            if (std::abs(term) <= eps * std::abs(sum) || std::abs(sum) < eps * m_volume)
-                ++count_return_condition;
-            else
-                count_return_condition = 0;
-            if (count_return_condition > 2)
-                return m_volume + sum; // regular exit
-            n_fac = m_sym_Ci ? -n_fac : mul_I(n_fac);
-        }
-        throw std::runtime_error("Series F(q) not converged");
-    }
-
-    // direct evaluation of analytic formula (coefficients may involve series)
-#ifdef ALGORITHM_DIAGNOSTIC
-    polyhedralDiagnosis.algo = 200;
-#endif
-    complex_t sum = 0;
-    for (const PolyhedralFace& Gk : m_faces) {
-        complex_t qn = Gk.normalProjectionConj(q); // conj(q)*normal
-        if (std::abs(qn) < eps * q.mag())
-            continue;
-        complex_t term = qn * Gk.ff(q, m_sym_Ci);
-#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
-        polyhedralDiagnosis.msg += boost::str(boost::format("  + face_ff = %23.17e+i*%23.17e\n")
-                                              % term.real() % term.imag());
-#endif
-        sum += term;
-    }
-#ifdef ALGORITHM_DIAGNOSTIC //_LEVEL2
-    polyhedralDiagnosis.msg +=
-        boost::str(boost::format(" -> raw sum = %23.17e+i*%23.17e; divisor = %23.17e\n")
-                   % sum.real() % sum.imag() % q.mag2());
-#endif
-    return sum / I / q.mag2();
-}
diff --git a/Sample/LibFF/Polyhedron.h b/Sample/LibFF/Polyhedron.h
deleted file mode 100644
index b603a930c7a9d66cc46d9ddd96c93a7a3d97921a..0000000000000000000000000000000000000000
--- a/Sample/LibFF/Polyhedron.h
+++ /dev/null
@@ -1,56 +0,0 @@
-//  ************************************************************************************************
-//
-//  BornAgain: simulate and fit reflection and scattering
-//
-//! @file      Sample/LibFF/Polyhedron.h
-//! @brief     Defines class Polyhedron.
-//!
-//! @homepage  http://www.bornagainproject.org
-//! @license   GNU General Public License v3 or higher (see COPYING)
-//! @copyright Forschungszentrum Jülich GmbH 2018
-//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
-//
-//  ************************************************************************************************
-
-#ifdef SWIG
-#error no need to expose this header to Swig
-#endif
-
-#ifndef USER_API
-#ifndef BORNAGAIN_SAMPLE_LIBFF_POLYHEDRON_H
-#define BORNAGAIN_SAMPLE_LIBFF_POLYHEDRON_H
-
-#include "Sample/LibFF/PolyhedralComponents.h"
-#include <memory>
-
-class PolyhedralTopology;
-
-//! A polyhedron, implementation class for use in IFormFactorPolyhedron
-
-class Polyhedron {
-public:
-    Polyhedron(const PolyhedralTopology& topology, double z_bottom,
-               const std::vector<R3>& vertices);
-    Polyhedron(const Polyhedron&) = delete;
-    ~Polyhedron();
-
-    void assert_platonic() const;
-    double volume() const;
-    double radius() const;
-
-    const std::vector<R3>& vertices() const; //! needed for topZ, bottomZ computation
-    complex_t evaluate_for_q(const C3& q) const;
-    complex_t evaluate_centered(const C3& q) const;
-
-private:
-    double m_z_bottom;
-    bool m_sym_Ci; //!< if true, then faces obtainable by inversion are not provided
-
-    std::vector<PolyhedralFace> m_faces;
-    double m_radius;
-    double m_volume;
-    std::vector<R3> m_vertices; //! for topZ, bottomZ computation only
-};
-
-#endif // BORNAGAIN_SAMPLE_LIBFF_POLYHEDRON_H
-#endif // USER_API
diff --git a/Tests/SingleUse/ff-tetrahedron.cpp b/Tests/SingleUse/ff-tetrahedron.cpp
index 94f467860d271dbfc0192f020e51ca3e1dabab1b..a9d6b0a46f2ab020c9cc37136844fc361e09dd47 100644
--- a/Tests/SingleUse/ff-tetrahedron.cpp
+++ b/Tests/SingleUse/ff-tetrahedron.cpp
@@ -1,6 +1,6 @@
 #include "Base/Math/Constants.h"
 #include "Sample/HardParticle/HardParticles.h"
-#include "Sample/LibFF/PolyhedralComponents.h"
+<ff/PolyhedralComponents.h>
 #include "Tests/GTestWrapper/google_test.h"
 #include <iostream>
 
diff --git a/Tests/SingleUse/ff-triangle.cpp b/Tests/SingleUse/ff-triangle.cpp
deleted file mode 100644
index 06ebb61bc7299c8731bcde33fa8ad99634c2d9ac..0000000000000000000000000000000000000000
--- a/Tests/SingleUse/ff-triangle.cpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#include "Base/Math/Constants.h"
-#include "Sample/LibFF/PolyhedralComponents.h"
-#include "Tests/GTestWrapper/google_test.h"
-#include <cstdio>
-
-//! Ad-hoc test of triangle form factor.
-//!
-//! Used while preparing polyhedral form factor manuscript for publication - JWu, dec 2020.
-//!
-//! To reactivate this code, just move in to ../UnitTest/Numeric/.
-
-class FFTriangleTest : public ::testing::Test {
-};
-
-TEST_F(FFTriangleTest, Triangle)
-{
-    const double a = 1.;
-    const double as = a / 2;
-    const double ac = a / sqrt(3) / 2;
-    const double ah = a / sqrt(3);
-    const std::vector<R3> V{{-ac, as, 0.}, {-ac, -as, 0.}, {ah, 0., 0.}};
-
-    const PolyhedralFace T(V, false);
-    EXPECT_NEAR(T.area(), sqrt(3) / 4, 1e-15);
-
-    int failures = 0;
-    const int M = 37;
-    for (int j = 0; j < M; ++j) {
-        const double phi = M_PI_2 * j / (M - 1);
-        const C3 uQ{sin(phi), cos(phi), 0.};
-        const int N = 2800 + j;
-        for (int i = 0; i < N; ++i) {
-            const double q = 1e-17 * pow(1.7e17, i / (N - 1.));
-            const C3 Q = uQ * q;
-            const double f1 = std::abs(T.ff_2D_direct(Q));
-            const double f2 = std::abs(T.ff_2D_expanded(Q));
-            const double relerr = std::abs(f1 - f2) / f2;
-            if (relerr > 7e-16) {
-                printf("ERR1 %9.6f %16.11e %21.16e %21.16e %10.4e\n", phi, q, f1, f2, relerr);
-                ++failures;
-            }
-            if (q > 1e-7)
-                continue;
-            const double relerr2 = std::abs(f1 - T.area()) / f2;
-            if (relerr2 > 7e-16) {
-                printf("ERR2 %9.6f %16.11e %21.16e %21.16e %10.4e\n", phi, q, f1, f2, relerr2);
-                ++failures;
-            }
-        }
-    }
-    EXPECT_EQ(failures, 0);
-}
diff --git a/Tests/Unit/Numeric/BisectFF.cpp b/Tests/Unit/Numeric/BisectFF.cpp
index 01781b883cee79badd5ff377a99090318cd54978..521eda3e4a99bbe57ca7f43373541532fb1a858f 100644
--- a/Tests/Unit/Numeric/BisectFF.cpp
+++ b/Tests/Unit/Numeric/BisectFF.cpp
@@ -2,7 +2,7 @@
 
 #include "Base/Math/Constants.h"
 #include "Sample/HardParticle/HardParticles.h"
-#include "Sample/LibFF/PolyhedralComponents.h" // for diagnostic
+<ff/PolyhedralComponents.h> // for diagnostic
 #include "Tests/GTestWrapper/google_test.h"
 #include <cassert>
 #include <complex>
diff --git a/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp b/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
index 6516b1ab44c60a1ea094f6d8a3bfde4813969f75..a17fba8eef1c393b5cdd9df1f3ca5d4970757c6a 100644
--- a/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
+++ b/Tests/Unit/Numeric/FormFactorSpecializationTest.cpp
@@ -1,5 +1,5 @@
 #include "Sample/HardParticle/HardParticles.h"
-#include "Sample/LibFF/PolyhedralComponents.h" // for diagnostic
+#include <ff/PolyhedralComponents.h> // for diagnostic
 #include "Tests/GTestWrapper/google_test.h"
 #include "Tests/Unit/Numeric/MultiQTestbed.h"
 
diff --git a/Tests/Unit/Numeric/FormFactorSymmetryTest.cpp b/Tests/Unit/Numeric/FormFactorSymmetryTest.cpp
index a9d4b679563f5510e96a31c9c5c6b8d7fb09c406..0245f534a9d673297c0bb9d88a39ff276a32cdca 100644
--- a/Tests/Unit/Numeric/FormFactorSymmetryTest.cpp
+++ b/Tests/Unit/Numeric/FormFactorSymmetryTest.cpp
@@ -1,6 +1,6 @@
 #include "Base/Math/Constants.h"
 #include "Sample/HardParticle/HardParticles.h"
-#include "Sample/LibFF/PolyhedralComponents.h" // for diagnostic
+#include <ff/PolyhedralComponents.h> // for diagnostic
 #include "Tests/GTestWrapper/google_test.h"
 #include "Tests/Unit/Numeric/MultiQTestbed.h"
 
diff --git a/auto/Doc/man/bornagain.1 b/auto/Doc/man/bornagain.1
index 3a43708bef531db48ce991af8c1a9a5bb66b0865..1abbe2bfc4582f0cc481eb032118e34b8b5eff1e 100644
--- a/auto/Doc/man/bornagain.1
+++ b/auto/Doc/man/bornagain.1
@@ -133,7 +133,7 @@
 .\" ========================================================================
 .\"
 .IX Title "BORNAGAIN 1"
-.TH BORNAGAIN 1 "2021-03-09" "perl v5.32.1" "BornAgain manual"
+.TH BORNAGAIN 1 "2020-11-23" "perl v5.32.1" "BornAgain manual"
 .\" For nroff, turn off justification.  Always turn off hyphenation; it makes
 .\" way too many mistakes in technical documents.
 .if n .ad l
diff --git a/auto/Wrap/doxygenSample.i b/auto/Wrap/doxygenSample.i
index 796fc91041af1c7ff78a32ae41093271db31195a..a604e8a6599d30573fd008ddbf94d784ef8015f9 100644
--- a/auto/Wrap/doxygenSample.i
+++ b/auto/Wrap/doxygenSample.i
@@ -5093,162 +5093,6 @@ Sets the relative weight of this layout.
 ";
 
 
-// File: classPolygonalTopology.xml
-%feature("docstring") PolygonalTopology "
-
-For internal use in  PolyhedralFace.
-
-C++ includes: PolyhedralTopology.h
-";
-
-
-// File: classPolyhedralEdge.xml
-%feature("docstring") PolyhedralEdge "
-
-One edge of a polygon, for form factor computation.
-
-C++ includes: PolyhedralComponents.h
-";
-
-%feature("docstring")  PolyhedralEdge::PolyhedralEdge "PolyhedralEdge::PolyhedralEdge(R3 Vlow, R3 Vhig)
-";
-
-%feature("docstring")  PolyhedralEdge::E "R3 PolyhedralEdge::E() const
-";
-
-%feature("docstring")  PolyhedralEdge::R "R3 PolyhedralEdge::R() const
-";
-
-%feature("docstring")  PolyhedralEdge::qE "complex_t PolyhedralEdge::qE(C3 q) const
-";
-
-%feature("docstring")  PolyhedralEdge::qR "complex_t PolyhedralEdge::qR(C3 q) const
-";
-
-%feature("docstring")  PolyhedralEdge::contrib "complex_t PolyhedralEdge::contrib(int m, C3 qpa, complex_t qrperp) const
-
-Returns sum_l=0^M/2 u^2l v^(M-2l) / (2l+1)!(M-2l)! - vperp^M/M! 
-";
-
-
-// File: classPolyhedralFace.xml
-%feature("docstring") PolyhedralFace "
-
-A polygon, for form factor computation.
-
-C++ includes: PolyhedralComponents.h
-";
-
-%feature("docstring")  PolyhedralFace::PolyhedralFace "PolyhedralFace::PolyhedralFace(const std::vector< R3 > &_V=std::vector< R3 >(), bool _sym_S2=false)
-
-Sets internal variables for given vertex chain.
-
-Parameters:
------------
-
-V: 
-oriented vertex list
-
-_sym_S2: 
-true if face has a perpedicular two-fold symmetry axis 
-";
-
-%feature("docstring")  PolyhedralFace::area "double PolyhedralFace::area() const
-";
-
-%feature("docstring")  PolyhedralFace::pyramidalVolume "double PolyhedralFace::pyramidalVolume() const
-";
-
-%feature("docstring")  PolyhedralFace::radius3d "double PolyhedralFace::radius3d() const
-";
-
-%feature("docstring")  PolyhedralFace::normalProjectionConj "complex_t PolyhedralFace::normalProjectionConj(C3 q) const
-
-Returns conj(q)*normal [Vec3::dot is antilinear in 'this' argument]. 
-";
-
-%feature("docstring")  PolyhedralFace::ff_n "complex_t PolyhedralFace::ff_n(int n, C3 q) const
-
-Returns contribution qn*f_n [of order q^(n+1)] from this face to the polyhedral form factor. 
-";
-
-%feature("docstring")  PolyhedralFace::ff "complex_t PolyhedralFace::ff(C3 q, bool sym_Ci) const
-
-Returns the contribution ff(q) of this face to the polyhedral form factor. 
-";
-
-%feature("docstring")  PolyhedralFace::ff_2D "complex_t PolyhedralFace::ff_2D(C3 qpa) const
-
-Returns the two-dimensional form factor of this face, for use in a prism. 
-";
-
-%feature("docstring")  PolyhedralFace::ff_2D_direct "complex_t PolyhedralFace::ff_2D_direct(C3 qpa) const
-
-Two-dimensional form factor, for use in prism, from sum over edge form factors. 
-";
-
-%feature("docstring")  PolyhedralFace::ff_2D_expanded "complex_t PolyhedralFace::ff_2D_expanded(C3 qpa) const
-
-Two-dimensional form factor, for use in prism, from power series. 
-";
-
-%feature("docstring")  PolyhedralFace::assert_Ci "void PolyhedralFace::assert_Ci(const PolyhedralFace &other) const
-
-Throws if deviation from inversion symmetry is detected. Does not check vertices. 
-";
-
-
-// File: classPolyhedralTopology.xml
-%feature("docstring") PolyhedralTopology "
-
-For internal use in  IFormFactorPolyhedron.
-
-C++ includes: PolyhedralTopology.h
-";
-
-
-// File: classPolyhedron.xml
-%feature("docstring") Polyhedron "
-
-A polyhedron, implementation class for use in  IFormFactorPolyhedron.
-
-C++ includes: Polyhedron.h
-";
-
-%feature("docstring")  Polyhedron::Polyhedron "Polyhedron::Polyhedron(const PolyhedralTopology &topology, double z_bottom, const std::vector< R3 > &vertices)
-";
-
-%feature("docstring")  Polyhedron::Polyhedron "Polyhedron::Polyhedron(const Polyhedron &)=delete
-";
-
-%feature("docstring")  Polyhedron::~Polyhedron "Polyhedron::~Polyhedron()
-";
-
-%feature("docstring")  Polyhedron::assert_platonic "void Polyhedron::assert_platonic() const
-";
-
-%feature("docstring")  Polyhedron::volume "double Polyhedron::volume() const
-";
-
-%feature("docstring")  Polyhedron::radius "double Polyhedron::radius() const
-";
-
-%feature("docstring")  Polyhedron::vertices "const std::vector< R3 > & Polyhedron::vertices() const
-";
-
-%feature("docstring")  Polyhedron::evaluate_for_q "complex_t Polyhedron::evaluate_for_q(const C3 &q) const
-
-needed for topZ, bottomZ computation
-
-Returns the form factor F(q) of this polyhedron, respecting the offset z_bottom. 
-";
-
-%feature("docstring")  Polyhedron::evaluate_centered "complex_t Polyhedron::evaluate_centered(const C3 &q) const
-
-Returns the form factor F(q) of this polyhedron, with origin at z=0. 
-";
-
-
 // File: classPrism.xml
 %feature("docstring") Prism "";
 
@@ -5622,49 +5466,43 @@ C++ includes: ZLimits.h
 ";
 
 
-// File: namespace_0d120.xml
-
-
-// File: namespace_0d123.xml
-
+// File: namespace_0d124.xml
 
-// File: namespace_0d129.xml
 
+// File: namespace_0d128.xml
 
-// File: namespace_0d133.xml
 
+// File: namespace_0d132.xml
 
-// File: namespace_0d137.xml
 
+// File: namespace_0d142.xml
 
-// File: namespace_0d147.xml
 
-
-// File: namespace_0d149.xml
+// File: namespace_0d144.xml
 
 
 // File: namespace_0d16.xml
 
 
-// File: namespace_0d180.xml
+// File: namespace_0d175.xml
 
 
 // File: namespace_0d2.xml
 
 
-// File: namespace_0d208.xml
+// File: namespace_0d203.xml
 
 
-// File: namespace_0d221.xml
+// File: namespace_0d216.xml
 
 
-// File: namespace_0d231.xml
+// File: namespace_0d226.xml
 
 
-// File: namespace_0d253.xml
+// File: namespace_0d248.xml
 
 
-// File: namespace_0d266.xml
+// File: namespace_0d261.xml
 
 
 // File: namespace_0d36.xml
@@ -5926,6 +5764,9 @@ Returns a body-centered cubic (cI) lattice with edge length a. TODO: Clarify mea
 ";
 
 
+// File: namespaceff.xml
+
+
 // File: namespaceMaterialUtils.xml
 %feature("docstring")  MaterialUtils::ScalarReducedPotential "complex_t MaterialUtils::ScalarReducedPotential(complex_t n, R3 k, double n_ref)
 
@@ -6406,21 +6247,6 @@ Used by the hard sphere and by several soft sphere classes.
 // File: Lattice3D_8h.xml
 
 
-// File: PolyhedralComponents_8cpp.xml
-
-
-// File: PolyhedralComponents_8h.xml
-
-
-// File: PolyhedralTopology_8h.xml
-
-
-// File: Polyhedron_8cpp.xml
-
-
-// File: Polyhedron_8h.xml
-
-
 // File: SomeFormFactors_8cpp.xml
 
 
diff --git a/auto/Wrap/libBornAgainBase.py b/auto/Wrap/libBornAgainBase.py
index 5bcc8ab05ebc70b5f876497d7b344deebce0f33f..789add81281aeedb162be1ac13df00fe87de2d37 100644
--- a/auto/Wrap/libBornAgainBase.py
+++ b/auto/Wrap/libBornAgainBase.py
@@ -2982,7 +2982,7 @@ class FixedBinAxis(IAxis):
 # Register FixedBinAxis in _libBornAgainBase:
 _libBornAgainBase.FixedBinAxis_swigregister(FixedBinAxis)
 
-class R3(arrayR3_t):
+class R3(object):
     r"""Proxy of C++ Vec3< double > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -2990,8 +2990,8 @@ class R3(arrayR3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(R3 self, double const x_, double const y_, double const z_) -> R3
         __init__(R3 self) -> R3
-        __init__(R3 self, double const x, double const y, double const z) -> R3
         """
         _libBornAgainBase.R3_swiginit(self, _libBornAgainBase.new_R3(*args))
 
@@ -3075,6 +3075,14 @@ class R3(arrayR3_t):
         r"""real(R3 self) -> R3"""
         return _libBornAgainBase.R3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(R3 self, R3 other) -> bool"""
+        return _libBornAgainBase.R3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(R3 self, R3 other) -> bool"""
+        return _libBornAgainBase.R3___ne__(self, other)
+
     def angle(self, v):
         r"""angle(R3 self, R3 v) -> double"""
         return _libBornAgainBase.R3_angle(self, v)
@@ -3254,7 +3262,7 @@ class vector_R3(object):
 # Register vector_R3 in _libBornAgainBase:
 _libBornAgainBase.vector_R3_swigregister(vector_R3)
 
-class C3(arrayC3_t):
+class C3(object):
     r"""Proxy of C++ Vec3< std::complex< double > > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -3262,8 +3270,8 @@ class C3(arrayC3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(C3 self, std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_) -> C3
         __init__(C3 self) -> C3
-        __init__(C3 self, std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3
         """
         _libBornAgainBase.C3_swiginit(self, _libBornAgainBase.new_C3(*args))
 
@@ -3327,6 +3335,14 @@ class C3(arrayC3_t):
         r"""real(C3 self) -> R3"""
         return _libBornAgainBase.C3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(C3 self, C3 other) -> bool"""
+        return _libBornAgainBase.C3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(C3 self, C3 other) -> bool"""
+        return _libBornAgainBase.C3___ne__(self, other)
+
     def project(self, v):
         r"""project(C3 self, C3 v) -> C3"""
         return _libBornAgainBase.C3_project(self, v)
diff --git a/auto/Wrap/libBornAgainBase_wrap.cpp b/auto/Wrap/libBornAgainBase_wrap.cpp
index 99d2de6b3d62e31d2d81da0047e5aceeddddc483..b4d1acfa18494306daa7c063b114c0e432e0ff0a 100644
--- a/auto/Wrap/libBornAgainBase_wrap.cpp
+++ b/auto/Wrap/libBornAgainBase_wrap.cpp
@@ -30306,20 +30306,7 @@ SWIGINTERN PyObject *FixedBinAxis_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObj
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< double > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< double > *)new Vec3< double >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   double arg1 ;
   double arg2 ;
@@ -30356,6 +30343,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< double > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< double > *)new Vec3< double >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -30365,7 +30365,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_R3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_R3__SWIG_0(self, argc, argv);
+    return _wrap_new_R3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -30384,7 +30384,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_R3__SWIG_1(self, argc, argv);
+          return _wrap_new_R3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -30393,8 +30393,8 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_R3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< double >::Vec3()\n"
-    "    Vec3< double >::Vec3(double const,double const,double const)\n");
+    "    Vec3< double >::Vec3(double const,double const,double const)\n"
+    "    Vec3< double >::Vec3()\n");
   return 0;
 }
 
@@ -30903,6 +30903,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_R3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___eq__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator ==((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___ne__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator !=((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_R3_angle(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
@@ -32824,20 +32894,7 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< std::complex< double > > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::complex< double > arg1 ;
   std::complex< double > arg2 ;
@@ -32874,6 +32931,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -32883,7 +32953,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_C3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_C3__SWIG_0(self, argc, argv);
+    return _wrap_new_C3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -32902,7 +32972,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_C3__SWIG_1(self, argc, argv);
+          return _wrap_new_C3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -32911,8 +32981,8 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_C3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< std::complex< double > >::Vec3()\n"
-    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n");
+    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n"
+    "    Vec3< std::complex< double > >::Vec3()\n");
   return 0;
 }
 
@@ -33306,6 +33376,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_C3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___eq__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator ==((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___ne__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator !=((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_C3_project(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
@@ -36488,8 +36628,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "FixedBinAxis_swigregister", FixedBinAxis_swigregister, METH_O, NULL},
 	 { "FixedBinAxis_swiginit", FixedBinAxis_swiginit, METH_VARARGS, NULL},
 	 { "new_R3", _wrap_new_R3, METH_VARARGS, "\n"
-		"R3()\n"
-		"new_R3(double const x, double const y, double const z) -> R3\n"
+		"R3(double const x_, double const y_, double const z_)\n"
+		"new_R3() -> R3\n"
 		""},
 	 { "R3_x", _wrap_R3_x, METH_O, "R3_x(R3 self) -> double"},
 	 { "R3_y", _wrap_R3_y, METH_O, "R3_y(R3 self) -> double"},
@@ -36511,6 +36651,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
+	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
+	 { "R3___ne__", _wrap_R3___ne__, METH_VARARGS, "R3___ne__(R3 self, R3 other) -> bool"},
 	 { "R3_angle", _wrap_R3_angle, METH_VARARGS, "R3_angle(R3 self, R3 v) -> double"},
 	 { "R3_project", _wrap_R3_project, METH_VARARGS, "R3_project(R3 self, R3 v) -> R3"},
 	 { "delete_R3", _wrap_delete_R3, METH_O, "delete_R3(R3 self)"},
@@ -36579,8 +36721,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
 	 { "new_C3", _wrap_new_C3, METH_VARARGS, "\n"
-		"C3()\n"
-		"new_C3(std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3\n"
+		"C3(std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_)\n"
+		"new_C3() -> C3\n"
 		""},
 	 { "C3_x", _wrap_C3_x, METH_O, "C3_x(C3 self) -> std::complex< double >"},
 	 { "C3_y", _wrap_C3_y, METH_O, "C3_y(C3 self) -> std::complex< double >"},
@@ -36597,6 +36739,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
 	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
+	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
+	 { "C3___ne__", _wrap_C3___ne__, METH_VARARGS, "C3___ne__(C3 self, C3 other) -> bool"},
 	 { "C3_project", _wrap_C3_project, METH_VARARGS, "C3_project(C3 self, C3 v) -> C3"},
 	 { "delete_C3", _wrap_delete_C3, METH_O, "delete_C3(C3 self)"},
 	 { "C3_swigregister", C3_swigregister, METH_O, NULL},
@@ -36691,12 +36835,6 @@ static void *_p_CustomBinAxisTo_p_IAxis(void *x, int *SWIGUNUSEDPARM(newmemory))
 static void *_p_FixedBinAxisTo_p_IAxis(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((IAxis *)  ((FixedBinAxis *) x));
 }
-static void *_p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< std::complex< double >,3 > *)  ((Vec3< std::complex< double > > *) x));
-}
-static void *_p_Vec3T_double_tTo_p_std__arrayT_double_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< double,3 > *)  ((Vec3< double > *) x));
-}
 static swig_type_info _swigt__p_Bin1D = {"_p_Bin1D", "Bin1D *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_Bin1DCVector = {"_p_Bin1DCVector", "Bin1DCVector *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_Bin1DKVector = {"_p_Bin1DKVector", "Bin1DKVector *", 0, 0, (void*)0, 0};
@@ -36869,8 +37007,8 @@ static swig_cast_info _swigc__p_std__allocatorT_std__string_t[] = {  {&_swigt__p
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_unsigned_long_t[] = {  {&_swigt__p_std__allocatorT_unsigned_long_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_Vec3T_double_t, _p_Vec3T_double_tTo_p_std__arrayT_double_3_t, 0, 0},  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},  {&_swigt__p_Vec3T_std__complexT_double_t_t, _p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__complexT_double_t[] = {  {&_swigt__p_std__complexT_double_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__invalid_argument[] = {  {&_swigt__p_std__invalid_argument, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__lessT_std__string_t[] = {  {&_swigt__p_std__lessT_std__string_t, 0, 0, 0},{0, 0, 0, 0}};
diff --git a/auto/Wrap/libBornAgainCore.py b/auto/Wrap/libBornAgainCore.py
index 12c2ffae9285b80faa325289d3edbf77d5159ed6..c0d072b71fe22963690b0fa1d4b7231559f22712 100644
--- a/auto/Wrap/libBornAgainCore.py
+++ b/auto/Wrap/libBornAgainCore.py
@@ -1904,7 +1904,7 @@ _libBornAgainCore.vector_pvacuum_double_t_swigregister(vector_pvacuum_double_t)
 
 import libBornAgainFit
 import libBornAgainBase
-class R3(arrayR3_t):
+class R3(object):
     r"""Proxy of C++ Vec3< double > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -1912,8 +1912,8 @@ class R3(arrayR3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(R3 self, double const x_, double const y_, double const z_) -> R3
         __init__(R3 self) -> R3
-        __init__(R3 self, double const x, double const y, double const z) -> R3
         """
         _libBornAgainCore.R3_swiginit(self, _libBornAgainCore.new_R3(*args))
 
@@ -1997,6 +1997,14 @@ class R3(arrayR3_t):
         r"""real(R3 self) -> R3"""
         return _libBornAgainCore.R3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(R3 self, R3 other) -> bool"""
+        return _libBornAgainCore.R3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(R3 self, R3 other) -> bool"""
+        return _libBornAgainCore.R3___ne__(self, other)
+
     def angle(self, v):
         r"""angle(R3 self, R3 v) -> double"""
         return _libBornAgainCore.R3_angle(self, v)
@@ -2192,7 +2200,7 @@ class vector_R3(object):
 # Register vector_R3 in _libBornAgainCore:
 _libBornAgainCore.vector_R3_swigregister(vector_R3)
 
-class C3(arrayC3_t):
+class C3(object):
     r"""Proxy of C++ Vec3< std::complex< double > > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -2200,8 +2208,8 @@ class C3(arrayC3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(C3 self, std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_) -> C3
         __init__(C3 self) -> C3
-        __init__(C3 self, std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3
         """
         _libBornAgainCore.C3_swiginit(self, _libBornAgainCore.new_C3(*args))
 
@@ -2265,6 +2273,14 @@ class C3(arrayC3_t):
         r"""real(C3 self) -> R3"""
         return _libBornAgainCore.C3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(C3 self, C3 other) -> bool"""
+        return _libBornAgainCore.C3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(C3 self, C3 other) -> bool"""
+        return _libBornAgainCore.C3___ne__(self, other)
+
     def project(self, v):
         r"""project(C3 self, C3 v) -> C3"""
         return _libBornAgainCore.C3_project(self, v)
diff --git a/auto/Wrap/libBornAgainCore_wrap.cpp b/auto/Wrap/libBornAgainCore_wrap.cpp
index cb4071560ed2f929268f29be68c8d05a1c40d907..e5b0a49ec000e239af4de7dd3d0b8b619dc84122 100644
--- a/auto/Wrap/libBornAgainCore_wrap.cpp
+++ b/auto/Wrap/libBornAgainCore_wrap.cpp
@@ -26941,20 +26941,7 @@ SWIGINTERN PyObject *vector_pvacuum_double_t_swiginit(PyObject *SWIGUNUSEDPARM(s
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< double > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< double > *)new Vec3< double >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   double arg1 ;
   double arg2 ;
@@ -26991,6 +26978,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< double > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< double > *)new Vec3< double >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -27000,7 +27000,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_R3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_R3__SWIG_0(self, argc, argv);
+    return _wrap_new_R3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -27019,7 +27019,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_R3__SWIG_1(self, argc, argv);
+          return _wrap_new_R3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -27028,8 +27028,8 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_R3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< double >::Vec3()\n"
-    "    Vec3< double >::Vec3(double const,double const,double const)\n");
+    "    Vec3< double >::Vec3(double const,double const,double const)\n"
+    "    Vec3< double >::Vec3()\n");
   return 0;
 }
 
@@ -27538,6 +27538,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_R3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___eq__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator ==((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___ne__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator !=((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_R3_angle(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
@@ -29583,20 +29653,7 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< std::complex< double > > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::complex< double > arg1 ;
   std::complex< double > arg2 ;
@@ -29633,6 +29690,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -29642,7 +29712,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_C3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_C3__SWIG_0(self, argc, argv);
+    return _wrap_new_C3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -29661,7 +29731,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_C3__SWIG_1(self, argc, argv);
+          return _wrap_new_C3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -29670,8 +29740,8 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_C3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< std::complex< double > >::Vec3()\n"
-    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n");
+    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n"
+    "    Vec3< std::complex< double > >::Vec3()\n");
   return 0;
 }
 
@@ -30065,6 +30135,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_C3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___eq__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator ==((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___ne__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator !=((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_C3_project(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
@@ -46397,8 +46537,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_pvacuum_double_t_swigregister", vector_pvacuum_double_t_swigregister, METH_O, NULL},
 	 { "vector_pvacuum_double_t_swiginit", vector_pvacuum_double_t_swiginit, METH_VARARGS, NULL},
 	 { "new_R3", _wrap_new_R3, METH_VARARGS, "\n"
-		"R3()\n"
-		"new_R3(double const x, double const y, double const z) -> R3\n"
+		"R3(double const x_, double const y_, double const z_)\n"
+		"new_R3() -> R3\n"
 		""},
 	 { "R3_x", _wrap_R3_x, METH_O, "R3_x(R3 self) -> double"},
 	 { "R3_y", _wrap_R3_y, METH_O, "R3_y(R3 self) -> double"},
@@ -46420,6 +46560,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
+	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
+	 { "R3___ne__", _wrap_R3___ne__, METH_VARARGS, "R3___ne__(R3 self, R3 other) -> bool"},
 	 { "R3_angle", _wrap_R3_angle, METH_VARARGS, "R3_angle(R3 self, R3 v) -> double"},
 	 { "R3_project", _wrap_R3_project, METH_VARARGS, "R3_project(R3 self, R3 v) -> R3"},
 	 { "R3___add__", _wrap_R3___add__, METH_VARARGS, "R3___add__(R3 self, R3 rhs) -> R3"},
@@ -46492,8 +46634,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
 	 { "new_C3", _wrap_new_C3, METH_VARARGS, "\n"
-		"C3()\n"
-		"new_C3(std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3\n"
+		"C3(std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_)\n"
+		"new_C3() -> C3\n"
 		""},
 	 { "C3_x", _wrap_C3_x, METH_O, "C3_x(C3 self) -> std::complex< double >"},
 	 { "C3_y", _wrap_C3_y, METH_O, "C3_y(C3 self) -> std::complex< double >"},
@@ -46510,6 +46652,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
 	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
+	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
+	 { "C3___ne__", _wrap_C3___ne__, METH_VARARGS, "C3___ne__(C3 self, C3 other) -> bool"},
 	 { "C3_project", _wrap_C3_project, METH_VARARGS, "C3_project(C3 self, C3 v) -> C3"},
 	 { "delete_C3", _wrap_delete_C3, METH_O, "delete_C3(C3 self)"},
 	 { "C3_swigregister", C3_swigregister, METH_O, NULL},
@@ -47930,12 +48074,6 @@ static void *_p_ISampleNodeTo_p_INode(void *x, int *SWIGUNUSEDPARM(newmemory)) {
 static void *_p_IBornFFTo_p_INode(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((INode *) (ISampleNode *)(IFormFactor *) ((IBornFF *) x));
 }
-static void *_p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< std::complex< double >,3 > *)  ((Vec3< std::complex< double > > *) x));
-}
-static void *_p_Vec3T_double_tTo_p_std__arrayT_double_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< double,3 > *)  ((Vec3< double > *) x));
-}
 static void *_p_GISASSimulationTo_p_ISimulation2D(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((ISimulation2D *)  ((GISASSimulation *) x));
 }
@@ -48272,8 +48410,8 @@ static swig_cast_info _swigc__p_std__allocatorT_std__string_t[] = {  {&_swigt__p
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_unsigned_long_t[] = {  {&_swigt__p_std__allocatorT_unsigned_long_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_Vec3T_double_t, _p_Vec3T_double_tTo_p_std__arrayT_double_3_t, 0, 0},  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},  {&_swigt__p_Vec3T_std__complexT_double_t_t, _p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__complexT_double_t[] = {  {&_swigt__p_std__complexT_double_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__functionT_bool_fsize_tF_t[] = {  {&_swigt__p_std__functionT_bool_fsize_tF_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__invalid_argument[] = {  {&_swigt__p_std__invalid_argument, 0, 0, 0},{0, 0, 0, 0}};
diff --git a/auto/Wrap/libBornAgainDevice.py b/auto/Wrap/libBornAgainDevice.py
index d5d2e91d9b0f6c879848ff0d496d58d0f4b6c4ae..550059af8d20c51790a38699a4c364e3812887cd 100644
--- a/auto/Wrap/libBornAgainDevice.py
+++ b/auto/Wrap/libBornAgainDevice.py
@@ -1904,7 +1904,7 @@ _libBornAgainDevice.vector_pvacuum_double_t_swigregister(vector_pvacuum_double_t
 
 import libBornAgainFit
 import libBornAgainBase
-class R3(arrayR3_t):
+class R3(object):
     r"""Proxy of C++ Vec3< double > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -1912,8 +1912,8 @@ class R3(arrayR3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(R3 self, double const x_, double const y_, double const z_) -> R3
         __init__(R3 self) -> R3
-        __init__(R3 self, double const x, double const y, double const z) -> R3
         """
         _libBornAgainDevice.R3_swiginit(self, _libBornAgainDevice.new_R3(*args))
 
@@ -1997,6 +1997,14 @@ class R3(arrayR3_t):
         r"""real(R3 self) -> R3"""
         return _libBornAgainDevice.R3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(R3 self, R3 other) -> bool"""
+        return _libBornAgainDevice.R3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(R3 self, R3 other) -> bool"""
+        return _libBornAgainDevice.R3___ne__(self, other)
+
     def angle(self, v):
         r"""angle(R3 self, R3 v) -> double"""
         return _libBornAgainDevice.R3_angle(self, v)
@@ -2182,7 +2190,7 @@ class vector_R3(object):
 # Register vector_R3 in _libBornAgainDevice:
 _libBornAgainDevice.vector_R3_swigregister(vector_R3)
 
-class C3(arrayC3_t):
+class C3(object):
     r"""Proxy of C++ Vec3< std::complex< double > > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -2190,8 +2198,8 @@ class C3(arrayC3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(C3 self, std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_) -> C3
         __init__(C3 self) -> C3
-        __init__(C3 self, std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3
         """
         _libBornAgainDevice.C3_swiginit(self, _libBornAgainDevice.new_C3(*args))
 
@@ -2255,6 +2263,14 @@ class C3(arrayC3_t):
         r"""real(C3 self) -> R3"""
         return _libBornAgainDevice.C3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(C3 self, C3 other) -> bool"""
+        return _libBornAgainDevice.C3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(C3 self, C3 other) -> bool"""
+        return _libBornAgainDevice.C3___ne__(self, other)
+
     def project(self, v):
         r"""project(C3 self, C3 v) -> C3"""
         return _libBornAgainDevice.C3_project(self, v)
diff --git a/auto/Wrap/libBornAgainDevice_wrap.cpp b/auto/Wrap/libBornAgainDevice_wrap.cpp
index a4349e17b27fd96a9944bb94d09bfbf448ae7055..0dfdcc980f9e765ee85477ae4091c907ade9491a 100644
--- a/auto/Wrap/libBornAgainDevice_wrap.cpp
+++ b/auto/Wrap/libBornAgainDevice_wrap.cpp
@@ -26459,20 +26459,7 @@ SWIGINTERN PyObject *vector_pvacuum_double_t_swiginit(PyObject *SWIGUNUSEDPARM(s
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< double > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< double > *)new Vec3< double >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   double arg1 ;
   double arg2 ;
@@ -26509,6 +26496,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< double > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< double > *)new Vec3< double >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -26518,7 +26518,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_R3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_R3__SWIG_0(self, argc, argv);
+    return _wrap_new_R3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -26537,7 +26537,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_R3__SWIG_1(self, argc, argv);
+          return _wrap_new_R3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -26546,8 +26546,8 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_R3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< double >::Vec3()\n"
-    "    Vec3< double >::Vec3(double const,double const,double const)\n");
+    "    Vec3< double >::Vec3(double const,double const,double const)\n"
+    "    Vec3< double >::Vec3()\n");
   return 0;
 }
 
@@ -27056,6 +27056,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_R3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___eq__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator ==((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___ne__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator !=((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_R3_angle(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
@@ -28977,20 +29047,7 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< std::complex< double > > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::complex< double > arg1 ;
   std::complex< double > arg2 ;
@@ -29027,6 +29084,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -29036,7 +29106,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_C3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_C3__SWIG_0(self, argc, argv);
+    return _wrap_new_C3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -29055,7 +29125,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_C3__SWIG_1(self, argc, argv);
+          return _wrap_new_C3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -29064,8 +29134,8 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_C3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< std::complex< double > >::Vec3()\n"
-    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n");
+    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n"
+    "    Vec3< std::complex< double > >::Vec3()\n");
   return 0;
 }
 
@@ -29459,6 +29529,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_C3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___eq__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator ==((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___ne__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator !=((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_C3_project(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
@@ -47255,8 +47395,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_pvacuum_double_t_swigregister", vector_pvacuum_double_t_swigregister, METH_O, NULL},
 	 { "vector_pvacuum_double_t_swiginit", vector_pvacuum_double_t_swiginit, METH_VARARGS, NULL},
 	 { "new_R3", _wrap_new_R3, METH_VARARGS, "\n"
-		"R3()\n"
-		"new_R3(double const x, double const y, double const z) -> R3\n"
+		"R3(double const x_, double const y_, double const z_)\n"
+		"new_R3() -> R3\n"
 		""},
 	 { "R3_x", _wrap_R3_x, METH_O, "R3_x(R3 self) -> double"},
 	 { "R3_y", _wrap_R3_y, METH_O, "R3_y(R3 self) -> double"},
@@ -47278,6 +47418,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
+	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
+	 { "R3___ne__", _wrap_R3___ne__, METH_VARARGS, "R3___ne__(R3 self, R3 other) -> bool"},
 	 { "R3_angle", _wrap_R3_angle, METH_VARARGS, "R3_angle(R3 self, R3 v) -> double"},
 	 { "R3_project", _wrap_R3_project, METH_VARARGS, "R3_project(R3 self, R3 v) -> R3"},
 	 { "delete_R3", _wrap_delete_R3, METH_O, "delete_R3(R3 self)"},
@@ -47352,8 +47494,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
 	 { "new_C3", _wrap_new_C3, METH_VARARGS, "\n"
-		"C3()\n"
-		"new_C3(std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3\n"
+		"C3(std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_)\n"
+		"new_C3() -> C3\n"
 		""},
 	 { "C3_x", _wrap_C3_x, METH_O, "C3_x(C3 self) -> std::complex< double >"},
 	 { "C3_y", _wrap_C3_y, METH_O, "C3_y(C3 self) -> std::complex< double >"},
@@ -47370,6 +47512,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
 	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
+	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
+	 { "C3___ne__", _wrap_C3___ne__, METH_VARARGS, "C3___ne__(C3 self, C3 other) -> bool"},
 	 { "C3_project", _wrap_C3_project, METH_VARARGS, "C3_project(C3 self, C3 v) -> C3"},
 	 { "delete_C3", _wrap_delete_C3, METH_O, "delete_C3(C3 self)"},
 	 { "C3_swigregister", C3_swigregister, METH_O, NULL},
@@ -49631,12 +49775,6 @@ static void *_p_FootprintGaussTo_p_INode(void *x, int *SWIGUNUSEDPARM(newmemory)
 static void *_p_IDetector2DTo_p_INode(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((INode *) (IDetector *) ((IDetector2D *) x));
 }
-static void *_p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< std::complex< double >,3 > *)  ((Vec3< std::complex< double > > *) x));
-}
-static void *_p_Vec3T_double_tTo_p_std__arrayT_double_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< double,3 > *)  ((Vec3< double > *) x));
-}
 static void *_p_ResolutionFunction2DGaussianTo_p_IResolutionFunction2D(void *x, int *SWIGUNUSEDPARM(newmemory)) {
     return (void *)((IResolutionFunction2D *)  ((ResolutionFunction2DGaussian *) x));
 }
@@ -49996,8 +50134,8 @@ static swig_cast_info _swigc__p_std__allocatorT_std__string_t[] = {  {&_swigt__p
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_unsigned_long_t[] = {  {&_swigt__p_std__allocatorT_unsigned_long_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_Vec3T_double_t, _p_Vec3T_double_tTo_p_std__arrayT_double_3_t, 0, 0},  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},  {&_swigt__p_Vec3T_std__complexT_double_t_t, _p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__complexT_double_t[] = {  {&_swigt__p_std__complexT_double_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__functionT_void_fSimulationAreaIterator_const_RF_t[] = {  {&_swigt__p_std__functionT_void_fSimulationAreaIterator_const_RF_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__invalid_argument[] = {  {&_swigt__p_std__invalid_argument, 0, 0, 0},{0, 0, 0, 0}};
diff --git a/auto/Wrap/libBornAgainResample.py b/auto/Wrap/libBornAgainResample.py
index 87fccde149a8f1adc478e1b2342efe95e11f2dd4..9c1de200acd8387fc97f8137d559fa5ae707698f 100644
--- a/auto/Wrap/libBornAgainResample.py
+++ b/auto/Wrap/libBornAgainResample.py
@@ -2060,7 +2060,7 @@ _libBornAgainResample.SimulationOptions_swigregister(SimulationOptions)
 
 import libBornAgainFit
 import libBornAgainBase
-class R3(arrayR3_t):
+class R3(object):
     r"""Proxy of C++ Vec3< double > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -2068,8 +2068,8 @@ class R3(arrayR3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(R3 self, double const x_, double const y_, double const z_) -> R3
         __init__(R3 self) -> R3
-        __init__(R3 self, double const x, double const y, double const z) -> R3
         """
         _libBornAgainResample.R3_swiginit(self, _libBornAgainResample.new_R3(*args))
 
@@ -2153,6 +2153,14 @@ class R3(arrayR3_t):
         r"""real(R3 self) -> R3"""
         return _libBornAgainResample.R3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(R3 self, R3 other) -> bool"""
+        return _libBornAgainResample.R3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(R3 self, R3 other) -> bool"""
+        return _libBornAgainResample.R3___ne__(self, other)
+
     def angle(self, v):
         r"""angle(R3 self, R3 v) -> double"""
         return _libBornAgainResample.R3_angle(self, v)
@@ -2348,7 +2356,7 @@ class vector_R3(object):
 # Register vector_R3 in _libBornAgainResample:
 _libBornAgainResample.vector_R3_swigregister(vector_R3)
 
-class C3(arrayC3_t):
+class C3(object):
     r"""Proxy of C++ Vec3< std::complex< double > > class."""
 
     thisown = property(lambda x: x.this.own(), lambda x, v: x.this.own(v), doc="The membership flag")
@@ -2356,8 +2364,8 @@ class C3(arrayC3_t):
 
     def __init__(self, *args):
         r"""
+        __init__(C3 self, std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_) -> C3
         __init__(C3 self) -> C3
-        __init__(C3 self, std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3
         """
         _libBornAgainResample.C3_swiginit(self, _libBornAgainResample.new_C3(*args))
 
@@ -2421,6 +2429,14 @@ class C3(arrayC3_t):
         r"""real(C3 self) -> R3"""
         return _libBornAgainResample.C3_real(self)
 
+    def __eq__(self, other):
+        r"""__eq__(C3 self, C3 other) -> bool"""
+        return _libBornAgainResample.C3___eq__(self, other)
+
+    def __ne__(self, other):
+        r"""__ne__(C3 self, C3 other) -> bool"""
+        return _libBornAgainResample.C3___ne__(self, other)
+
     def project(self, v):
         r"""project(C3 self, C3 v) -> C3"""
         return _libBornAgainResample.C3_project(self, v)
diff --git a/auto/Wrap/libBornAgainResample_wrap.cpp b/auto/Wrap/libBornAgainResample_wrap.cpp
index dc0b8e481e5df76610943e68530d6aa582299f56..1dd3d81da1bf343b4db7dea141a40e5731cfabaf 100644
--- a/auto/Wrap/libBornAgainResample_wrap.cpp
+++ b/auto/Wrap/libBornAgainResample_wrap.cpp
@@ -26902,20 +26902,7 @@ SWIGINTERN PyObject *SimulationOptions_swiginit(PyObject *SWIGUNUSEDPARM(self),
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< double > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< double > *)new Vec3< double >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   double arg1 ;
   double arg2 ;
@@ -26952,6 +26939,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_R3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< double > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< double > *)new Vec3< double >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_double_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -26961,7 +26961,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_R3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_R3__SWIG_0(self, argc, argv);
+    return _wrap_new_R3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -26980,7 +26980,7 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_R3__SWIG_1(self, argc, argv);
+          return _wrap_new_R3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -26989,8 +26989,8 @@ SWIGINTERN PyObject *_wrap_new_R3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_R3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< double >::Vec3()\n"
-    "    Vec3< double >::Vec3(double const,double const,double const)\n");
+    "    Vec3< double >::Vec3(double const,double const,double const)\n"
+    "    Vec3< double >::Vec3()\n");
   return 0;
 }
 
@@ -27499,6 +27499,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_R3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___eq__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___eq__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator ==((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_R3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< double > *arg1 = (Vec3< double > *) 0 ;
+  Vec3< double > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "R3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_double_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "R3___ne__" "', argument " "1"" of type '" "Vec3< double > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< double > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_double_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "R3___ne__" "', argument " "2"" of type '" "Vec3< double > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< double > * >(argp2);
+  result = (bool)((Vec3< double > const *)arg1)->operator !=((Vec3< double > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_R3_angle(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< double > *arg1 = (Vec3< double > *) 0 ;
@@ -29544,20 +29614,7 @@ SWIGINTERN PyObject *vector_R3_swiginit(PyObject *SWIGUNUSEDPARM(self), PyObject
   return SWIG_Python_InitShadowInstance(args);
 }
 
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
-  PyObject *resultobj = 0;
-  Vec3< std::complex< double > > *result = 0 ;
-  
-  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
-  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
-  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
-SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_0(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **swig_obj) {
   PyObject *resultobj = 0;
   std::complex< double > arg1 ;
   std::complex< double > arg2 ;
@@ -29594,6 +29651,19 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_new_C3__SWIG_1(PyObject *SWIGUNUSEDPARM(self), Py_ssize_t nobjs, PyObject **SWIGUNUSEDPARM(swig_obj)) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *result = 0 ;
+  
+  if ((nobjs < 0) || (nobjs > 0)) SWIG_fail;
+  result = (Vec3< std::complex< double > > *)new Vec3< std::complex< double > >();
+  resultobj = SWIG_NewPointerObj(SWIG_as_voidptr(result), SWIGTYPE_p_Vec3T_std__complexT_double_t_t, SWIG_POINTER_NEW |  0 );
+  return resultobj;
+fail:
+  return NULL;
+}
+
+
 SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   Py_ssize_t argc;
   PyObject *argv[4] = {
@@ -29603,7 +29673,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
   if (!(argc = SWIG_Python_UnpackTuple(args, "new_C3", 0, 3, argv))) SWIG_fail;
   --argc;
   if (argc == 0) {
-    return _wrap_new_C3__SWIG_0(self, argc, argv);
+    return _wrap_new_C3__SWIG_1(self, argc, argv);
   }
   if (argc == 3) {
     int _v;
@@ -29622,7 +29692,7 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
           _v = SWIG_CheckState(res);
         }
         if (_v) {
-          return _wrap_new_C3__SWIG_1(self, argc, argv);
+          return _wrap_new_C3__SWIG_0(self, argc, argv);
         }
       }
     }
@@ -29631,8 +29701,8 @@ SWIGINTERN PyObject *_wrap_new_C3(PyObject *self, PyObject *args) {
 fail:
   SWIG_Python_RaiseOrModifyTypeError("Wrong number or type of arguments for overloaded function 'new_C3'.\n"
     "  Possible C/C++ prototypes are:\n"
-    "    Vec3< std::complex< double > >::Vec3()\n"
-    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n");
+    "    Vec3< std::complex< double > >::Vec3(std::complex< double > const,std::complex< double > const,std::complex< double > const)\n"
+    "    Vec3< std::complex< double > >::Vec3()\n");
   return 0;
 }
 
@@ -30026,6 +30096,76 @@ fail:
 }
 
 
+SWIGINTERN PyObject *_wrap_C3___eq__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___eq__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___eq__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___eq__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator ==((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
+SWIGINTERN PyObject *_wrap_C3___ne__(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
+  PyObject *resultobj = 0;
+  Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
+  Vec3< std::complex< double > > *arg2 = 0 ;
+  void *argp1 = 0 ;
+  int res1 = 0 ;
+  void *argp2 = 0 ;
+  int res2 = 0 ;
+  PyObject *swig_obj[2] ;
+  bool result;
+  
+  if (!SWIG_Python_UnpackTuple(args, "C3___ne__", 2, 2, swig_obj)) SWIG_fail;
+  res1 = SWIG_ConvertPtr(swig_obj[0], &argp1,SWIGTYPE_p_Vec3T_std__complexT_double_t_t, 0 |  0 );
+  if (!SWIG_IsOK(res1)) {
+    SWIG_exception_fail(SWIG_ArgError(res1), "in method '" "C3___ne__" "', argument " "1"" of type '" "Vec3< std::complex< double > > const *""'"); 
+  }
+  arg1 = reinterpret_cast< Vec3< std::complex< double > > * >(argp1);
+  res2 = SWIG_ConvertPtr(swig_obj[1], &argp2, SWIGTYPE_p_Vec3T_std__complexT_double_t_t,  0  | 0);
+  if (!SWIG_IsOK(res2)) {
+    SWIG_exception_fail(SWIG_ArgError(res2), "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  if (!argp2) {
+    SWIG_exception_fail(SWIG_ValueError, "invalid null reference " "in method '" "C3___ne__" "', argument " "2"" of type '" "Vec3< std::complex< double > > const &""'"); 
+  }
+  arg2 = reinterpret_cast< Vec3< std::complex< double > > * >(argp2);
+  result = (bool)((Vec3< std::complex< double > > const *)arg1)->operator !=((Vec3< std::complex< double > > const &)*arg2);
+  resultobj = SWIG_From_bool(static_cast< bool >(result));
+  return resultobj;
+fail:
+  PyErr_Clear();
+  Py_INCREF(Py_NotImplemented);
+  return Py_NotImplemented;
+}
+
+
 SWIGINTERN PyObject *_wrap_C3_project(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   Vec3< std::complex< double > > *arg1 = (Vec3< std::complex< double > > *) 0 ;
@@ -32885,8 +33025,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "SimulationOptions_swigregister", SimulationOptions_swigregister, METH_O, NULL},
 	 { "SimulationOptions_swiginit", SimulationOptions_swiginit, METH_VARARGS, NULL},
 	 { "new_R3", _wrap_new_R3, METH_VARARGS, "\n"
-		"R3()\n"
-		"new_R3(double const x, double const y, double const z) -> R3\n"
+		"R3(double const x_, double const y_, double const z_)\n"
+		"new_R3() -> R3\n"
 		""},
 	 { "R3_x", _wrap_R3_x, METH_O, "R3_x(R3 self) -> double"},
 	 { "R3_y", _wrap_R3_y, METH_O, "R3_y(R3 self) -> double"},
@@ -32908,6 +33048,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "R3_unit", _wrap_R3_unit, METH_O, "R3_unit(R3 self) -> R3"},
 	 { "R3_complex", _wrap_R3_complex, METH_O, "R3_complex(R3 self) -> C3"},
 	 { "R3_real", _wrap_R3_real, METH_O, "R3_real(R3 self) -> R3"},
+	 { "R3___eq__", _wrap_R3___eq__, METH_VARARGS, "R3___eq__(R3 self, R3 other) -> bool"},
+	 { "R3___ne__", _wrap_R3___ne__, METH_VARARGS, "R3___ne__(R3 self, R3 other) -> bool"},
 	 { "R3_angle", _wrap_R3_angle, METH_VARARGS, "R3_angle(R3 self, R3 v) -> double"},
 	 { "R3_project", _wrap_R3_project, METH_VARARGS, "R3_project(R3 self, R3 v) -> R3"},
 	 { "R3___add__", _wrap_R3___add__, METH_VARARGS, "R3___add__(R3 self, R3 rhs) -> R3"},
@@ -32980,8 +33122,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "vector_R3_swigregister", vector_R3_swigregister, METH_O, NULL},
 	 { "vector_R3_swiginit", vector_R3_swiginit, METH_VARARGS, NULL},
 	 { "new_C3", _wrap_new_C3, METH_VARARGS, "\n"
-		"C3()\n"
-		"new_C3(std::complex< double > const x, std::complex< double > const y, std::complex< double > const z) -> C3\n"
+		"C3(std::complex< double > const x_, std::complex< double > const y_, std::complex< double > const z_)\n"
+		"new_C3() -> C3\n"
 		""},
 	 { "C3_x", _wrap_C3_x, METH_O, "C3_x(C3 self) -> std::complex< double >"},
 	 { "C3_y", _wrap_C3_y, METH_O, "C3_y(C3 self) -> std::complex< double >"},
@@ -32998,6 +33140,8 @@ static PyMethodDef SwigMethods[] = {
 	 { "C3_magxy", _wrap_C3_magxy, METH_O, "C3_magxy(C3 self) -> double"},
 	 { "C3_unit", _wrap_C3_unit, METH_O, "C3_unit(C3 self) -> C3"},
 	 { "C3_real", _wrap_C3_real, METH_O, "C3_real(C3 self) -> R3"},
+	 { "C3___eq__", _wrap_C3___eq__, METH_VARARGS, "C3___eq__(C3 self, C3 other) -> bool"},
+	 { "C3___ne__", _wrap_C3___ne__, METH_VARARGS, "C3___ne__(C3 self, C3 other) -> bool"},
 	 { "C3_project", _wrap_C3_project, METH_VARARGS, "C3_project(C3 self, C3 v) -> C3"},
 	 { "delete_C3", _wrap_delete_C3, METH_O, "delete_C3(C3 self)"},
 	 { "C3_swigregister", C3_swigregister, METH_O, NULL},
@@ -33098,12 +33242,6 @@ static PyMethodDef SwigMethods_proxydocs[] = {
 
 /* -------- TYPE CONVERSION AND EQUIVALENCE RULES (BEGIN) -------- */
 
-static void *_p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< std::complex< double >,3 > *)  ((Vec3< std::complex< double > > *) x));
-}
-static void *_p_Vec3T_double_tTo_p_std__arrayT_double_3_t(void *x, int *SWIGUNUSEDPARM(newmemory)) {
-    return (void *)((std::array< double,3 > *)  ((Vec3< double > *) x));
-}
 static swig_type_info _swigt__p_MultiLayer = {"_p_MultiLayer", "MultiLayer *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_SimulationOptions = {"_p_SimulationOptions", "SimulationOptions *", 0, 0, (void*)0, 0};
 static swig_type_info _swigt__p_ThreadInfo = {"_p_ThreadInfo", "ThreadInfo *", 0, 0, (void*)0, 0};
@@ -33244,8 +33382,8 @@ static swig_cast_info _swigc__p_std__allocatorT_std__string_t[] = {  {&_swigt__p
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_double_std__allocatorT_double_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t[] = {  {&_swigt__p_std__allocatorT_std__vectorT_int_std__allocatorT_int_t_t_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__allocatorT_unsigned_long_t[] = {  {&_swigt__p_std__allocatorT_unsigned_long_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_Vec3T_double_t, _p_Vec3T_double_tTo_p_std__arrayT_double_3_t, 0, 0},  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
-static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},  {&_swigt__p_Vec3T_std__complexT_double_t_t, _p_Vec3T_std__complexT_double_t_tTo_p_std__arrayT_std__complexT_double_t_3_t, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_double_3_t[] = {  {&_swigt__p_std__arrayT_double_3_t, 0, 0, 0},{0, 0, 0, 0}};
+static swig_cast_info _swigc__p_std__arrayT_std__complexT_double_t_3_t[] = {  {&_swigt__p_std__arrayT_std__complexT_double_t_3_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__complexT_double_t[] = {  {&_swigt__p_std__complexT_double_t, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__invalid_argument[] = {  {&_swigt__p_std__invalid_argument, 0, 0, 0},{0, 0, 0, 0}};
 static swig_cast_info _swigc__p_std__lessT_std__string_t[] = {  {&_swigt__p_std__lessT_std__string_t, 0, 0, 0},{0, 0, 0, 0}};
diff --git a/cmake/BornAgain/Dependences.cmake b/cmake/BornAgain/Dependences.cmake
index 875a1e2e03cc59865e7418cedd58eec881fb4b92..583d35375c7ae19ac83ebd9774723a5c777b0d09 100644
--- a/cmake/BornAgain/Dependences.cmake
+++ b/cmake/BornAgain/Dependences.cmake
@@ -8,6 +8,9 @@ find_package(LibHeinz REQUIRED)
 message(STATUS "LibHeinz: found=${LibHeinz_FOUND}, include_dirs=${LibHeinz_INCLUDE_DIR}, "
     "version=${LibHeinz_VERSION}")
 
+find_package(formfactor REQUIRED)
+message(STATUS "formfactor: found=${formfactor_FOUND}, version=${formfactor_VERSION}")
+
 find_package(Threads REQUIRED)
 find_package(FFTW3 REQUIRED)