diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index e2c42bf9e8f3925735b15edd27b9e339ef6aaa67..110e51e90cbadde1bd6ac7eca0a4564b95c8a310 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -410,6 +410,12 @@ class TEST(_code):
         if not os.path.exists(os.path.join(self.work_dir, self.simname + '.h5')):
             self.write_par(
                     particle_ic = particle_initial_condition)
+            if self.dns_type == 'test_tracer_set':
+                ofile = h5py.File(os.path.join(self.work_dir, 'test_particle_sample.h5'), 'w')
+                ofile.create_group('tracers0')
+                ofile.create_group('tracers0/position')
+                ofile.create_group('tracers0/velocity')
+                ofile.close()
             if self.dns_type == 'test_interpolation':
                 if type(particle_initial_condition) == type(None):
                     pbase_shape = (self.parameters['nparticles'],)
diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index a9cefdc1e2d3601f5679787f1dfd4db8bf5fd71c..094582bc9bc747611c8c30b3341ffca7339fc50e 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -102,26 +102,59 @@ int test_tracer_set<rnumber>::do_work(void)
     pset.init(bla);
 
     // initialize particle output object
-    particles_output_hdf5<long long int, double, 3> *particles_output_writer_mpi = new particles_output_hdf5<
+    particles_output_hdf5<long long int, double, 3> *particle_output_writer = new particles_output_hdf5<
         long long int, double, 3>(
                 MPI_COMM_WORLD,
                 "tracers0",
                 pset.getTotalNumberOfParticles(),
                 0);
+    particles_output_sampling_hdf5<long long int, double, 3> *particle_sample_writer = new particles_output_sampling_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                pset.getTotalNumberOfParticles(),
+                "test_particle_sample.h5",
+                "tracers0",
+                "position/0");
 
     fti->set_field(vec_field);
     trhs->setVelocity(fti);
 
     particle_solver psolver(pset, 0);
 
+    particle_output_writer->open_file("test_particle_checkpoint.h5");
+    psolver.template writeCheckpoint<3>(particle_output_writer);
+    pset.writeSample(
+            vec_field,
+            particle_sample_writer,
+            "tracers0",
+            "velocity",
+            psolver.getIteration());
+    pset.writeStateTriplet(
+            0,
+            particle_sample_writer,
+            "tracers0",
+            "position",
+            psolver.getIteration());
     psolver.Euler(0.001, *trhs);
     psolver.setIteration(1);
-    particles_output_writer_mpi->open_file("test_particle_checkpoint.h5");
-    psolver.template writeCheckpoint<3>(particles_output_writer_mpi);
-    particles_output_writer_mpi->close_file();
+    psolver.template writeCheckpoint<3>(particle_output_writer);
+    pset.writeSample(
+            vec_field,
+            particle_sample_writer,
+            "tracers0",
+            "velocity",
+            psolver.getIteration());
+    pset.writeStateTriplet(
+            0,
+            particle_sample_writer,
+            "tracers0",
+            "position",
+            psolver.getIteration());
+    particle_output_writer->close_file();
 
     // deallocate
-    delete particles_output_writer_mpi;
+    delete particle_sample_writer;
+    delete particle_output_writer;
     delete trhs;
     delete fti;
     delete kk;
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index 4a54e6e6587e098d5dd26c198ef9cc53cc116165..2361b142fdaedbe1da65a496e3e2ef403733008c 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -159,7 +159,7 @@ class abstract_particle_set
         template <typename field_rnumber,
                   field_backend be,
                   field_components fc>
-        int sample(
+        int writeSample(
                 field<field_rnumber, be, fc> *field_to_sample,
                 particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                 const std::string species_name,
@@ -190,7 +190,7 @@ class abstract_particle_set
             return EXIT_SUCCESS;
         }
 
-        int sample_state_triple(
+        int writeStateTriplet(
                 const int i0,
                 particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
                 const std::string species_name,