diff --git a/Resample/FFCompute/ComputeBA.cpp b/Resample/FFCompute/ComputeBA.cpp
index 182396945fdec94fd6a2bb6855c21b68d87a70c4..d52c177c3c714bc27ca2bcf5591ba3cec37d7cd3 100644
--- a/Resample/FFCompute/ComputeBA.cpp
+++ b/Resample/FFCompute/ComputeBA.cpp
@@ -25,7 +25,9 @@ ComputeBA* ComputeBA::clone() const
     return new ComputeBA(*m_ff, m_iLayer);
 }
 
-complex_t ComputeBA::computeFF(const WavevectorInfo& wavevectors) const
+complex_t ComputeBA::computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const
 {
     return m_ff->theFF(wavevectors);
 }
diff --git a/Resample/FFCompute/ComputeBA.h b/Resample/FFCompute/ComputeBA.h
index 27ffcbbe2a2bc8730c87942a301e594702ff4cfd..a74cc5ffc10643e8c106cb04439cf6123f64beac 100644
--- a/Resample/FFCompute/ComputeBA.h
+++ b/Resample/FFCompute/ComputeBA.h
@@ -35,7 +35,9 @@ public:
     ComputeBA* clone() const override;
 
     //! Calculates and returns a form factor calculation in BA
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 };
 
 #endif // BORNAGAIN_RESAMPLE_FFCOMPUTE_COMPUTEBA_H
diff --git a/Resample/FFCompute/ComputeBAPol.cpp b/Resample/FFCompute/ComputeBAPol.cpp
index 44c6fe2f054d879bbc2a001fa3edf761431a74f6..1e3184e9df1b9b6ce941cbc26c112022d4924af6 100644
--- a/Resample/FFCompute/ComputeBAPol.cpp
+++ b/Resample/FFCompute/ComputeBAPol.cpp
@@ -25,7 +25,9 @@ ComputeBAPol* ComputeBAPol::clone() const
     return new ComputeBAPol(*m_ff, m_iLayer);
 }
 
-complex_t ComputeBAPol::computeFF(const WavevectorInfo&) const
+complex_t ComputeBAPol::computeFF(const WavevectorInfo&,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const
 {
     throw std::runtime_error("ComputeBAPol::evaluate: "
                              "should never be called for matrix interactions");
diff --git a/Resample/FFCompute/ComputeBAPol.h b/Resample/FFCompute/ComputeBAPol.h
index 408e3009c97032dce976cf95a4392dc22cf99b03..d77f2dfd917527b52436301451bc651f24a6d45a 100644
--- a/Resample/FFCompute/ComputeBAPol.h
+++ b/Resample/FFCompute/ComputeBAPol.h
@@ -36,7 +36,9 @@ public:
     ComputeBAPol* clone() const override;
 
     //! Throws not-implemented exception
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 
     //! Calculates and returns a polarized form factor calculation in BA
     Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const override;
diff --git a/Resample/FFCompute/ComputeDWBA.cpp b/Resample/FFCompute/ComputeDWBA.cpp
index f62687126cd4a01ad2d491a11599f897ca1e8d19..4567a0b2b2d74bb07365d3268b333d9c7c0b85bf 100644
--- a/Resample/FFCompute/ComputeDWBA.cpp
+++ b/Resample/FFCompute/ComputeDWBA.cpp
@@ -23,24 +23,22 @@ ComputeDWBA::~ComputeDWBA() = default;
 
 ComputeDWBA* ComputeDWBA::clone() const
 {
-    ComputeDWBA* result = new ComputeDWBA(*m_ff, m_iLayer);
-    auto in_coefs = std::unique_ptr<const IFlux>(m_inFlux ? m_inFlux->clone() : nullptr);
-    auto out_coefs = std::unique_ptr<const IFlux>(m_outFlux ? m_outFlux->clone() : nullptr);
-    result->setFlux(std::move(in_coefs), std::move(out_coefs));
-    return result;
+    return new ComputeDWBA(*m_ff, m_iLayer);
 }
 
-complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
+complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const
 {
     // Retrieve the two different incoming wavevectors in the layer
     cvector_t k_i_T = wavevectors.getKi();
-    k_i_T.setZ(-m_inFlux->getScalarKz());
+    k_i_T.setZ(-inFlux->getScalarKz());
     cvector_t k_i_R = k_i_T;
     k_i_R.setZ(-k_i_T.z());
 
     // Retrieve the two different outgoing wavevector bins in the layer
     cvector_t k_f_T = wavevectors.getKf();
-    k_f_T.setZ(m_outFlux->getScalarKz());
+    k_f_T.setZ(outFlux->getScalarKz());
     cvector_t k_f_R = k_f_T;
     k_f_R.setZ(-k_f_T.z());
 
@@ -52,10 +50,10 @@ complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
     WavevectorInfo k_RR(k_i_R, k_f_R, wavelength);
 
     // Get the four R,T coefficients
-    complex_t T_in = m_inFlux->getScalarT();
-    complex_t R_in = m_inFlux->getScalarR();
-    complex_t T_out = m_outFlux->getScalarT();
-    complex_t R_out = m_outFlux->getScalarR();
+    complex_t T_in = inFlux->getScalarT();
+    complex_t R_in = inFlux->getScalarR();
+    complex_t T_out = outFlux->getScalarT();
+    complex_t R_out = outFlux->getScalarR();
 
     // The four different scattering contributions; S stands for scattering
     // off the particle, R for reflection off the layer interface
@@ -66,9 +64,3 @@ complex_t ComputeDWBA::computeFF(const WavevectorInfo& wavevectors) const
 
     return term_S + term_RS + term_SR + term_RSR;
 }
-
-void ComputeDWBA::setFlux(std::unique_ptr<const IFlux> inFlux, std::unique_ptr<const IFlux> outFlux)
-{
-    m_inFlux = std::move(inFlux);
-    m_outFlux = std::move(outFlux);
-}
diff --git a/Resample/FFCompute/ComputeDWBA.h b/Resample/FFCompute/ComputeDWBA.h
index 108ae3b0c26bb7adbd8f7be92a9dc39aaa137ade..73936d9ad2dc8a789deeca343d317746f7116f5c 100644
--- a/Resample/FFCompute/ComputeDWBA.h
+++ b/Resample/FFCompute/ComputeDWBA.h
@@ -37,10 +37,9 @@ public:
     ComputeDWBA* clone() const override;
 
     //! Returns the coherent sum of the four DWBA terms for scalar scattering.
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
-
-    void setFlux(std::unique_ptr<const IFlux> inFlux,
-                 std::unique_ptr<const IFlux> outFlux) override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 
     friend class TestPolarizedDWBATerms;
 
diff --git a/Resample/FFCompute/ComputeDWBAPol.cpp b/Resample/FFCompute/ComputeDWBAPol.cpp
index d0be17ca072b5eb300aadc2ad59ec3dbe7cd1d2a..bacc851fc1f686f17ff767d532ab66891b118baa 100644
--- a/Resample/FFCompute/ComputeDWBAPol.cpp
+++ b/Resample/FFCompute/ComputeDWBAPol.cpp
@@ -43,7 +43,9 @@ ComputeDWBAPol* ComputeDWBAPol::clone() const
     return result;
 }
 
-complex_t ComputeDWBAPol::computeFF(const WavevectorInfo&) const
+complex_t ComputeDWBAPol::computeFF(const WavevectorInfo&,
+                        const std::unique_ptr<const IFlux>&,
+                        const std::unique_ptr<const IFlux>&) const
 {
     throw std::runtime_error("Bug: forbidden call of ComputeDWBAPol::evaluate");
 }
diff --git a/Resample/FFCompute/ComputeDWBAPol.h b/Resample/FFCompute/ComputeDWBAPol.h
index eeb05f46aa44e4895b1dbcc114c5222657a0c42c..20144d2a0b754d5d499d2c96b65d267745793ff7 100644
--- a/Resample/FFCompute/ComputeDWBAPol.h
+++ b/Resample/FFCompute/ComputeDWBAPol.h
@@ -37,7 +37,9 @@ public:
     ComputeDWBAPol* clone() const override;
 
     //! Throws not-implemented exception.
-    complex_t computeFF(const WavevectorInfo& wavevectors) const override;
+    complex_t computeFF(const WavevectorInfo& wavevectors,
+                        const std::unique_ptr<const IFlux>& inFlux,
+                        const std::unique_ptr<const IFlux>& outFlux) const override;
 
     //! Returns the coherent sum of the four DWBA terms for polarized scattering.
     Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const override;
diff --git a/Resample/FFCompute/IComputeFF.cpp b/Resample/FFCompute/IComputeFF.cpp
index e7301f36ddcdcad46aa4b264b5cce573aeae344a..818317ab77ed99e0b9a84eea304707ca4e207580 100644
--- a/Resample/FFCompute/IComputeFF.cpp
+++ b/Resample/FFCompute/IComputeFF.cpp
@@ -60,8 +60,7 @@ complex_t IComputeFF::coherentFF(const SimulationElement& ele)
     ASSERT(fresnel_map);
     auto inFlux = fresnel_map->getInFlux(ele.getKi(), m_iLayer);
     auto outFlux = fresnel_map->getOutFlux(ele.getMeanKf(), m_iLayer);
-    setFlux(std::move(inFlux), std::move(outFlux));
-    return computeFF({ele});
+    return computeFF({ele}, inFlux, outFlux);
 }
 
 Eigen::Matrix2cd IComputeFF::coherentPolFF(const SimulationElement& ele)
diff --git a/Resample/FFCompute/IComputeFF.h b/Resample/FFCompute/IComputeFF.h
index 6d0cfffd2ca3c2fcfe42b931c879116806d0161f..c8e7e4e963cb636f787cb7078d664d59c95745d1 100644
--- a/Resample/FFCompute/IComputeFF.h
+++ b/Resample/FFCompute/IComputeFF.h
@@ -48,7 +48,9 @@ public:
     virtual double radialExtension() const;
     virtual double bottomZ(const IRotation& rotation) const;
     virtual double topZ(const IRotation& rotation) const;
-    virtual complex_t computeFF(const WavevectorInfo& wavevectors) const = 0;
+    virtual complex_t computeFF(const WavevectorInfo& wavevectors,
+                                const std::unique_ptr<const IFlux>& inFlux,
+                                const std::unique_ptr<const IFlux>& outFlux) const = 0;
     //! Returns scattering amplitude for matrix interactions
     virtual Eigen::Matrix2cd computePolFF(const WavevectorInfo& wavevectors) const;
     //! Sets reflection/transmission info