diff --git a/Resample/Element/DiffuseElement.cpp b/Resample/Element/DiffuseElement.cpp
index 2f8b3b4acebcff2986814b42cde79cac9da86253..bc35ca3f69f58b1a74b37832836014f6d4b7bfd6 100644
--- a/Resample/Element/DiffuseElement.cpp
+++ b/Resample/Element/DiffuseElement.cpp
@@ -61,12 +61,12 @@ void DiffuseElement::setFluxes(const Fluxes* fluxes_in, const Fluxes* fluxes_out
 
 const IFlux* DiffuseElement::fluxIn(size_t i_layer) const
 {
-    return (*m_fluxes_in)[i_layer].get();
+    return (*m_fluxes_in)[i_layer];
 }
 
 const IFlux* DiffuseElement::fluxOut(size_t i_layer) const
 {
-    return (*m_fluxes_out)[i_layer].get();
+    return (*m_fluxes_out)[i_layer];
 }
 
 DiffuseElement DiffuseElement::pointElement(double x, double y) const
diff --git a/Resample/Element/DiffuseElement.h b/Resample/Element/DiffuseElement.h
index d110d3122d30268280bde8f4168c6d392a99e5c7..6b088549e6bbf6c8025f5cbe9ad57d51c7efa7c9 100644
--- a/Resample/Element/DiffuseElement.h
+++ b/Resample/Element/DiffuseElement.h
@@ -25,12 +25,11 @@
 #include <memory>
 #include <vector>
 
+class Fluxes;
 class IFlux;
 class IPixel;
 class WavevectorInfo;
 
-using Fluxes = std::vector<std::unique_ptr<const IFlux>>;
-
 //! Data stucture containing both input and output of a single detector cell.
 //! @ingroup simulation
 
diff --git a/Resample/Flux/IFlux.h b/Resample/Flux/IFlux.h
index 27c0dfc399b51e94320bfcea459653456ed72a4a..dd8e96a71c2655222ed370ada1f5762fc6b82ddf 100644
--- a/Resample/Flux/IFlux.h
+++ b/Resample/Flux/IFlux.h
@@ -21,8 +21,7 @@
 #define BORNAGAIN_RESAMPLE_FLUX_IFLUX_H
 
 #include "Base/Spin/Spinor.h"
-#include <memory>
-#include <vector>
+#include "Base/Types/OwningVector.h"
 
 //! Interface to access reflection/transmission coefficients.
 //! Realized by ScalarFlux and MatrixFlux.
@@ -45,5 +44,7 @@ public:
     virtual Spinor getKz() const = 0;
 };
 
+class Fluxes : public OwningVector<IFlux> {};
+
 #endif // BORNAGAIN_RESAMPLE_FLUX_IFLUX_H
 #endif // USER_API
diff --git a/Resample/Processed/ReSample.h b/Resample/Processed/ReSample.h
index a9596c5ede56aec42568e00f132dfb1e4f613b0f..2b71f47c6c1f8d321c5e101516488909ff63b765 100644
--- a/Resample/Processed/ReSample.h
+++ b/Resample/Processed/ReSample.h
@@ -25,14 +25,12 @@
 #include <memory>
 #include <vector>
 
-class IFlux;
+class Fluxes;
 class Material;
 class MultiLayer;
 class reLayout;
 class SimulationOptions;
 
-using Fluxes = std::vector<std::unique_ptr<const IFlux>>;
-
 //! Data structure that contains all the necessary data for scattering calculations.
 //!
 //! If the usage of average materials is requested, layers and particles are sliced into multiple
diff --git a/Resample/Specular/ComputeFluxMagnetic.cpp b/Resample/Specular/ComputeFluxMagnetic.cpp
index 470c9c2e1421f4bd2e675ca54d0136fed79c6482..e15980dcd9d95dc0b375e7d5d3d10d1332547026 100644
--- a/Resample/Specular/ComputeFluxMagnetic.cpp
+++ b/Resample/Specular/ComputeFluxMagnetic.cpp
@@ -183,7 +183,7 @@ Fluxes Compute::SpecularMagnetic::fluxes(const SliceStack& slices, const R3& k,
 
     Fluxes result;
     for (auto& coeff : computeFlux(slices, kz, slices.roughnessModel(), forward))
-        result.emplace_back(std::make_unique<const MatrixFlux>(coeff));
+        result.emplace_back(new MatrixFlux(coeff));
 
     return result;
 }
diff --git a/Resample/Specular/ComputeFluxMagnetic.h b/Resample/Specular/ComputeFluxMagnetic.h
index da3201d69ba39f16f9e04563df882c424cebf44f..256cac3e3b7930e9916599c08e0acfd1f64529fb 100644
--- a/Resample/Specular/ComputeFluxMagnetic.h
+++ b/Resample/Specular/ComputeFluxMagnetic.h
@@ -21,15 +21,12 @@
 #define BORNAGAIN_RESAMPLE_SPECULAR_COMPUTEFLUXMAGNETIC_H
 
 #include "Base/Spin/SpinMatrix.h"
-#include "Sample/Multilayer/RoughnessModels.h"
 #include <heinz/Complex.h>
 #include <heinz/Vectors3D.h>
-#include <memory>
 #include <vector>
 
-class IFlux;
+class Fluxes;
 class SliceStack;
-using Fluxes = std::vector<std::unique_ptr<const IFlux>>;
 
 //! Methods to compute polarized propagation directions and fluxes as function of slice.
 //!
diff --git a/Resample/Specular/ComputeFluxScalar.cpp b/Resample/Specular/ComputeFluxScalar.cpp
index 450fd80e76a17c2dea73e0e23aa2c3a192bbc6b9..8ee5f1d983361489009a6c35e60110b49b733083 100644
--- a/Resample/Specular/ComputeFluxScalar.cpp
+++ b/Resample/Specular/ComputeFluxScalar.cpp
@@ -129,7 +129,7 @@ Fluxes Compute::SpecularScalar::fluxes(const SliceStack& slices, const R3& k)
 
     Fluxes result;
     for (size_t i = 0; i < kz.size(); ++i)
-        result.emplace_back(std::make_unique<const ScalarFlux>(kz[i], TR[i]));
+        result.emplace_back(new ScalarFlux(kz[i], TR[i]));
     return result;
 }
 
diff --git a/Resample/Specular/ComputeFluxScalar.h b/Resample/Specular/ComputeFluxScalar.h
index f434dfd8a5cc2aec2e9829dff854eb8e7dc32b15..49683c1821669fe0e09722e35c320606908b393b 100644
--- a/Resample/Specular/ComputeFluxScalar.h
+++ b/Resample/Specular/ComputeFluxScalar.h
@@ -23,12 +23,10 @@
 #include "Sample/Multilayer/RoughnessModels.h"
 #include <heinz/Complex.h>
 #include <heinz/Vectors3D.h>
-#include <memory>
 #include <vector>
 
-class IFlux;
+class Fluxes;
 class SliceStack;
-using Fluxes = std::vector<std::unique_ptr<const IFlux>>;
 
 //! Methods to compute scalar propagation directions and fluxes as function of slice.
 
diff --git a/Sim/Computation/DepthProbeComputation.cpp b/Sim/Computation/DepthProbeComputation.cpp
index 75ae8e4b58e8e8618403dd552bf250ec80b93566..0b37ab54655bdb45356940627d6302f201f3d06a 100644
--- a/Sim/Computation/DepthProbeComputation.cpp
+++ b/Sim/Computation/DepthProbeComputation.cpp
@@ -55,7 +55,7 @@ void DepthProbeComputation::runProtected()
             z_layer_top = i_layer ? m_re_sample.sliceTopZ(i_layer) : 0;
 
             // get R & T coefficients for current layer
-            const auto* flux = dynamic_cast<const ScalarFlux*>(fluxes[i_layer].get());
+            const auto* flux = dynamic_cast<const ScalarFlux*>(fluxes[i_layer]);
             ASSERT(flux);
             const complex_t R = flux->getScalarR();
             const complex_t T = flux->getScalarT();
diff --git a/Tests/Unit/Resample/RTTest.cpp b/Tests/Unit/Resample/RTTest.cpp
index 12a3d9a1a33aface6a917b26b30bca5cd53ac54b..1f8cb77fa4f8d735e8d4c6626d26b782a0583310 100644
--- a/Tests/Unit/Resample/RTTest.cpp
+++ b/Tests/Unit/Resample/RTTest.cpp
@@ -33,7 +33,7 @@ protected:
     {
         std::vector<ScalarFlux> result;
         for (auto& flux : fluxes)
-            result.push_back(*dynamic_cast<const ScalarFlux*>(flux.get()));
+            result.push_back(*dynamic_cast<const ScalarFlux*>(flux));
         return result;
     }
     const Material air = RefractiveMaterial("Air", 1e-8, 1e-8);
diff --git a/auto/Wrap/doxygenResample.i b/auto/Wrap/doxygenResample.i
index 91bde99fccab2c50e69bdb1777165d30e05306b7..34aea93e9609666f76abea53af011eaedaa01c2c 100644
--- a/auto/Wrap/doxygenResample.i
+++ b/auto/Wrap/doxygenResample.i
@@ -130,6 +130,10 @@ Tells if simulation element corresponds to a specular peak.
 ";
 
 
+// File: classFluxes.xml
+%feature("docstring") Fluxes "";
+
+
 // File: classIFlux.xml
 %feature("docstring") IFlux "