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..10dd1d93ee92d5a48cbcfbc5ef73fbad13182766 100644
--- a/Sample/HardParticle/IFormFactorPolyhedron.cpp
+++ b/Sample/HardParticle/IFormFactorPolyhedron.cpp
@@ -17,7 +17,7 @@
 //! "Form factor (Fourier shape transform) of polygon and polyhedron."
 
 #include "Sample/HardParticle/IFormFactorPolyhedron.h"
-#include "Sample/LibFF/Polyhedron.h"
+#include <ff/Polyhedron.h>
 
 // #ifdef ALGORITHM_DIAGNOSTIC // TODO restore
 // void IFormFactorPolyhedron::setLimits(double _q, int _n)
@@ -37,10 +37,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.h b/Sample/HardParticle/Prism.h
index 4d1ac331dd281f3a56f7b58d576f9232e74850bb..e464053386c56aa9d458a51c815aeabf2500193c 100644
--- a/Sample/HardParticle/Prism.h
+++ b/Sample/HardParticle/Prism.h
@@ -21,7 +21,7 @@
 #define BORNAGAIN_SAMPLE_HARDPARTICLE_PRISM_H
 
 #include "Sample/LibFF/PolyhedralComponents.h"
-#include "Sample/LibFF/PolyhedralTopology.h"
+#include <ff/PolyhedralTopology.h>
 #include <memory>
 
 class Prism {
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/cmake/BornAgain/Dependences.cmake b/cmake/BornAgain/Dependences.cmake
index 875a1e2e03cc59865e7418cedd58eec881fb4b92..f48b75045cc0c094e31f26895299fff098b62d62 100644
--- a/cmake/BornAgain/Dependences.cmake
+++ b/cmake/BornAgain/Dependences.cmake
@@ -8,6 +8,10 @@ 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}, include_dirs=${formfactor_INCLUDE_DIR}, "
+    "version=${formfactor_VERSION}")
+
 find_package(Threads REQUIRED)
 find_package(FFTW3 REQUIRED)