diff --git a/src/Hipace.cpp b/src/Hipace.cpp
index 50f34474e7..79e06e3649 100644
--- a/src/Hipace.cpp
+++ b/src/Hipace.cpp
@@ -466,8 +466,8 @@ Hipace::Evolve ()
         // exit loop over time steps, if max time is exceeded
         if (m_physical_time > m_max_time) break;
 
-        m_multi_beam.InSituWriteToFile(step, m_physical_time, m_3D_geom[0]);
-        m_multi_plasma.InSituWriteToFile(step, m_physical_time, m_3D_geom[0]);
+        m_multi_beam.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
+        m_multi_plasma.InSituWriteToFile(step, m_physical_time, m_3D_geom[0], m_max_step, m_max_time);
 
         // printing and resetting predictor corrector loop diagnostics
         if (m_verbose>=2 && !m_explicit) amrex::AllPrint() << "Rank " << rank <<
@@ -501,8 +501,9 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step)
     // for function, the parallelcontext is the transverse communicator
     amrex::ParallelContext::push(m_comm_xy);
 
-    m_multi_beam.InSituComputeDiags(step, islice, islice_local);
-    m_multi_plasma.InSituComputeDiags(step, islice);
+    m_multi_beam.InSituComputeDiags(step, islice, islice_local, m_max_step,
+                                    m_physical_time, m_max_time);
+    m_multi_plasma.InSituComputeDiags(step, islice, m_max_step, m_physical_time, m_max_time);
 
     // Get this laser slice from the 3D array
     m_multi_laser.Copy(islice, false);
diff --git a/src/diagnostics/Diagnostic.H b/src/diagnostics/Diagnostic.H
index 9350613651..1e4559df0f 100644
--- a/src/diagnostics/Diagnostic.H
+++ b/src/diagnostics/Diagnostic.H
@@ -8,6 +8,8 @@
 #ifndef DIAGNOSTIC_H_
 #define DIAGNOSTIC_H_
 
+#include "utils/IOUtil.H"
+
 #include <AMReX_MultiFab.H>
 #include <AMReX_Vector.H>
 
@@ -74,10 +76,9 @@ public:
                          int output_step, int max_step,
                          amrex::Real output_time, amrex::Real max_time)
     {
-        return fd.m_output_period > 0 && fd.m_nfields_with_laser > 0 && (
-            (output_time == max_time) ||
-            (output_step == max_step) ||
-            (output_step % fd.m_output_period == 0) );
+        return  fd.m_nfields_with_laser > 0 &&
+            utils::doDiagnostics(fd.m_output_period, output_step, max_step,
+                            output_time, max_time);
     }
 
     /** \brief determines if any field diagnostic has any output on in this time step
@@ -106,10 +107,9 @@ public:
     bool hasBeamOutput (int output_step, int max_step,
                         amrex::Real output_time, amrex::Real max_time) const
     {
-        return m_beam_output_period > 0 && m_output_beam_names.size() > 0 && (
-            (output_time == max_time) ||
-            (output_step == max_step) ||
-            (output_step % m_beam_output_period == 0) );
+        return m_output_beam_names.size() > 0 &&
+            utils::doDiagnostics(m_beam_output_period, output_step, max_step,
+                            output_time, max_time);
     }
 
     /** \brief determines if any field or beam diagnostic has any output on in this time step
diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H
index 5558901324..0e87f192ed 100644
--- a/src/particles/beam/BeamParticleContainer.H
+++ b/src/particles/beam/BeamParticleContainer.H
@@ -185,9 +185,9 @@ public:
                      const int which_slice, const int islice_local, const int nghost);
 
     amrex::Long TotalNumberOfParticles (bool only_valid=true, bool only_local=false) const;
-
-    bool doInSitu (int step);
-
+    /** How often the insitu beam diagnostics should be computed and written
+     * Default is 0, meaning no output */
+    int m_insitu_period {0};
 private:
     std::string m_name; /**< name of the species */
     amrex::Real m_zmin; /**< Min longitudinal particle position of the beam */
@@ -249,9 +249,6 @@ private:
     amrex::Vector<int> m_insitu_sum_idata;
     /** Prefix/path for the output files */
     std::string m_insitu_file_prefix = "diags/insitu";
-    /** How often the insitu beam diagnostics should be computed and written
-     * Default is 0, meaning no output */
-    int m_insitu_period {0};
 };
 
 #endif
diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp
index 2d02c83a0d..745f5b6faf 100644
--- a/src/particles/beam/BeamParticleContainer.cpp
+++ b/src/particles/beam/BeamParticleContainer.cpp
@@ -304,12 +304,6 @@ void BeamParticleContainer::TagByLevel (const int current_N_level,
     );
 }
 
-bool
-BeamParticleContainer::doInSitu (int step)
-{
-    return (m_insitu_period > 0 && step % m_insitu_period == 0);
-}
-
 void
 BeamParticleContainer::InSituComputeDiags (int islice, int islice_local)
 {
diff --git a/src/particles/beam/MultiBeam.H b/src/particles/beam/MultiBeam.H
index f6f40ceba4..ebd1b5a4f9 100644
--- a/src/particles/beam/MultiBeam.H
+++ b/src/particles/beam/MultiBeam.H
@@ -83,9 +83,22 @@ public:
      * \param[in] step time step of simulation
      * \param[in] islice current slice, on which diags are computed.
      * \param[in] islice_local local index of the slice
+     * \param[in] max_step maximum time step of simulation
+     * \param[in] physical_time physical time at the given step
+     * \param[in] max_time maximum time of simulation
+     */
+    void InSituComputeDiags (int step, int islice, int islice_local,
+                             int max_step, amrex::Real physical_time,
+                             amrex::Real max_time);
+    /** Write reduced beam diagnostics to file
+     * \param[in] step time step of simulation
+     * \param[in] time physical time at the given step
+     * \param[in] geom Simulation geometry
+     * \param[in] max_step maximum time step of simulation
+     * \param[in] max_time maximum time of simulation
      */
-    void InSituComputeDiags (int step, int islice, int islice_local);
-    void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom);
+    void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom,
+                            int max_step, amrex::Real max_time);
     /** Loop over species and init them
      * \param[in] geom Simulation geometry
      * \return physical time at which the simulation will start
diff --git a/src/particles/beam/MultiBeam.cpp b/src/particles/beam/MultiBeam.cpp
index 2b523072b3..f25587b451 100644
--- a/src/particles/beam/MultiBeam.cpp
+++ b/src/particles/beam/MultiBeam.cpp
@@ -10,6 +10,7 @@
 #include "particles/sorting/SliceSort.H"
 #include "particles/pusher/BeamParticleAdvance.H"
 #include "utils/DeprecatedInput.H"
+#include "utils/IOUtil.H"
 #include "utils/HipaceProfilerWrapper.H"
 
 MultiBeam::MultiBeam ()
@@ -215,20 +216,25 @@ MultiBeam::PackLocalGhostParticles (int it)
 }
 
 void
-MultiBeam::InSituComputeDiags (int step, int islice, int islice_local)
+MultiBeam::InSituComputeDiags (int step, int islice, int islice_local,
+                               int max_step, amrex::Real physical_time,
+                               amrex::Real max_time)
 {
-    for (int i = 0; i < m_nbeams; ++i) {
-        if (m_all_beams[i].doInSitu(step)) {
-            m_all_beams[i].InSituComputeDiags(islice, islice_local);
+    for (auto& beam : m_all_beams) {
+        if (utils::doDiagnostics(beam.m_insitu_period, step,
+                            max_step, physical_time, max_time)) {
+            beam.InSituComputeDiags(islice, islice_local);
         }
     }
 }
 
 void
-MultiBeam::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom)
+MultiBeam::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom,
+                              int max_step, amrex::Real max_time)
 {
     for (auto& beam : m_all_beams) {
-        if (beam.doInSitu(step)) {
+        if (utils::doDiagnostics(beam.m_insitu_period, step,
+                            max_step, time, max_time)) {
             beam.InSituWriteToFile(step, time, geom);
         }
     }
diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H
index 791db0614e..2ac504ea80 100644
--- a/src/particles/plasma/MultiPlasma.H
+++ b/src/particles/plasma/MultiPlasma.H
@@ -148,9 +148,21 @@ public:
     /** Compute reduced plasma diagnostics of current slice, store in member variable.
      * \param[in] step time step of simulation
      * \param[in] islice current slice, on which diags are computed.
+     * \param[in] max_step maximum time step of simulation
+     * \param[in] physical_time physical time at the given step
+     * \param[in] max_time maximum time of simulation
      */
-    void InSituComputeDiags (int step, int islice);
-    void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom);
+    void InSituComputeDiags (int step, int islice, int max_step,
+                             amrex::Real physical_time, amrex::Real max_time);
+    /** Write reduced beam diagnostics to file
+     * \param[in] step time step of simulation
+     * \param[in] time physical time at the given step
+     * \param[in] geom Simulation geometry
+     * \param[in] max_step maximum time step of simulation
+     * \param[in] max_time maximum time of simulation
+     */
+    void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom,
+                            int max_step, amrex::Real max_time);
 
     int m_sort_bin_size {32}; /**< Tile size to sort plasma particles */
 
diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp
index 674d904a97..46233015a2 100644
--- a/src/particles/plasma/MultiPlasma.cpp
+++ b/src/particles/plasma/MultiPlasma.cpp
@@ -12,6 +12,7 @@
 #include "particles/sorting/TileSort.H"
 #include "utils/HipaceProfilerWrapper.H"
 #include "utils/DeprecatedInput.H"
+#include "utils/IOUtil.H"
 #include "Hipace.H"
 
 MultiPlasma::MultiPlasma ()
@@ -203,20 +204,24 @@ MultiPlasma::TagByLevel (const int current_N_level, amrex::Vector<amrex::Geometr
 }
 
 void
-MultiPlasma::InSituComputeDiags (int step, int islice)
+MultiPlasma::InSituComputeDiags (int step, int islice, int max_step,
+                                amrex::Real physical_time, amrex::Real max_time)
 {
     for (auto& plasma : m_all_plasmas) {
-        if (plasma.doInSitu(step)) {
+        if (utils::doDiagnostics(plasma.m_insitu_period, step,
+                            max_step, physical_time, max_time)) {
             plasma.InSituComputeDiags(islice);
         }
     }
 }
 
 void
-MultiPlasma::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom)
+MultiPlasma::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom,
+                                int max_step, amrex::Real max_time)
 {
     for (auto& plasma : m_all_plasmas) {
-        if (plasma.doInSitu(step)) {
+        if (utils::doDiagnostics(plasma.m_insitu_period, step,
+                            max_step, time, max_time)) {
             plasma.InSituWriteToFile(step, time, geom);
         }
     }
diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H
index 8ad760b55e..1f7e1af664 100644
--- a/src/particles/plasma/PlasmaParticleContainer.H
+++ b/src/particles/plasma/PlasmaParticleContainer.H
@@ -150,7 +150,6 @@ public:
      * \param[in] geom Geometry object for the whole domain
      */
     void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom);
-    bool doInSitu (int step);
 
     amrex::Parser m_parser; /**< owns data for m_density_func */
     amrex::ParserExecutor<3> m_density_func; /**< Density function for the plasma */
@@ -186,7 +185,9 @@ public:
     int m_reorder_period = 0;
     /** 2D reordering index type. 0: cell, 1: node, 2: both */
     amrex::IntVect m_reorder_idx_type = {0, 0, 0};
-
+    /** How often the insitu plasma diagnostics should be computed and written
+     * Default is 0, meaning no output */
+    int m_insitu_period {0};
 private:
     std::string m_name; /**< name of the species */
     int m_nslices; /**< number of z slices of the domain */
@@ -215,9 +216,6 @@ private:
     amrex::Vector<int> m_insitu_sum_idata;
     /** Prefix/path for the output files */
     std::string m_insitu_file_prefix = "diags/plasma_insitu";
-    /** How often the insitu plasma diagnostics should be computed and written
-     * Default is 0, meaning no output */
-    int m_insitu_period {0};
 };
 
 /** \brief Iterator over boxes in a particle container */
diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp
index c3e159827f..081386500f 100644
--- a/src/particles/plasma/PlasmaParticleContainer.cpp
+++ b/src/particles/plasma/PlasmaParticleContainer.cpp
@@ -431,12 +431,6 @@ IonizationModule (const int lev,
     }
 }
 
-bool
-PlasmaParticleContainer::doInSitu (int step)
-{
-    return (m_insitu_period > 0 && step % m_insitu_period == 0);
-}
-
 void
 PlasmaParticleContainer::InSituComputeDiags (int islice)
 {
diff --git a/src/utils/IOUtil.H b/src/utils/IOUtil.H
index 7c32bca1ab..fdfd93c82a 100644
--- a/src/utils/IOUtil.H
+++ b/src/utils/IOUtil.H
@@ -39,7 +39,7 @@ namespace utils
      * (used for compatibility with the openPMD API)
      */
     std::vector<std::uint64_t>
-    getReversedVec( const amrex::IntVect& v );
+    getReversedVec ( const amrex::IntVect& v );
 
     /** \brief
      * Convert Real* pointer to a std::vector<double>,
@@ -47,11 +47,22 @@ namespace utils
      * (used for compatibility with the openPMD API)
      */
     std::vector<double>
-    getReversedVec( const amrex::Real* v );
+    getReversedVec ( const amrex::Real* v );
+
+    /** \brief
+     * returns whether output should be writen to file
+     * \param[in] output_period period of the output
+     * \param[in] output_step current step
+     * \param[in] max_step maximum step of simulation
+     * \param[in] output_time physical time of the current step
+     * \param[in] max_time maximum simulation time
+     */
+     bool doDiagnostics (int output_period, int output_step, int max_step,
+                    amrex::Real output_time, amrex::Real max_time);
 
 #ifdef HIPACE_USE_OPENPMD
     std::pair< std::string, std::string >
-    name2openPMD( std::string const& fullName );
+    name2openPMD ( std::string const& fullName );
 
     /** Get the openPMD physical dimensionality of a record
      *
@@ -59,7 +70,7 @@ namespace utils
      * @return map with base quantities and power scaling
      */
     std::map< openPMD::UnitDimension, double >
-    getUnitDimension( std::string const & record_name );
+    getUnitDimension ( std::string const & record_name );
 #endif
 
 }
diff --git a/src/utils/IOUtil.cpp b/src/utils/IOUtil.cpp
index db8c145028..8906b7bd3a 100644
--- a/src/utils/IOUtil.cpp
+++ b/src/utils/IOUtil.cpp
@@ -14,7 +14,7 @@
 
 
 std::vector< double >
-utils::getRelativeCellPosition(amrex::FArrayBox const& fab)
+utils::getRelativeCellPosition (amrex::FArrayBox const& fab)
 {
     amrex::Box const box = fab.box();
     amrex::IndexType const idx_type = box.ixType();
@@ -31,7 +31,7 @@ utils::getRelativeCellPosition(amrex::FArrayBox const& fab)
 }
 
 std::vector<std::uint64_t>
-utils::getReversedVec( const amrex::IntVect& v )
+utils::getReversedVec ( const amrex::IntVect& v )
 {
   // Convert the IntVect v to and std::vector u
   std::vector<std::uint64_t> u = {
@@ -49,7 +49,7 @@ utils::getReversedVec( const amrex::IntVect& v )
 }
 
 std::vector<double>
-utils::getReversedVec( const amrex::Real* v )
+utils::getReversedVec ( const amrex::Real* v )
 {
   // Convert Real* v to and std::vector u
   std::vector<double> u = {
@@ -66,9 +66,19 @@ utils::getReversedVec( const amrex::Real* v )
   return u;
 }
 
+bool
+utils::doDiagnostics (int output_period, int output_step, int max_step,
+               amrex::Real output_time, amrex::Real max_time)
+{
+    return output_period > 0 && (
+        (output_time == max_time) ||
+        (output_step == max_step) ||
+        (output_step % output_period == 0) );
+}
+
 #ifdef HIPACE_USE_OPENPMD
 std::pair< std::string, std::string >
-utils::name2openPMD( std::string const& fullName )
+utils::name2openPMD ( std::string const& fullName )
 {
     std::string record_name = fullName;
     std::string component_name = openPMD::RecordComponent::SCALAR;
@@ -87,7 +97,7 @@ utils::name2openPMD( std::string const& fullName )
  * @return map with base quantities and power scaling
  */
 std::map< openPMD::UnitDimension, double >
-utils::getUnitDimension( std::string const & record_name )
+utils::getUnitDimension ( std::string const & record_name )
 {
 
     if( record_name == "position" ) return {