From 3ff473cca03cd7657a94ae5c44df57e406df21b5 Mon Sep 17 00:00:00 2001
From: chrisrogers1234 <chrisrogers1234@yahoo.co.uk>
Date: Mon, 20 Mar 2017 17:36:32 +0000
Subject: [PATCH] Fixes #56 - OpalRing components throw segmentation fault and
 #59 - Unit tests throw segmentation fault; partial implementation of #68
 DumpEmFields

---
 src/Algorithms/ParallelCyclotronTracker.cpp   |   4 +-
 src/BasicActions/CMakeLists.txt               |   4 +-
 src/BasicActions/DumpEMFields.cpp             | 226 +++++++++++
 src/BasicActions/DumpEMFields.h               | 117 ++++++
 src/BasicActions/DumpFields.h                 |   6 +-
 src/Classic/AbsBeamline/Component.cpp         |   2 +-
 src/Classic/AbsBeamline/Ring.cpp              |   2 +
 src/Classic/AbsBeamline/SBend3D.cpp           |   4 +-
 .../Fields/Interpolation/CMakeLists.txt       |   4 +-
 src/Classic/Fields/Interpolation/NDGrid.cpp   | 204 ++++++++++
 src/Classic/Fields/Interpolation/NDGrid.h     | 383 ++++++++++++++++++
 .../Fields/Interpolation/PPSolveFactory.cpp   |   6 +-
 src/Classic/Fields/SectorMagneticFieldMap.cpp |  22 +-
 src/Classic/Utilities/RingSection.cpp         |   5 +-
 src/OpalConfigure/Configure.cpp               |   3 +-
 tests/CMakeLists.txt                          |   4 +-
 tests/Main.cpp                                |  10 +-
 tests/classic_src/AbsBeamline/RingTest.cpp    | 198 +++------
 tests/classic_src/AbsBeamline/SBend3DTest.cpp |   4 -
 .../Fields/Interpolation/NDGridTest.cpp       |  36 ++
 .../Fields/Interpolation/ThreeDGridTest.cpp   |   3 +
 tests/opal_src/BasicActions/CMakeLists.txt    |   1 +
 .../BasicActions/DumpEMFieldsTest.cpp         | 180 ++++++++
 tests/opal_src/Distribution/GaussTest.cpp     |  65 +--
 tests/opal_test_utilities/CMakeLists.txt      |   8 +
 tests/opal_test_utilities/SilenceTest.h       |  56 +++
 26 files changed, 1338 insertions(+), 219 deletions(-)
 create mode 100644 src/BasicActions/DumpEMFields.cpp
 create mode 100644 src/BasicActions/DumpEMFields.h
 create mode 100644 src/Classic/Fields/Interpolation/NDGrid.cpp
 create mode 100644 src/Classic/Fields/Interpolation/NDGrid.h
 create mode 100644 tests/classic_src/Fields/Interpolation/NDGridTest.cpp
 create mode 100644 tests/opal_src/BasicActions/DumpEMFieldsTest.cpp
 create mode 100644 tests/opal_test_utilities/CMakeLists.txt
 create mode 100644 tests/opal_test_utilities/SilenceTest.h

diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 366159ca6..4d851fb9a 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -73,6 +73,7 @@
 #include "Utilities/OpalException.h"
 
 #include "BasicActions/DumpFields.h"
+#include "BasicActions/DumpEMFields.h"
 
 #include "Structure/H5PartWrapperForPC.h"
 #include "Structure/BoundaryGeometry.h"
@@ -1145,6 +1146,7 @@ void ParallelCyclotronTracker::execute() {
     extE_m = Vector_t(0.0, 0.0, 0.0);
     extB_m = Vector_t(0.0, 0.0, 0.0);
     DumpFields::writeFields((((*FieldDimensions.begin())->second).second));
+    DumpEMFields::writeFields((((*FieldDimensions.begin())->second).second));
 
     if(timeIntegrator_m == 0) {
         *gmsg << "* 4th order Runge-Kutta integrator" << endl;
@@ -5573,4 +5575,4 @@ void ParallelCyclotronTracker::evaluateSpaceChargeField() {
         localToGlobal(itsBunch->Bf, phi);
         localToGlobal(itsBunch->R, phi, meanR);
     }
-}
\ No newline at end of file
+}
diff --git a/src/BasicActions/CMakeLists.txt b/src/BasicActions/CMakeLists.txt
index 4d6e8d38f..aa7b58850 100644
--- a/src/BasicActions/CMakeLists.txt
+++ b/src/BasicActions/CMakeLists.txt
@@ -2,6 +2,7 @@ set (_SRCS
   Call.cpp
   Dump.cpp
   DumpFields.cpp
+  DumpEMFields.cpp
   Echo.cpp
   Help.cpp
   Option.cpp
@@ -26,6 +27,7 @@ add_opal_sources(${_SRCS})
 set (HDRS
     Call.h
     DumpFields.h
+    DumpEMFields.h
     Dump.h
     Echo.h
     Help.h
@@ -42,4 +44,4 @@ set (HDRS
     What.h
 )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/BasicActions")
\ No newline at end of file
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/BasicActions")
diff --git a/src/BasicActions/DumpEMFields.cpp b/src/BasicActions/DumpEMFields.cpp
new file mode 100644
index 000000000..9d0c97583
--- /dev/null
+++ b/src/BasicActions/DumpEMFields.cpp
@@ -0,0 +1,226 @@
+/*
+ *  Copyright (c) 2016, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are met:
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *  2. Redistributions in binary form must reproduce the above copyright notice,
+ *     this list of conditions and the following disclaimer in the documentation
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to
+ *     endorse or promote products derived from this software without specific
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+
+#include <fstream>
+
+#include "Fields/Interpolation/NDGrid.h"  // classic
+#include "AbsBeamline/Component.h"  // classic
+#include "Utilities/OpalException.h"
+#include "Attributes/Attributes.h"
+#include "BasicActions/DumpEMFields.h"
+
+// extern Inform *gmsg;
+
+std::unordered_set<DumpEMFields*> DumpEMFields::dumpsSet_m;
+
+std::string DumpEMFields::dumpemfields_docstring =
+std::string("The \"DUMPEMFIELDS\" statement dumps a field map to a user-defined")+
+std::string(" field file, for checking that fields are generated correctly.")+
+std::string(" The fields are written out on a grid in Cartesian space and time.");
+
+DumpEMFields::DumpEMFields() :
+                    Action(13, "DUMPEMFIELDS", dumpemfields_docstring.c_str()),
+                    grid_m(NULL), filename_m("") {
+    // would be nice if "steps" could be integer
+    itsAttr[0] = Attributes::makeReal
+                 ("X_START", "Start point in the grid in x [mm]");
+    itsAttr[1] = Attributes::makeReal
+                 ("DX", "Grid step size in x [mm]");
+    itsAttr[2] = Attributes::makeReal
+                 ("X_STEPS", "Number of steps in x");
+    itsAttr[3] = Attributes::makeReal
+                 ("Y_START", "Start point in the grid in y [mm]");
+    itsAttr[4] = Attributes::makeReal
+                 ("DY", "Grid step size in y [mm]");
+    itsAttr[5] = Attributes::makeReal
+                 ("Y_STEPS", "Number of steps in y");
+    itsAttr[6] = Attributes::makeReal
+                 ("Z_START", "Start point in the grid in z [mm]");
+    itsAttr[7] = Attributes::makeReal
+                 ("DZ", "Grid step size in z [mm]");
+    itsAttr[8] = Attributes::makeReal
+                 ("Z_STEPS", "Number of steps in z");
+    itsAttr[9] = Attributes::makeReal
+                 ("T_START", "Start point in the grid in time [ns]");
+    itsAttr[10] = Attributes::makeReal
+                 ("DT", "Grid step size in time [ns]");
+    itsAttr[11] = Attributes::makeReal
+                 ("T_STEPS", "Number of steps in time");
+    itsAttr[12] = Attributes::makeString
+                 ("FILE_NAME", "Name of the file to which field data is dumped");
+}
+
+DumpEMFields::~DumpEMFields() {
+    if (grid_m != NULL) {
+        delete grid_m;
+    }
+    if (dumpsSet_m.find(this) != dumpsSet_m.end()) {
+        dumpsSet_m.erase(this);
+    }
+}
+
+DumpEMFields* DumpEMFields::clone(const std::string &name) {
+    DumpEMFields* dumper = new DumpEMFields();
+    if (grid_m != NULL) {
+        dumper->grid_m = grid_m->clone();
+    }
+    dumper->filename_m = filename_m;
+    if (dumpsSet_m.find(this) != dumpsSet_m.end()) {
+        dumpsSet_m.insert(dumper);
+    }
+    return dumper;
+}
+
+void DumpEMFields::execute() {
+    buildGrid();
+    // the routine for action (OpalParser/OpalParser) calls execute and then
+    // deletes 'this'; so we must build a copy that lasts until the field maps
+    // are constructed and we are ready for tracking (which is when the field
+    // maps are written). Hence the clone call below.
+    dumpsSet_m.insert(this->clone(""));
+}
+
+void DumpEMFields::buildGrid() {
+    std::vector<double> spacing(4); 
+    std::vector<double> origin(4); 
+    std::vector<int> gridSize(4); 
+    origin[0] = Attributes::getReal(itsAttr[0]);
+    spacing[0] = Attributes::getReal(itsAttr[1]);
+    double nx = Attributes::getReal(itsAttr[2]);
+
+    origin[1] = Attributes::getReal(itsAttr[3]);
+    spacing[1] = Attributes::getReal(itsAttr[4]);
+    double ny = Attributes::getReal(itsAttr[5]);
+
+    origin[2] = Attributes::getReal(itsAttr[6]);
+    spacing[2] = Attributes::getReal(itsAttr[7]);
+    double nz = Attributes::getReal(itsAttr[8]);
+
+    origin[3] = Attributes::getReal(itsAttr[9]);
+    spacing[3] = Attributes::getReal(itsAttr[10]);
+    double nt = Attributes::getReal(itsAttr[11]);
+
+    checkInt(nx, "X_STEPS");
+    checkInt(ny, "Y_STEPS");
+    checkInt(nz, "Z_STEPS");
+    checkInt(nt, "T_STEPS");
+
+    gridSize[0] = nx;
+    gridSize[1] = ny;
+    gridSize[2] = nz;
+    gridSize[3] = nt;
+
+    if (grid_m != NULL) {
+        delete grid_m;
+    }
+
+    grid_m = new interpolation::NDGrid(4, &gridSize[0], &spacing[0], &origin[0]);
+
+    filename_m = Attributes::getString(itsAttr[12]);
+}
+
+void DumpEMFields::writeFields(Component* field) {
+    typedef std::unordered_set<DumpEMFields*>::iterator dump_iter;
+    for (dump_iter it = dumpsSet_m.begin(); it != dumpsSet_m.end(); ++it) {
+        (*it)->writeFieldThis(field);
+    }
+}
+
+void DumpEMFields::checkInt(double real, std::string name, double tolerance) {
+    if (fabs(floor(real) - real) > tolerance) {
+        throw OpalException("DumpEMFields::checkInt",
+                            "Value for "+name+
+                            " should be an integer but a real value was found");
+    }
+    if (floor(real) < 0.5) {
+        throw OpalException("DumpEMFields::checkInt",
+                            "Value for "+name+" should be 1 or more");
+    }
+}
+
+void DumpEMFields::writeFieldThis(Component* field) {
+    if (grid_m == NULL) {
+        throw OpalException("DumpEMFields::writeFieldThis",
+                            "The grid was NULL; there was a problem with the DumpEMFields initialisation.");
+    }
+
+    if (field == NULL) {
+        throw OpalException("DumpEMFields::writeFieldThis",
+                            "The field to be written was NULL.");
+    }
+    double time = 0.;
+    std::vector<double> point_std(4);
+    Vector_t point(0., 0., 0.);
+    Vector_t centroid(0., 0., 0.);
+    std::ofstream fout;
+    try {
+        fout.open(filename_m.c_str(), std::ofstream::out);
+    } catch (std::exception& exc) {
+        throw OpalException("DumpEMFields::writeFieldThis",
+                            "Failed to open DumpEMFields file "+filename_m);
+    }
+    if (!fout.good()) {
+        throw OpalException("DumpEMFields::writeFieldThis",
+                            "Failed to open DumpEMFields file "+filename_m);
+    }
+    // set precision
+    fout << grid_m->end().toInteger() << "\n";
+    fout << 1 << " x [mm]\n";
+    fout << 2 << " y [mm]\n";
+    fout << 3 << " z [mm]\n";
+    fout << 4 << " t [ns]\n";
+    fout << 5 << " Bx [kGauss]\n";
+    fout << 6 << " By [kGauss]\n";
+    fout << 7 << " Bz [kGauss]\n";
+    fout << 8 << " Ex [?]\n";
+    fout << 9 << " Ey [?]\n";
+    fout << 10 << " Ez [?]\n";
+    fout << 0 << std::endl;
+    for (interpolation::Mesh::Iterator it = grid_m->begin();
+         it < grid_m->end();
+         ++it) {
+        Vector_t E(0., 0., 0.);
+        Vector_t B(0., 0., 0.);
+        it.getPosition(&point_std[0]);
+        for (size_t i = 0; i < 3; ++i) {
+            point[i] = point_std[i];
+        }
+        time = point_std[3];
+
+        field->apply(point, centroid, time, E, B);
+        fout << point[0] << " " << point[1] << " " << point[2] << " " << time << " ";
+        fout << B[0] << " " << B[1] << " " << B[2] << " ";
+        fout << E[0] << " " << E[1] << " " << E[2] << "\n";
+    }
+    if (!fout.good()) {
+        throw OpalException("DumpEMFields::writeFieldThis",
+                            "Something went wrong during writing "+filename_m);
+    }
+    fout.close();
+}
+
diff --git a/src/BasicActions/DumpEMFields.h b/src/BasicActions/DumpEMFields.h
new file mode 100644
index 000000000..b26b8dfb9
--- /dev/null
+++ b/src/BasicActions/DumpEMFields.h
@@ -0,0 +1,117 @@
+/*
+ *  Copyright (c) 2017, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are met:
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *  2. Redistributions in binary form must reproduce the above copyright notice,
+ *     this list of conditions and the following disclaimer in the documentation
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to
+ *     endorse or promote products derived from this software without specific
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef OPAL_BASICACTIONS_DUMPEMFIELDS_HH
+#define OPAL_BASICACTIONS_DUMPEMFIELDS_HH
+
+#include <unordered_set>
+#include <string>
+#include "AbstractObjects/Action.h"
+
+namespace interpolation {
+  class NDGrid;
+}
+
+class Component;
+
+/** DumpEMFields dumps the dynamically changing fields of a Ring in a user-
+ *  defined grid.
+ *
+ *  The idea is to print out the field map across a 4D grid in space-time for
+ *  debugging purposes. The problem is to manage the DumpEMFields object through
+ *  three phases of program execution; initial construction, parsing and then
+ *  actual field map writing (where we need to somehow let DumpFields know what
+ *  the field maps are). So for each  DumpFields object created, we store in a
+ *  set. When the execute() method is called, DumpFields builds a grid using
+ *  the parsed information.
+ *
+ *  When the ParallelCyclotronTracker is about to start tracking, it calls
+ *  writeFields method which loops over the static set of DumpFields and writes
+ *  each one. It is not the cleanest implementation, but I can't see a better
+ *  way.
+ *
+ *  The DumpEMFields themselves operate by iterating over a NDGrid object
+ *  and looking up the field/writing it out on each grid point.
+ *
+ */
+class DumpEMFields : public Action {
+  public:
+    /** Constructor */
+    DumpEMFields();
+
+    /** Destructor deletes grid_m and if in the dumps set, take it out */
+    virtual ~DumpEMFields();
+
+    /** Make a clone (overloadable copy-constructor).
+     *    @param name not used
+     *  If this is in the dumpsSet_m, so will the clone. Not sure how the
+     *  itsAttr stuff works, so this may not get properly copied?
+     */
+    virtual DumpEMFields *clone(const std::string &name);
+
+    /** Builds the grid but does not write the field map
+     *
+     *  Builds a grid of points in x-y-z space using the NDGrid algorithm.
+     *  Checks that X_STEPS, Y_STEPS, Z_STEPS are integers or throws
+     *  OpalException.
+     */
+    virtual void execute();
+
+    /** Write the fields for all defined DumpEMFields objects
+     *    @param field borrowed reference to the Component object that holds the
+     *    field map; caller owns the memory.
+     *  Iterates over the DumpEMFields in the dumpsSet_m and calls writeFieldThis
+     *  on each DumpEMFields. This writes each field map in turn. Format is:
+     *  <number of rows>
+     *  <column 1> <units>
+     *  <column 2> <units>
+     *  <column 3> <units>
+     *  <column 4> <units>
+     *  <column 5> <units>
+     *  <column 6> <units>
+     *  0
+     *  <field map data>
+     */
+    static void writeFields(Component* field);
+
+  private:
+    virtual void writeFieldThis(Component* field);
+    virtual void buildGrid();
+    static void checkInt(double value, std::string name, double tolerance = 1e-9);
+
+    interpolation::NDGrid* grid_m;
+    std::string filename_m;
+
+    static std::unordered_set<DumpEMFields*> dumpsSet_m;
+    static std::string dumpemfields_docstring;
+
+    DumpEMFields(const DumpEMFields& dump);  // disabled
+    DumpEMFields& operator=(const DumpEMFields& dump);  // disabled
+};
+
+#endif // ifdef OPAL_DUMPFIELDS_HH
+
diff --git a/src/BasicActions/DumpFields.h b/src/BasicActions/DumpFields.h
index 739f482c5..54e47d600 100644
--- a/src/BasicActions/DumpFields.h
+++ b/src/BasicActions/DumpFields.h
@@ -38,9 +38,9 @@ namespace interpolation {
 
 class Component;
 
-/** DumpFields dumps the fields of a Ring in a user-defined grid
+/** DumpFields dumps the static magnetic field of a Ring in a user-defined grid
  *
- *  The idea is to print out the field map across a grid in space-time for
+ *  The idea is to print out the field map across a grid in space for
  *  debugging purposes. The problem is to manage the DumpFields object through
  *  three phases of program execution; initial construction, parsing and then
  *  actual field map writing (where we need to somehow let DumpFields know what
@@ -56,6 +56,8 @@ class Component;
  *  The DumpFields themselves operate by iterating over a ThreeDGrid object
  *  and looking up the field/writing it out on each grid point.
  *
+ *  In order to dump time dependent fields, for example RF, see the DumpEMFields
+ *  action.
  */
 class DumpFields : public Action {
   public:
diff --git a/src/Classic/AbsBeamline/Component.cpp b/src/Classic/AbsBeamline/Component.cpp
index 14b1d1aaf..94713338c 100644
--- a/src/Classic/AbsBeamline/Component.cpp
+++ b/src/Classic/AbsBeamline/Component.cpp
@@ -17,7 +17,7 @@
 //
 // $Date: 2000/03/27 09:32:31 $
 // $Author: fci $
-//
+//c
 // ------------------------------------------------------------------------
 
 #include "AbsBeamline/Component.h"
diff --git a/src/Classic/AbsBeamline/Ring.cpp b/src/Classic/AbsBeamline/Ring.cpp
index aaefd7506..1edafb35f 100644
--- a/src/Classic/AbsBeamline/Ring.cpp
+++ b/src/Classic/AbsBeamline/Ring.cpp
@@ -105,8 +105,10 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
                 << " out of the field map boundary" << endl;
         lossDS_m->addParticle(refPartBunch_m->R[id], refPartBunch_m->P[id], id);
         lossDS_m->save();
+
         refPartBunch_m->Bin[id] = -1;
     }
+
     return flagNeedUpdate;
 }
 
diff --git a/src/Classic/AbsBeamline/SBend3D.cpp b/src/Classic/AbsBeamline/SBend3D.cpp
index 4eb7ec876..83b563a2a 100644
--- a/src/Classic/AbsBeamline/SBend3D.cpp
+++ b/src/Classic/AbsBeamline/SBend3D.cpp
@@ -33,6 +33,7 @@ SBend3D::SBend3D(const std::string &name)
         : Component(name), map_m(NULL),
          planarArcGeometry_m(1., 1.), fieldUnits_m(1.), lengthUnits_m(1.),
          polyOrder_m(1), smoothOrder_m(1), dummy() {
+    setElType(isDrift);
 }
 
 SBend3D::SBend3D(const SBend3D &right)
@@ -44,6 +45,7 @@ SBend3D::SBend3D(const SBend3D &right)
     RefPartBunch_m = right.RefPartBunch_m;
     if (right.map_m != NULL)
         map_m = new SectorMagneticFieldMap(*right.map_m);
+    setElType(isDrift);
 }
 
 SBend3D::~SBend3D() {
@@ -117,4 +119,4 @@ void SBend3D::setFieldMapFileName(std::string name) {
 
 void SBend3D::accept(BeamlineVisitor& visitor) const {
     visitor.visitSBend3D(*this);
-}
\ No newline at end of file
+}
diff --git a/src/Classic/Fields/Interpolation/CMakeLists.txt b/src/Classic/Fields/Interpolation/CMakeLists.txt
index 7cf76fb82..abdbb1249 100644
--- a/src/Classic/Fields/Interpolation/CMakeLists.txt
+++ b/src/Classic/Fields/Interpolation/CMakeLists.txt
@@ -4,6 +4,7 @@ set (_SRCS
   Mesh.cpp
   MMatrix.cpp
   MVector.cpp
+  NDGrid.cpp
   PolynomialPatch.cpp
   PPSolveFactory.cpp
   SolveFactory.cpp
@@ -25,6 +26,7 @@ set (HDRS
     Mesh.h
     MMatrix.h
     MVector.h
+    NDGrid.cpp
     PolynomialCoefficient.h
     PolynomialPatch.h
     PPSolveFactory.h
@@ -35,4 +37,4 @@ set (HDRS
     VectorMap.h
 )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Interpolation")
\ No newline at end of file
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/Interpolation")
diff --git a/src/Classic/Fields/Interpolation/NDGrid.cpp b/src/Classic/Fields/Interpolation/NDGrid.cpp
new file mode 100644
index 000000000..2863371f0
--- /dev/null
+++ b/src/Classic/Fields/Interpolation/NDGrid.cpp
@@ -0,0 +1,204 @@
+// MAUS WARNING: THIS IS LEGACY CODE.
+#include <math.h>
+#include <iomanip>
+#include "Utilities/OpalException.h"
+#include "Fields/Interpolation/Mesh.h"
+#include "Fields/Interpolation/NDGrid.h"
+
+namespace interpolation {
+/// ////// NDGrid ///////
+NDGrid::NDGrid() : coord_m(), maps_m(), constantSpacing_m(false)  {
+}
+
+
+
+NDGrid::NDGrid(std::vector<int> size, std::vector<const double *> gridCoordinates)
+    : coord_m(), maps_m(), constantSpacing_m(false) { 
+    for (unsigned int i=0; i<size.size(); i++) {
+        if (size[i] < 1) {
+            throw(OpalException("NDGrid::NDGrid(...)",
+                                "ND Grid must be at least 1x1x...x1 grid"));
+        }
+        coord_m.push_back(std::vector<double>(gridCoordinates[i],
+                                              gridCoordinates[i] + size[i])); 
+    }
+    setConstantSpacing();
+}
+
+NDGrid::NDGrid(int nDimensions, int* size, double* spacing, double* min)
+    : coord_m(nDimensions), maps_m(), constantSpacing_m(true) {
+    for (int i=0; i<nDimensions; i++) {
+        if (size[i] < 1) {
+            throw(OpalException("NDGrid::NDGrid(...)",
+                                "ND Grid must be at least 1x1x...x1 grid"));
+        }
+        coord_m[i] = std::vector<double>(size[i]);
+        for (unsigned int j=0; j<coord_m[i].size(); j++) {
+            coord_m[i][j] = min[i] + j*spacing[i];
+        }
+    }
+}
+
+NDGrid::NDGrid(std::vector< std::vector<double> > gridCoordinates)
+    : coord_m(gridCoordinates), maps_m(), constantSpacing_m(false) {
+    for (unsigned int i=0; i<gridCoordinates.size(); i++) {
+        if (gridCoordinates[i].size() < 1) {
+            throw (OpalException("NDGrid::NDGrid(...)",
+                                "ND Grid must be at least 1x1x...x1 grid"));
+        }
+    }
+    setConstantSpacing();
+}
+
+
+double* NDGrid::newCoordArray  ( const int& dimension)  const {
+    double * array = new double[coord_m[dimension].size() ];
+    for (unsigned int i=0; i<coord_m[dimension].size(); i++) {
+        array[i] = coord_m[dimension][i];
+    }
+    return array;
+}
+
+//Mesh::Iterator wraps around a std::vector<int>
+//it[0] is least significant, it[max] is most signifcant
+Mesh::Iterator& NDGrid::addEquals(Mesh::Iterator& lhs, int difference) const {
+    if (difference < 0) {
+        return subEquals(lhs, -difference);
+    }
+    std::vector<int> index(coord_m.size(),0);
+    std::vector<int> content(coord_m.size(),1);
+    for (int i = int(index.size()-2); i >= 0; i--) {
+        content[i] = content[i+1]*coord_m[i+1].size(); //content could be member variable
+    }
+    for (int i = 0; i < int(index.size()); i++) {
+        index[i] = difference/content[i];
+        difference -= index[i] * content[i];
+    }
+    for (unsigned int i=0; i<index.size(); i++) {
+        lhs.state_m[i] += index[i];
+    }
+    for (int i = int(index.size())-1; i > 0; i--) {
+        if (lhs.state_m[i] > int(coord_m[i].size())) {
+            lhs.state_m[i-1]++;
+            lhs.state_m[i] -= coord_m[i].size();
+        }
+    }
+
+    return lhs;
+}
+
+Mesh::Iterator& NDGrid::subEquals(Mesh::Iterator& lhs, int difference) const {
+    if (difference < 0) {
+        return addEquals(lhs, -difference);
+    }
+    std::vector<int> index(coord_m.size(),0);
+    std::vector<int> content(coord_m.size(),1);
+    for (int i = int(index.size()-2); i >= 0; i--) {
+        content[i] = content[i+1]*coord_m[i+1].size(); //content could be member variable
+    }
+    for (int i = 0; i < int(index.size()); i++) {
+        index[i]    = difference/content[i];
+        difference -= index[i] * content[i];
+    }
+    for (unsigned int i = 0; i < index.size(); i++) {
+        lhs.state_m[i] -= index[i];
+    }
+    for (int i=int(index.size())-1; i>0; i--) {
+        if (lhs.state_m[i] < 1) {
+            lhs.state_m[i-1]--;
+            lhs.state_m[i] += coord_m[i].size();
+        }
+    }
+    return lhs;
+}
+
+Mesh::Iterator& NDGrid::addEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
+    return addEquals(lhs, rhs.toInteger());
+}
+
+Mesh::Iterator& NDGrid::subEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
+    return subEquals(lhs, rhs.toInteger());
+}
+
+Mesh::Iterator& NDGrid::addOne(Mesh::Iterator& lhs) const {
+    int i = coord_m.size()-1;
+    while (lhs[i] == int(coord_m[i].size()) && i>0) {
+        lhs[i]=1; i--;
+    }
+    lhs[i]++;
+    return lhs;
+}
+
+Mesh::Iterator& NDGrid::subOne(Mesh::Iterator& lhs) const {
+    lhs[coord_m.size()-1] -= 1;
+
+    int i = coord_m.size()-1;
+    while (lhs[i] == 0 && i>0) {
+        lhs.state_m[i] = coord_m[i].size();
+        i--;
+        lhs.state_m[i]--;
+    }
+    return lhs;
+}
+
+void NDGrid::setConstantSpacing() {
+    constantSpacing_m = true;
+    for (unsigned int i = 0; i < coord_m.size(); i++) {
+        for (unsigned int j = 0; j < coord_m[i].size()-1; j++) {
+            double coord_j1 = coord_m[i][j+1];
+            double coord_j0 = coord_m[i][j];
+            double coord_1 = coord_m[i][1];
+            double coord_0 = coord_m[i][0];
+            if (fabs(1-(coord_j1-coord_j0)/(coord_1-coord_0)) > tolerance_m) {
+                constantSpacing_m = false;
+                return;
+            }
+        }
+    }
+}
+
+bool NDGrid::isGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const {
+    unsigned int i = 0;
+    // if all equal; rhs[i] = rhs.last
+    while (lhs.state_m[i] == rhs.state_m[i] && i < rhs.state_m.size()-1) {
+        i++; 
+    }
+    return (lhs[i] > rhs[i]);
+}
+
+int NDGrid::toInteger(const Mesh::Iterator& lhs) const {
+    int difference = 0;
+    std::vector<int> index  (coord_m.size(),0);
+    std::vector<int> content(coord_m.size(),1);
+    for (int i = int(index.size()-2); i >= 0; i--) {
+        content[i]  = content[i+1]*coord_m[i+1].size();
+    }
+    for (int i = 0; i < int(index.size()); i++) {
+        difference += (lhs.state_m[i]-1) * (content[i]);
+    }
+    return difference;
+}
+
+Mesh::Iterator NDGrid::getNearest(const double* position) const {
+    std::vector<int>    index(coord_m.size());
+    std::vector<double> pos(position, position+coord_m.size());
+    lowerBound(pos, index);
+    for (unsigned int i = 0; i < coord_m.size(); i++)
+    {
+        if (index[i] < int(coord_m[i].size()-1) && index[i] >= 0) {
+            index[i] += (2*(position[i] - coord_m[i][index[i]])  > coord_m[i][index[i]+1]-coord_m[i][index[i]] ? 2 : 1);
+        } else {
+            index[i]++;
+        }
+        if (index[i] < 1) {
+            index[i] = 1; 
+        }
+        if (index[i] > int(coord_m[i].size())) {
+            index[i] = coord_m[i].size(); 
+        }
+    }
+    return Mesh::Iterator(index, this);
+}
+
+////// NDGrid END ////////
+}
diff --git a/src/Classic/Fields/Interpolation/NDGrid.h b/src/Classic/Fields/Interpolation/NDGrid.h
new file mode 100644
index 000000000..293fb95b1
--- /dev/null
+++ b/src/Classic/Fields/Interpolation/NDGrid.h
@@ -0,0 +1,383 @@
+/* 
+ *  Copyright (c) 2017, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without 
+ *  modification, are permitted provided that the following conditions are met: 
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer. 
+ *  2. Redistributions in binary form must reproduce the above copyright notice, 
+ *     this list of conditions and the following disclaimer in the documentation 
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to 
+ *     endorse or promote products derived from this software without specific 
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef _CLASSIC_FIELDS_NDGRID_HH_
+#define _CLASSIC_FIELDS_NDGRID_HH_
+
+#include <vector>
+#include <algorithm>
+#include <math.h>
+#include <ostream>
+
+#include "Fields/Interpolation/Mesh.h"
+
+namespace interpolation {
+
+/** NDGrid holds grid information for a rectilinear grid in some arbitrary
+ *  dimensional space.
+ *
+ *  NDGrid holds the necessary information for regular and irregular grid data
+ *  for interpolation routines (but it all has to be rectilinear).
+ *
+ *  Where an irregular grid is used, we use std::lower_bound to find mapping
+ *  from a generic point to the grid points.
+ *
+ *  NDGrid holds a list of pointers to VectorMaps that use the grid for
+ *  interpolation and functions to Add and Remove VectorMaps. If the Remove
+ *  function removes the last VectorMap, then the ThreeDGrid is deleted. It's a
+ *  sort of smart pointer.
+ *
+ *  In general NDGrid is not bound checked.
+ */
+
+class NDGrid : public Mesh {
+  public:
+    class Iterator;
+
+    /** Inheritable copy constructor */
+    inline NDGrid* clone();
+
+    /** Not implemented */
+    inline Mesh* dual() const;
+
+    /** Build a default, empty grid with zero dimension */
+    NDGrid();
+
+    /** Build a regular grid
+     *
+     *  @param nDimensions: number of dimensions of the grid
+     *  @param size: array of ints of length nDimensions. Gives the number of
+     *               grid points in each dimension. Caller owns the memory
+     *               pointed to by size. Must be at least 1 grid point in each
+     *               dimension.
+     *  @param spacing: array of doubles of length nDimensions. Gives the grid
+     *               spacing in each dimension. Caller owns memory pointed to by
+     *               spacing.
+     *  @param min: array of doubles of length nDimensions. Gives the coordinate
+     *              of the zeroth element (origin). Caller owns the memory
+     *              pointed to by min.
+     */
+    NDGrid(int nDimensions, int* size, double* spacing, double* min);
+
+    /** Build an irregular grid
+     *
+     *  @param size: The number of grid points in each dimension. This vector
+     *               must be the same length as gridCoordinates vector.
+     *  @param gridCoordinates: vector of doubles, with the same length as size.
+     *               Set the coordinates of the grid in each dimension. Caller
+     *               owns the memory pointed to by gridCoordinates.
+     */
+    NDGrid(std::vector<int> size, std::vector<const double *> gridCoordinates);
+
+    /** Build an irregular grid
+     *
+     *  @param gridCoordinates: Vector of size equal to the dimension of the
+     *                          grid. Element gridCoordinates[i][j] is the
+     *                          position of the jth grid coordinate in the ith
+     *                          dimension. i.e. to find the coordinates of the
+     *                          grid point (23, 17, 21) we would lookup
+     *                             gridCoordinate[0][23]
+     *                             gridCoordinate[1][17]
+     *                             gridCoordinate[2][21]
+     */
+    NDGrid(std::vector< std::vector<double> > gridCoordinates);
+
+    /** Destructor */
+    ~NDGrid() {;}
+
+    /** Get the coordinate of a set of grid points
+     *
+     *  @param index: the index of the grid point to be found; indexed from 1.
+     *         Should be in range 1 <= index <= length
+     *  @param dimension: the dimension of the grid point to be found. Should be
+     *         in range 0 <= dimension < grid dimension
+     *
+     *  Note the coord is NOT bound checked - so it is a memory error to go
+     *  out of the bounds.
+     */
+    inline double& coord(const int& index, const int& dimension);
+
+    /** Get the coordinate of a set of grid points
+     *
+     *  @param index: the index of the grid point to be found; indexed from 1.
+     *         Should be in range 1 <= index <= length
+     *  @param dimension: the dimension of the grid point to be found. Should be
+     *         in range 0 <= dimension < grid dimension
+     *
+     *  Note the call is NOT bound checked - so it is a memory error to go
+     *  out of the bounds.
+     */
+    inline const double& coord(const int& index, const int& dimension) const;
+
+    /** Get the size of the grid along a given dimension
+     *
+     *  @param dimension: the dimension of the grid point to be found. Should be
+     *         in range 0 <= dimension < grid dimension
+     *
+     *  Note the call is NOT bound checked - so it is a memory error to go
+     *  out of the bounds.
+     */
+    inline int size(const int& dimension) const;
+
+    /** Get the vector of grid points along a particular dimension
+     *
+     *  @param dimension: the dimension of the grid point to be found. Should be
+     *         in range 0 <= dimension < grid dimension
+     *
+     *  Note the call is NOT bound checked - so it is a memory error to go
+     *  out of the bounds.
+     */
+    inline std::vector<double> coordVector(const int& dimension) const;
+
+    /** Get an array of grid points along a particular dimension
+     *
+     *  @param dimension: the dimension of the grid point to be found. Should be
+     *         in range 0 <= dimension < grid dimension
+     *
+     *  Note the call is NOT bound checked - so it is a memory error to go
+     *  out of the bounds.
+     *  @returns an array of length this->size(dimension); caller owns the 
+     *  returned memory.
+     */
+    double* newCoordArray  (const int& dimension) const;
+
+    /** Get the grid point of the grid lower bound in a particular dimension
+     *
+     *  @param x: position in the grid coordinate
+     *  @param dimension: the grid coordinate to search along
+     *  @param xIndex: will be set with the resultant grid index
+     *
+     *  Sets xIndex to the index of the largest mesh position that is less than
+     *  x in the given dimension.
+     *
+     */
+    inline void coordLowerBound(const double& x, const int& dimension, int& xIndex) const;
+
+    /** Get the grid point of the grid lower bound
+     *
+     *  @param pos: position in the grid coordinate. Should have a length equal
+     *         to the dimension of the grid.
+     *  @param xIndex: will be set with the resultant grid index. Should have a
+     *         lenght equal to the dimension of the grid.
+     *
+     *  Sets xIndex to the index of the largest mesh position that is less than
+     *  pos in all dimension.
+     *
+     *  Note that pos and xIndex are not bound checked.
+     */
+    inline void lowerBound(const std::vector<double>& pos, std::vector<int>& xIndex) const;
+
+    /** Return the lowest grid coordinate in the given dimension */
+    inline double  min(const int& dimension) const;
+
+    /** Return the highest grid coordinate in the given dimension */
+    inline double  max(const int& dimension) const;
+
+    /** Reset the coordinates in the dimension to a new value
+     *
+     *  @param dimension: the dimension of the coordates that will be reset
+     *  @param nCoords: the length of array x
+     *  @param x: array of points. Caller owns the memory allocated to x.
+     */
+    inline void setCoord(int dimension, int nCoords, double * x);
+
+    /** Returns the origin of the mesh (lowest possible index in all dimensions)
+     */
+    inline Mesh::Iterator begin() const;
+
+    /** Returns the end of the mesh (highest possible index in all dimensions)
+     */
+    inline Mesh::Iterator end() const;
+
+    /** Get the position of the iterator
+     *  
+     *  @param it: iterator on this Mesh
+     *  @param position: double array of length equal to the mesh dimension.
+     *             Filled with the iterator position. Caller owns the memory
+     *             to which position points.
+     */
+    inline void getPosition(const Mesh::Iterator& it, double * position) const;
+
+    /** Get the dimension of the Mesh */
+    inline int getPositionDimension() const;
+
+    /** Returns true if the Mesh is a regular grid */
+    inline bool getConstantSpacing() const;
+
+    /** Set to true to flag the Mesh as a regular grid
+     *
+     *  @param spacing: flag indicating if the Mesh is regular.
+     *
+     *  Nb: If this is set for an irregular grid... I don't know what will
+     *  happen
+     */
+    inline void setConstantSpacing(bool spacing);
+
+    /** Automatically detector whether the Mesh is a regular grid
+     *
+     *  Loops over all grid points and checks the step size is regular, with
+     *  tolerance given by tolernance_m
+     */
+    void setConstantSpacing();
+ 
+    /** Return an integer corresponding to the iterator position
+     *
+     *  @param lhs: iterator whose position is returned.
+     *  
+     *  @returns an integer between 0 and number of grid points; 0 is
+     *  this->begin(), and each call to iterator++ increments the integer by 1
+     */
+    int toInteger(const Mesh::Iterator& lhs) const;
+
+    /** Return an iterator corresponding to nearest mesh point to position
+     *
+     *  @param position: Array of length dimension, giving a position near the
+     *                   grid.
+     *  
+     *  @returns iterator that points to the grid location closest to position.
+     */
+    Mesh::Iterator getNearest(const double* position) const;
+
+protected:
+
+    //Change position
+    virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, int difference) const;
+    virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, int difference) const;
+    virtual Mesh::Iterator& addEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
+    virtual Mesh::Iterator& subEquals(Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
+    virtual Mesh::Iterator& addOne   (Mesh::Iterator& lhs) const;
+    virtual Mesh::Iterator& subOne   (Mesh::Iterator& lhs) const;
+    //Check relative position
+    virtual bool isGreater(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs) const;
+
+private:
+    std::vector< std::vector<double> > coord_m;
+    std::vector<VectorMap*>            maps_m;
+    bool                               constantSpacing_m;
+    const double tolerance_m = 1e-9;
+
+    friend Mesh::Iterator  operator++(Mesh::Iterator& lhs, int);
+    friend Mesh::Iterator  operator--(Mesh::Iterator& lhs, int);
+    friend Mesh::Iterator& operator++(Mesh::Iterator& lhs);
+    friend Mesh::Iterator& operator--(Mesh::Iterator& lhs);
+    friend Mesh::Iterator  operator- (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend Mesh::Iterator  operator+ (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend Mesh::Iterator& operator-=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
+    friend Mesh::Iterator& operator+=(Mesh::Iterator& lhs,       const Mesh::Iterator& rhs);
+
+    friend bool operator==(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend bool operator!=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend bool operator>=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend bool operator<=(const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend bool operator< (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+    friend bool operator> (const Mesh::Iterator& lhs, const Mesh::Iterator& rhs);
+};
+
+double& NDGrid::coord(const int& index, const int& dimension) {
+    return coord_m[dimension][index-1];
+}
+
+const double& NDGrid::coord(const int& index, const int& dimension) const {
+    return coord_m[dimension][index-1];
+}
+
+int NDGrid::size(const int& dimension) const {
+    return coord_m[dimension].size();
+}
+
+std::vector<double> NDGrid::coordVector(const int& dimension) const {
+    return coord_m[dimension];
+}
+
+void NDGrid::coordLowerBound(const double& x, const int& dimension, int& xIndex) const {
+    if(constantSpacing_m) {
+        double x0 = coord_m[dimension][0];
+        double x1 = coord_m[dimension][1];
+        xIndex = static_cast<int>(floor((x - x0)/(x1-x0))); 
+    } else {
+        const std::vector<double>& c_t(coord_m[dimension]);
+        xIndex = std::lower_bound(c_t.begin(), c_t.end(), x)-c_t.begin()-1;
+    }
+}
+
+void NDGrid::lowerBound(const std::vector<double>& pos, std::vector<int>& xIndex) const {
+    xIndex = std::vector<int>(pos.size());
+    for(unsigned int i=0; i<pos.size(); i++) {
+        coordLowerBound(pos[i], i, xIndex[i]);
+    }
+}
+
+double NDGrid::min(const int& dimension) const {
+    return coord_m[dimension][0];
+}
+
+double NDGrid::max(const int& dimension) const {
+    return coord_m[dimension][coord_m[dimension].size()-1];
+}
+
+NDGrid* NDGrid::clone() {
+    return new NDGrid(*this);
+}
+
+Mesh* NDGrid::dual() const {
+    return NULL;
+}
+
+void NDGrid::setCoord(int dimension, int nCoords, double * x) {
+    coord_m[dimension] = std::vector<double>(x, x+nCoords);
+}
+
+Mesh::Iterator NDGrid::begin() const {
+    return Mesh::Iterator(std::vector<int>(coord_m.size(), 1), this);
+}
+
+Mesh::Iterator NDGrid::end() const {
+    std::vector<int> end(coord_m.size(), 1);
+    end[0] = coord_m[0].size()+1;
+    return Mesh::Iterator(end, this);
+}
+
+void NDGrid::getPosition(const Mesh::Iterator& it, double * position) const {
+    for (unsigned int i=0; i<it.getState().size(); i++)
+        position[i] = coord_m[i][it[i]-1];
+}
+
+int NDGrid::getPositionDimension() const {
+    return coord_m.size();
+}
+
+bool NDGrid::getConstantSpacing() const {
+    return constantSpacing_m;
+}
+
+void NDGrid::setConstantSpacing(bool spacing) {
+    constantSpacing_m = spacing;
+}
+
+} // namespace interpolation
+
+#endif // _CLASSIC_FIELDS_NDGRID_HH_
diff --git a/src/Classic/Fields/Interpolation/PPSolveFactory.cpp b/src/Classic/Fields/Interpolation/PPSolveFactory.cpp
index e4ecccf10..0f7d9aa9c 100644
--- a/src/Classic/Fields/Interpolation/PPSolveFactory.cpp
+++ b/src/Classic/Fields/Interpolation/PPSolveFactory.cpp
@@ -33,12 +33,16 @@
 #include <gsl/gsl_sf_pow_int.h>
 #include <gsl/gsl_sf_gamma.h>
 
+#include "Utility/Inform.h" // ippl
+
 #include "Utilities/GeneralClassicException.h"
 
 #include "Fields/Interpolation/SolveFactory.h"
 #include "Fields/Interpolation/PolynomialPatch.h"
 #include "Fields/Interpolation/PPSolveFactory.h"
 
+extern Inform* gmsg;
+
 namespace interpolation {
 
 typedef PolynomialCoefficient Coeff;
@@ -308,7 +312,7 @@ PolynomialPatch* PPSolveFactory::solve() {
          --it) {
         double newPercentage = (total-it.toInteger())/double(total)*100.;
         if (newPercentage - oldPercentage > 10.) {
-            std::cout << "    Done " << newPercentage << " %" << std::endl;
+            *gmsg << "    Done " << newPercentage << " %" << endl;
             oldPercentage = newPercentage;
         }
         // find the set of values and derivatives for this point in the mesh
diff --git a/src/Classic/Fields/SectorMagneticFieldMap.cpp b/src/Classic/Fields/SectorMagneticFieldMap.cpp
index 2d70e8440..c3d72116f 100644
--- a/src/Classic/Fields/SectorMagneticFieldMap.cpp
+++ b/src/Classic/Fields/SectorMagneticFieldMap.cpp
@@ -140,7 +140,7 @@ void SectorMagneticFieldMap::setInterpolator(VectorMap* interpolator) {
                                      grid->maxX(), grid->maxY(), grid->maxZ(),
                                      0., 0., 0.);
     }
-    print(std::cerr);
+    getInfo(gmsg);
 }
 
 std::string SectorMagneticFieldMap::getSymmetry() const {
@@ -152,8 +152,10 @@ void SectorMagneticFieldMap::setSymmetry(std::string name) {
 }
 
 void SectorMagneticFieldMap::readMap() {
-    setInterpolator(IO::readMap
-            (filename_m, units_m, symmetry_m, poly_order_m, smoothing_order_m));
+    VectorMap* interpolator = IO::readMap
+             (filename_m, units_m, symmetry_m, poly_order_m, smoothing_order_m);
+    setInterpolator(interpolator);
+
 }
 
 void SectorMagneticFieldMap::freeMap() {
@@ -306,14 +308,16 @@ VectorMap* SectorMagneticFieldMap::IO::readMap(
         ThreeDGrid* grid = generateGrid(field_points, sym);
         // build interpolator (convert grid to useful units)
         if (polynomial_order == 1 && smoothing_order == 1) {
-            return getInterpolator(field_points, grid, sym);
+            VectorMap* interpolator = getInterpolator(field_points, grid, sym);
+            return interpolator;
         } else {
-            return getInterpolatorPolyPatch(
+            VectorMap* interpolator = getInterpolatorPolyPatch(
                                   field_points,
                                   grid,
                                   sym,
                                   polynomial_order,
                                   smoothing_order);
+            return interpolator;
         }
     } catch(std::exception& exc) {
         throw(LogicalError(
@@ -346,14 +350,15 @@ VectorMap* SectorMagneticFieldMap::IO::getInterpolatorPolyPatch(
     }
     // symmetry is dipole
     try {
-        INFOMSG("Calculating polynomials..." << endl);
+        *gmsg << "Calculating polynomials..." << endl;
         PPSolveFactory solver(grid, data, polynomial_order, smoothing_order);
         PolynomialPatch* patch = solver.solve();
-        INFOMSG("                       ... done" << endl);
+        *gmsg << "                       ... done" << endl;
         return patch;
     } catch (GeneralClassicException& exc) {
         throw exc;
     }
+
 }
 
 VectorMap* SectorMagneticFieldMap::IO::getInterpolator
@@ -395,7 +400,8 @@ VectorMap* SectorMagneticFieldMap::IO::getInterpolator
                ));
     }
     Interpolator3dGridTo3d* interpolator = new Interpolator3dGridTo3d(grid, Bx, By, Bz);
-    return dynamic_cast<VectorMap*>(interpolator);
+    VectorMap* interpolatorMap = dynamic_cast<VectorMap*>(interpolator);
+    return interpolatorMap;
 }
 
 bool SectorMagneticFieldMap::IO::comparator(std::vector<double> field_item1,
diff --git a/src/Classic/Utilities/RingSection.cpp b/src/Classic/Utilities/RingSection.cpp
index fec374195..81a78d1b2 100644
--- a/src/Classic/Utilities/RingSection.cpp
+++ b/src/Classic/Utilities/RingSection.cpp
@@ -95,14 +95,11 @@ bool RingSection::isPastEndPlane(const Vector_t& pos) const {
 bool RingSection::getFieldValue(const Vector_t& pos,
                                 const Vector_t& centroid, const double& t,
                                 Vector_t& E, Vector_t& B) const {
-    // std::cerr << startOrientation_m << " " << componentPosition_m << " " << componentOrientation_m(2) << std::endl;
     // transform position into local coordinate system
     Vector_t pos_local = pos-componentPosition_m;
     rotate(pos_local);
-    // std::cerr << "RingSection::getFieldValue for " << component_m->getName() << " at " << startPosition_m << "\n    Global: " << pos << " Local: " << pos_local << std::endl;
     rotateToTCoordinates(pos_local);
     bool outOfBounds = component_m->apply(pos_local, Vector_t(0.0), t, E, B);
-    // std::cerr << "    B: " << B << std::endl;
     if (outOfBounds) {
         return true;
     }
@@ -194,4 +191,4 @@ void RingSection::rotate_back(Vector_t& vector) const {
     const Vector_t temp(vector);
     vector(0) = +cos2_m * temp(0) - sin2_m * temp(1);
     vector(1) = +sin2_m * temp(0) + cos2_m * temp(1);
-}
\ No newline at end of file
+}
diff --git a/src/OpalConfigure/Configure.cpp b/src/OpalConfigure/Configure.cpp
index e27759d63..463465538 100644
--- a/src/OpalConfigure/Configure.cpp
+++ b/src/OpalConfigure/Configure.cpp
@@ -26,6 +26,7 @@
 #include "BasicActions/Call.h"
 #include "BasicActions/Dump.h"
 #include "BasicActions/DumpFields.h"
+#include "BasicActions/DumpEMFields.h"
 #include "BasicActions/Echo.h"
 #include "BasicActions/Help.h"
 #include "BasicActions/Option.h"
@@ -271,4 +272,4 @@ namespace Configure {
         makeActions();
         Versions::fillChanges();
     }
-};
\ No newline at end of file
+};
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4841084b9..c44a2510b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -26,6 +26,8 @@ endmacro()
 
 add_subdirectory (opal_src)
 add_subdirectory (classic_src)
+add_subdirectory (opal_test_utilities)
+
 
 set(TEST_SRCS ${CMAKE_CURRENT_SOURCE_DIR}/Main.cpp ${TEST_SRCS_LOCAL})
 
@@ -60,4 +62,4 @@ INCLUDE_DIRECTORIES(${GTEST_INCLUDE_DIR}
 # Build the test exe. We don't do an install on the unit test exe as it is
 # assumed that this is internal to opal
 ADD_EXECUTABLE(${TEST_EXE} ${TEST_SRCS}) # the opal and classic sources are not needed again if we link agains libOPAL and libCLASSIC!
-TARGET_LINK_LIBRARIES(${TEST_EXE} OPALib ${OPAL_LIBS} ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES} ${CCSE_LIBRARIES} -lgfortran ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${GTEST_BOTH_LIBRARIES} -lpthread)
\ No newline at end of file
+TARGET_LINK_LIBRARIES(${TEST_EXE} OPALib ${OPAL_LIBS} ${Trilinos_LIBRARIES} ${Trilinos_TPL_LIBRARIES} ${CCSE_LIBRARIES} -lgfortran ${OTHER_CMAKE_EXE_LINKER_FLAGS} ${GTEST_BOTH_LIBRARIES} -lpthread)
diff --git a/tests/Main.cpp b/tests/Main.cpp
index b265ef6ac..bab09f5b8 100644
--- a/tests/Main.cpp
+++ b/tests/Main.cpp
@@ -2,11 +2,19 @@
 
 #include "mpi.h"
 
+#include "Utility/Inform.h" // ippl
+
+extern Inform* gmsg;
+
 int main(int argc, char **argv) {
     ::testing::InitGoogleTest(&argc, argv);
+    gmsg = new Inform("UnitTests: ", std::cerr);
+    if (!gmsg) {
+        return 1;
+    }
     MPI_Init(&argc, &argv);
     int test_out = RUN_ALL_TESTS();
     MPI_Finalize();
 
     return test_out;
-}
\ No newline at end of file
+}
diff --git a/tests/classic_src/AbsBeamline/RingTest.cpp b/tests/classic_src/AbsBeamline/RingTest.cpp
index 7588f0d3f..7543e89ec 100644
--- a/tests/classic_src/AbsBeamline/RingTest.cpp
+++ b/tests/classic_src/AbsBeamline/RingTest.cpp
@@ -27,6 +27,8 @@
 
 #include "gtest/gtest.h"
 
+#include "Algorithms/PartData.h"
+#include "Algorithms/PartBunch.h"
 #include "AbsBeamline/Offset.h"
 #include "opal_src/Utilities/MockComponent.h"
 #include "AbsBeamline/Ring.h"
@@ -110,21 +112,12 @@ class OffsetFactory {
     std::vector<MockComponent*> mockVec_m;
 };
 
-
 TEST(RingTest, TestConstructDestruct) {
     // something here? someday...
 }
 
 TEST(RingTest, TestAppend1) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 5.;
         Ring ring("my_ring");
         ring.setLatticeRInit(radius);
@@ -135,43 +128,29 @@ TEST(RingTest, TestAppend1) {
         Offset off = Offset::localCylindricalOffset("cyl1", 0., Physics::pi/6., 1.);
         ring.appendElement(off);
         for (int i = 0; i < 3; ++i) {
-            EXPECT_NEAR(ring.getNextPosition()(i), Vector_t(5., -1., 0.)(i), 1e-6)  << ::burnAfterReading(debugOutput);
+            EXPECT_NEAR(ring.getNextPosition()(i), Vector_t(5., -1., 0.)(i), 1e-6);
             EXPECT_NEAR(ring.getNextNormal()(i), Vector_t(-sin(Physics::pi/6.),
                                                           -cos(Physics::pi/6.),
-                                                          0.)(i), 1e-6) << ::burnAfterReading(debugOutput);
+                                                          0.)(i), 1e-6);
         }
         ring.appendElement(off);
         for (int i = 0; i < 3; ++i) {
             EXPECT_NEAR(ring.getNextPosition()(i),
                         Vector_t(5.-sin(Physics::pi/6.),
-                                 -1.-cos(Physics::pi/6.), 0.)(i), 1e-6)  << ::burnAfterReading(debugOutput);
-            EXPECT_NEAR(ring.getNextNormal()(i), Vector_t(-sin(Physics::pi/3.),
-                                                         -cos(Physics::pi/3.),
-                                                          0.)(i), 1e-6)  << ::burnAfterReading(debugOutput);
+                                 -1.-cos(Physics::pi/6.), 0.)(i), 1e-6);
+            Vector_t expected(-sin(Physics::pi/3.), -cos(Physics::pi/3.), 0.);
+            EXPECT_NEAR(ring.getNextNormal()(i),
+                        expected(i),
+                        1e-6);
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cerr << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n" << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n";
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 TEST(RingTest, TestAppend2) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 5.;
         Ring ring("my_ring");
         ring.setLatticeRInit(radius);
@@ -185,44 +164,26 @@ TEST(RingTest, TestAppend2) {
             EXPECT_NEAR(ring.getNextPosition()(i),
                         Vector_t(cos(Physics::pi/24.),
                                  5.-sin(Physics::pi/24.), 0.)(i), 1e-6)
-                << i << "\n"
-                << ::burnAfterReading(debugOutput);
+                << i << "\n";
             EXPECT_NEAR(ring.getNextNormal()(i), Vector_t(cos(Physics::pi/6.),
                                                          -sin(Physics::pi/6.),
                                                           0.)(i), 1e-6)
-                << i << "\n"
-                << ::burnAfterReading(debugOutput);
+                << i << "\n";
         }
         ring.appendElement(off);
         ring.appendElement(off);
         for (int i = 0; i < 3; ++i) {
             EXPECT_NEAR(ring.getNextNormal()(i), Vector_t(0., -1., 0.)(i), 1e-6)
-                << i << "\n"
-                << ::burnAfterReading(debugOutput);
+                << i << "\n";
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cerr << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n" << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n";
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 TEST(RingTest, TestAppend3) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 5.;
         Ring ring("my_ring");
         ring.setLatticeRInit(radius);
@@ -235,26 +196,13 @@ TEST(RingTest, TestAppend3) {
             ring.appendElement(off);
         }
         ring.lockRing();
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cerr << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n" << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n";
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 TEST(RingTest, TestLatticeRInitPhiInit) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     for (double phi = -2.*Physics::pi;
          phi < 2.*Physics::pi;
          phi += Physics::pi/6.) {
@@ -271,38 +219,28 @@ TEST(RingTest, TestLatticeRInitPhiInit) {
                 for (size_t i = 0; i < 3; ++i) {
                     EXPECT_EQ(pos(i), refPos(i))
                         << i << " f: " << phi
-                        << " t: " << theta << " r: " << radius << "\n"
-                        << ::burnAfterReading(debugOutput);
+                        << " t: " << theta << " r: " << radius << "\n";
                 }
                 Vector_t norm = ring.getNextNormal();
                 Vector_t refNorm(cos(phi+theta), -sin(phi+theta), 0.);
                 for (size_t i = 0; i < 3; ++i) {
                     EXPECT_EQ(norm(i), refNorm(i))
                         << i << " f: " << phi
-                        << " t: " << theta << " r: " << radius << "\n"
-                        << ::burnAfterReading(debugOutput);
+                        << " t: " << theta << " r: " << radius << "\n";
                 }
             }
         }
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 // Check that we get the bounding box and rotation correct
 TEST(RingTest, TestApply) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     Ring ring("my_ring");
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 2.*(2.*sin(Physics::pi/6.)+1.*sin(Physics::pi/3.)+1.0);
+        PartData data;
+        PartBunch bunch(&data);
+        ring.setRefPartBunch(&bunch);
         ring.setLatticeRInit(radius-2.);
         ring.setLatticePhiInit(0.);
         ring.setLatticeThetaInit(0.);
@@ -317,49 +255,41 @@ TEST(RingTest, TestApply) {
         for (double x = -1.0001; x < 2.; x += 0.1) {
             double y = -2.2;
             Vector_t pos(x, y, -0.5);
-            Vector_t centroid, B, E;
+            Vector_t zero(0.0);
+            Vector_t centroid(0., 0., 0), B(0., 0., 0), E(0., 0., 0);
+            double t = 0;
             // std::cout << pos << " ** " << std::flush;
-            EXPECT_FALSE(ring.apply(pos, Vector_t(0.0), 0., E, B)) << " for pos " << pos << ::burnAfterReading(debugOutput);
+            EXPECT_FALSE(ring.apply(pos, zero, t, E, B));
             // std::cout << B << " " << E << std::endl;
             Vector_t BRef(0.0, 0.0, 0.0);
             if (x > 0. and x < 1.)
                 BRef = Vector_t(x-1., y+2., -0.5);
             for (int i = 0; i < 3; ++i) {
-                EXPECT_NEAR(B(i), BRef(i), 1e-6)   << " for pos " << pos << ::burnAfterReading(debugOutput);
-                EXPECT_NEAR(E(i), -BRef(i), 1e-6)  << " for pos " << pos << ::burnAfterReading(debugOutput);
+                EXPECT_NEAR(B(i), BRef(i), 1e-6)   << " for pos " << pos;
+                EXPECT_NEAR(E(i), -BRef(i), 1e-6)  << " for pos " << pos;
             }
         }
         // check that we get something reasonable for all phi
         for (double phi = 0.; phi < 2.*Physics::pi+0.1; phi += Physics::pi/100.) {
             Vector_t pos(radius/2.*sin(phi), radius/2.+radius/2.*cos(phi), 0.5);
             Vector_t centroid, B, E;
-            EXPECT_FALSE(ring.apply(pos, Vector_t(0.0), 0., E, B)) << ::burnAfterReading(debugOutput); // check we don't throw for all angles
+            EXPECT_FALSE(ring.apply(pos, Vector_t(0.0), 0., E, B)); // check we don't throw for all angles
             // std::cout << phi << " " << pos << " " << B << std::endl;
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cout << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n" << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n";
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 // Check that we get the bounding box correct - for exact sector geometry
 TEST(RingTest, TestApply2) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     Ring ring("my_ring");
     try {
         double radius = 1.5;
+        PartData data;
+        PartBunch bunch(&data);
+        ring.setRefPartBunch(&bunch);
         ring.setLatticeRInit(radius);
         ring.setLatticePhiInit(7.*Physics::pi/4.);
         ring.setLatticeThetaInit(0.);
@@ -374,29 +304,26 @@ TEST(RingTest, TestApply2) {
             Vector_t pos((radius+0.5)*sin(phi), (radius+0.5)*cos(phi), -0.5);
             Vector_t centroid, B, E;
             std::vector<RingSection*> sections = ring.getSectionsAt(pos);
-            EXPECT_FALSE(ring.apply(pos, Vector_t(0.0), 0., E, B)) << ::burnAfterReading(debugOutput);
+            EXPECT_FALSE(ring.apply(pos, Vector_t(0.0), 0., E, B));
             // check we don't throw for all angles
             // a few are coming out with Bz = 1. instead of Bz = 0.5; looks like
             // floating point precision issue? It's okay, Ring is not
             // responsible for bounding the field, Components are.
-            EXPECT_GE(-B(2), 0.1) << ::burnAfterReading(debugOutput);
-            EXPECT_LE(-B(2), 1.1) << ::burnAfterReading(debugOutput);
+            EXPECT_GE(-B(2), 0.1);
+            EXPECT_LE(-B(2), 1.1);
             // std::cout << phi << " " << pos << " " << B << " " << sections.size() << std::endl;
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cout << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n"  << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n" ;
     }
     // Now apply symmetry 2x10 fields instead of 20x1
     Ring ring2("my_ring");
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 1.5;
+        PartData data;
+        PartBunch bunch(&data);
+        ring2.setRefPartBunch(&bunch);
         ring2.setLatticeRInit(radius);
         ring2.setLatticePhiInit(7.*Physics::pi/4.);
         ring2.setLatticeThetaInit(0.);
@@ -413,24 +340,21 @@ TEST(RingTest, TestApply2) {
             std::vector<RingSection*> sections = ring2.getSectionsAt(pos);
             ring.apply(pos, Vector_t(0.0), 0., E, B1);
             ring2.apply(pos, Vector_t(0.0), 0., E, B2);
-            EXPECT_NEAR(B1(2), B2(2), 1e-6) << ::burnAfterReading(debugOutput);
+            EXPECT_NEAR(B1(2), B2(2), 1e-6);
             // std::cout << phi << " " << pos << " " << B << " " << sections.size() << std::endl;
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cout << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n"  << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n" ;
     }
     // Now overlapping - we have two elements in each position, should get twice
     // the field
     Ring ring3("my_ring");
     try {
-        std::cout.setstate(std::ios::failbit);
-        std::cerr.setstate(std::ios::failbit);
         double radius = 1.5;
+        PartData data;
+        PartBunch bunch(&data);
+        ring3.setRefPartBunch(&bunch);
         ring3.setLatticeRInit(radius);
         ring3.setLatticePhiInit(7.*Physics::pi/4.);
         ring3.setLatticeThetaInit(0.);
@@ -446,32 +370,22 @@ TEST(RingTest, TestApply2) {
             std::vector<RingSection*> sections = ring3.getSectionsAt(pos);
             ring.apply(pos, Vector_t(0.0), 0., E, B1);
             ring3.apply(pos, Vector_t(0.0), 0., E, B2);
-            EXPECT_NEAR(2.*B1(2), B2(2), 1e-6) << ::burnAfterReading(debugOutput);
+            EXPECT_NEAR(2.*B1(2), B2(2), 1e-6);
             // std::cout << phi << " " << pos << " " << B << " " << sections.size() << std::endl;
         }
-        std::cout.clear();
-        std::cerr.clear();
     } catch (OpalException& exc) {
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         std::cout << exc.what() << std::endl;
-        EXPECT_TRUE(false) << "Threw an exception\n" << ::burnAfterReading(debugOutput);
+        EXPECT_TRUE(false) << "Threw an exception\n";
     }
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 void testField(double s, double r, double y, double phi,
                double bx, double by, double bz, double tol) {
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
     double radius = 2.;
     Ring ring("test");
+    PartData data;
+    PartBunch bunch(&data);
+    ring.setRefPartBunch(&bunch);
     ring.setLatticeRInit(radius);
     ring.setLatticePhiInit(phi);
     ring.setLatticeThetaInit(0.);
@@ -487,18 +401,13 @@ void testField(double s, double r, double y, double phi,
                  radius*cos(phi)-s*sin(phi)+r*cos(phi),
                  y);
     ring.apply(pos, Vector_t(0.0), 0., E, B);
-    EXPECT_NEAR(B(0), bx, 1e-6) << ::burnAfterReading(debugOutput);
-    EXPECT_NEAR(B(1), by, 1e-6) << ::burnAfterReading(debugOutput);
-    EXPECT_NEAR(B(2), bz, 1e-6) << ::burnAfterReading(debugOutput);
+    EXPECT_NEAR(B(0), bx, 1e-6);
+    EXPECT_NEAR(B(1), by, 1e-6);
+    EXPECT_NEAR(B(2), bz, 1e-6);
     // std::cout << pos << " ** " << B << " ** " << Vector_t(bx, by, bz) << std::endl;
-
-    std::cout.rdbuf(defaultCout);
-    std::cerr.rdbuf(defaultCerr);
 }
 
 TEST(RingTest, TestApply3) {
-    std::cout.setstate(std::ios::failbit);
-    std::cerr.setstate(std::ios::failbit);
     testField(0.1, 0., 0.2, 0., 3., 1., 2., 1e-6);
     testField(0.1, 0., 0.2, Physics::pi, -3., -1., 2., 1e-6);
     testField(0.1, 0., 0.2, Physics::pi/2., 1., -3., 2., 1e-6);
@@ -506,7 +415,4 @@ TEST(RingTest, TestApply3) {
     testField(0.1, 0.15, 0.2, Physics::pi/6.,
               3.*cos(Physics::pi/6)+1.*sin(Physics::pi/6),
               -3.*sin(Physics::pi/6)+1.*cos(Physics::pi/6), 2., 1e-6);
-
-    std::cout.clear();
-    std::cerr.clear();
-}
\ No newline at end of file
+}
diff --git a/tests/classic_src/AbsBeamline/SBend3DTest.cpp b/tests/classic_src/AbsBeamline/SBend3DTest.cpp
index 59489a009..f41ffc0f9 100644
--- a/tests/classic_src/AbsBeamline/SBend3DTest.cpp
+++ b/tests/classic_src/AbsBeamline/SBend3DTest.cpp
@@ -245,13 +245,9 @@ TEST(SBend3DTest, SBend3DPolyPatchTest) {
                 for (double phi =-dphi/2.; phi < dphi*3.1; phi += dphi/2.) {
                     Vector_t B, E, centroid;
                     Vector_t pos(r*cos(phi)-radius, y, r*sin(phi));
-                    std::cerr << "pos: " << pos << " r: " << r << " phi: " << phi/Physics::pi*180. << std::endl;
                     field1->apply(pos, Vector_t(0.0), 0, E, B);
-                    std::cerr << "   field1: " << B << std::endl;
                     field2->apply(pos, Vector_t(0.0), 0, E, B);
-                    std::cerr << "   field2: " << B << std::endl;
                     field3->apply(pos, Vector_t(0.0), 0, E, B);
-                    std::cerr << "   field3: " << B << std::endl;
                 }
     } catch (LogicalError& err) {
         std::cerr << err.what() << std::endl;
diff --git a/tests/classic_src/Fields/Interpolation/NDGridTest.cpp b/tests/classic_src/Fields/Interpolation/NDGridTest.cpp
new file mode 100644
index 000000000..c7f6a6d6d
--- /dev/null
+++ b/tests/classic_src/Fields/Interpolation/NDGridTest.cpp
@@ -0,0 +1,36 @@
+/* 
+ *  Copyright (c) 2017, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without 
+ *  modification, are permitted provided that the following conditions are met: 
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer. 
+ *  2. Redistributions in binary form must reproduce the above copyright notice, 
+ *     this list of conditions and the following disclaimer in the documentation 
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to 
+ *     endorse or promote products derived from this software without specific 
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include "gtest/gtest.h"
+#include "Fields/Interpolation/NDGrid.h"
+
+using interpolation::ThreeDGrid;
+
+namespace ndgridtest {
+
+}
+
diff --git a/tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp b/tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp
index f7af14d28..643628863 100644
--- a/tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp
+++ b/tests/classic_src/Fields/Interpolation/ThreeDGridTest.cpp
@@ -30,6 +30,8 @@
 
 using interpolation::ThreeDGrid;
 
+namespace threegridtest {
+
 TEST(ThreeDGridTest, LowerBoundTest) {
     std::vector<double> xVar(7), yVar(8), zVar(9);
     for (size_t i = 0; i < xVar.size(); ++i) {
@@ -82,3 +84,4 @@ TEST(ThreeDGridTest, LowerBoundTest) {
         }
     }
 }
+}
diff --git a/tests/opal_src/BasicActions/CMakeLists.txt b/tests/opal_src/BasicActions/CMakeLists.txt
index 25628c4b5..17f19d8b6 100644
--- a/tests/opal_src/BasicActions/CMakeLists.txt
+++ b/tests/opal_src/BasicActions/CMakeLists.txt
@@ -1,5 +1,6 @@
 set (_SRCS
     DumpFieldsTest.cpp
+    DumpEMFieldsTest.cpp
 )
 
 include_directories (
diff --git a/tests/opal_src/BasicActions/DumpEMFieldsTest.cpp b/tests/opal_src/BasicActions/DumpEMFieldsTest.cpp
new file mode 100644
index 000000000..90065f124
--- /dev/null
+++ b/tests/opal_src/BasicActions/DumpEMFieldsTest.cpp
@@ -0,0 +1,180 @@
+/*
+ *  Copyright (c) 2014, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are met:
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *  2. Redistributions in binary form must reproduce the above copyright notice,
+ *     this list of conditions and the following disclaimer in the documentation
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to
+ *     endorse or promote products derived from this software without specific
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <fstream>
+#include <iostream>
+
+#include "gtest/gtest.h"
+
+#include "opal_src/Utilities/MockComponent.h"
+
+#include "Attributes/Attributes.h"
+#include "Utilities/OpalException.h"
+#include "BasicActions/DumpEMFields.h"
+
+namespace DumpEMFieldsTest {
+
+void setOneAttribute(DumpEMFields* dump, std::string name, double value) {
+    Attributes::setReal(*dump->findAttribute(name), value);
+}
+
+void setAttributes(DumpEMFields* dump,
+                   double x0, double dx, double nx,
+                   double y0, double dy, double ny,
+                   double z0, double dz, double nz,
+                   double t0, double dt, double nt,
+                   std::string filename) {
+    setOneAttribute(dump, "X_START", x0);
+    setOneAttribute(dump, "DX", dx);
+    setOneAttribute(dump, "X_STEPS", nx);
+    setOneAttribute(dump, "Y_START", y0);
+    setOneAttribute(dump, "DY", dy);
+    setOneAttribute(dump, "Y_STEPS", ny);
+    setOneAttribute(dump, "Z_START", z0);
+    setOneAttribute(dump, "DZ", dz);
+    setOneAttribute(dump, "Z_STEPS", nz);
+    setOneAttribute(dump, "T_START", t0);
+    setOneAttribute(dump, "DT", dt);
+    setOneAttribute(dump, "T_STEPS", nt);
+    Attributes::setString(*dump->findAttribute("FILE_NAME"), filename);
+}
+
+TEST(DumpEMFieldsTest, ConstructorDestructor) {
+    // neither in the set and grid is null
+    DumpEMFields* dump1 = new DumpEMFields();
+    delete dump1;
+    // grid is not null and it is in the set
+    DumpEMFields* dump2 = new DumpEMFields();
+    setAttributes(dump2, 1., 1., 1.,   1., 1., 1.,   1., 1., 1.,   1., 1., 1.,  "/dev/null");
+    dump2->execute();
+    delete dump2;
+}
+
+void execute_throws(DumpEMFields* dump, std::string reason) {
+    try {
+        dump->execute();
+        EXPECT_TRUE(false) << reason;
+    } catch (OpalException& exc) {
+        // pass;
+    }
+}
+
+TEST(DumpEMFieldsTest, executeTest) {
+    // dump the fields
+    DumpEMFields dump1;
+    execute_throws(&dump1, "should throw due to nsteps < 1");
+    setAttributes(&dump1, 1., 1., 1.,   1., 1., 1.,   1., 1., 1.,   1., 1., 1.,  "/dev/null");
+    dump1.execute();  // should be okay (normal)
+    setAttributes(&dump1, -1., -1., 1.,   -1., -1., 1.,   -1., -1., 1.,   1., 1., 1.,  "/dev/null");
+    dump1.execute();  // should be okay (-ve step is okay)
+    setAttributes(&dump1, -1., -1., 0.,   -1., -1., 1.,   -1., -1., 1.,   1., 1., 1.,  "/dev/null");
+    execute_throws(&dump1, "should throw due to nsteps x < 1");
+    setAttributes(&dump1, -1., -1., 1.,   -1., -1., 0.,   -1., -1., 1.,   1., 1., 1.,  "/dev/null");
+    execute_throws(&dump1, "should throw due to nsteps y < 1");
+    setAttributes(&dump1, -1., -1., 1.,   -1., -1., 1.,   -1., -1., 0.,   1., 1., 1.,  "/dev/null");
+    execute_throws(&dump1, "should throw due to nsteps z < 1");
+    setAttributes(&dump1, -1., -1., 1.,   -1., -1., 1.,   -1., -1., 1.,   1., 1., 0.,  "/dev/null");
+    execute_throws(&dump1, "should throw due to nsteps t < 1");
+    setAttributes(&dump1, -1., -1., 1.,   -1., -1., 1.,   -1., -1., 1.5,   1., 1., 1.,  "/dev/null");
+    execute_throws(&dump1, "should throw due to nsteps not integer");
+}
+
+void clear_files() {
+    size_t n_str_array = 4;
+    std::string str_array[4] = {"test1", "test2", "test3", "test4"};
+    for (size_t i = 0; i < n_str_array; ++i) {
+        if (fopen(str_array[i].c_str(), "r") != NULL) {
+            remove(str_array[i].c_str());
+        }
+    }
+}
+
+TEST(DumpEMFieldsTest, writeFieldsTest) {
+    clear_files();
+    DumpEMFields dump1;
+    setAttributes(&dump1, 1., 1., 1.,   1., 1., 1.,   1., 1., 1.,   1., 1., 1.,  "test1");
+    dump1.execute();
+    DumpEMFields dump2;
+    setAttributes(&dump2, 1., 1., 1.,   1., 1., 1.,   1., 1., 1.,   1., 1., 1.,  "test2");
+    dump2.execute();
+    DumpEMFields dump3;
+    setAttributes(&dump3, 1., 1., 1.,   1., 1., 1.,   1., 1., 1.,   1., 1., 1.,  "test3");
+    // note we don't execute dump3; so it should not be written
+    DumpEMFields dump4;
+    setAttributes(&dump4, 0.1, 0.1, 3.,   -0.1, 0.2, 2.,   0.2, 0.3, 2.,   1., 1., 2.,  "test4");
+    dump4.execute();
+    MockComponent comp;
+    try {
+        DumpEMFields::writeFields(&comp);
+    } catch (OpalException& exc) {
+        EXPECT_TRUE(false) << "Threw OpalException on writefields: " << exc.what() << std::endl;;
+    }
+    std::ifstream fin1("test1");
+    EXPECT_TRUE(fin1.good());
+    std::ifstream fin2("test2");
+    EXPECT_TRUE(fin2.good());
+    std::ifstream fin3("test3");
+    EXPECT_FALSE(fin3.good());  // does not exist
+    std::ifstream fin4("test4");
+    EXPECT_TRUE(fin4.good());
+    int n_lines;
+    fin4 >> n_lines;
+    EXPECT_EQ(n_lines, 24);
+    std::string test_line;
+    for (size_t i = 0; i < 12; ++i) {
+        std::getline(fin4, test_line);
+    }
+    std::vector<double> line(10, 0.);
+    double tol = 1e-9;
+    for (size_t line_index = 0; line_index < 24; ++line_index) {
+        for (size_t i = 0; i < 10; ++i) {
+            fin4 >> line[i];
+        }
+        if (line_index == 0) {
+            EXPECT_NEAR(line[0], 0.1, tol);
+            EXPECT_NEAR(line[1], -0.1, tol);
+            EXPECT_NEAR(line[2], 0.2, tol);
+            EXPECT_NEAR(line[3], 1., tol);
+        }
+        if (line[1] < 0.) {
+            EXPECT_NEAR(line[4], line[0], tol);
+            EXPECT_NEAR(line[5], line[1], tol);
+            EXPECT_NEAR(line[6], line[2], tol);
+        } else {
+            EXPECT_NEAR(line[4], 0., tol);
+            EXPECT_NEAR(line[5], 0., tol);
+            EXPECT_NEAR(line[6], 0., tol);
+        }
+        EXPECT_NEAR(line[7], -line[4], tol);
+        EXPECT_NEAR(line[8], -line[5], tol);
+        EXPECT_NEAR(line[9], -line[6], tol);
+    }
+    clear_files();
+}
+
+}
+
diff --git a/tests/opal_src/Distribution/GaussTest.cpp b/tests/opal_src/Distribution/GaussTest.cpp
index 9749aef33..ae89e60b7 100644
--- a/tests/opal_src/Distribution/GaussTest.cpp
+++ b/tests/opal_src/Distribution/GaussTest.cpp
@@ -14,6 +14,8 @@ Inform *gmsg;
 #include "Parser/FileStream.h"
 #include "Physics/Physics.h"
 
+#include "opal_test_utilities/SilenceTest.h"
+
 #include "gsl/gsl_statistics_double.h"
 
 #include <cstdio>
@@ -33,13 +35,7 @@ namespace {
 }
 
 TEST(GaussTest, FullSigmaTest1) {
-
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
+    OpalTestUtilities::SilenceTest silencer(true);
     char inputFileName[] = "GaussDistributionTest.in";
     std::string input = "OPTION, ECHO=FALSE;\n"
         "OPTION, CZERO=FALSE;\n"
@@ -85,8 +81,6 @@ TEST(GaussTest, FullSigmaTest1) {
         is = new FileStream(inputFileName);
     } catch(...) {
         is = 0;
-        std::cout.rdbuf(defaultCout);
-        std::cerr.rdbuf(defaultCerr);
         throw new OpalException("FullSigmaTest", "Could not read string");
     }
 
@@ -95,8 +89,6 @@ TEST(GaussTest, FullSigmaTest1) {
         try {
             parser->run(is);
         } catch (...) {
-            std::cout.rdbuf(defaultCout);;
-            std::cerr.rdbuf(defaultCerr);;
             throw new OpalException("FullSigmaTest", "Could not parse input");
         }
     }
@@ -106,8 +98,6 @@ TEST(GaussTest, FullSigmaTest1) {
         distObj = OPAL->find("DIST1");
     } catch(...) {
         distObj = 0;
-        std::cout.rdbuf(defaultCout);;
-        std::cerr.rdbuf(defaultCerr);;
         throw new OpalException("FullSigmaTest", "Could not find distribution");
     }
 
@@ -138,17 +128,15 @@ TEST(GaussTest, FullSigmaTest1) {
         const double expectedR61 = 1.362;
         const double expectedR62 = -0.2685;
 
-        EXPECT_LT(std::abs(expectedR11 - R11),  0.05 * expectedR11) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR21 - R21), -0.05 * expectedR21) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR22 - R22),  0.05 * expectedR22) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR51 - R51),  0.05 * expectedR51) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR52 - R52),  0.05 * expectedR52) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR61 - R61),  0.05 * expectedR61) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR62 - R62), -0.05 * expectedR62) << ::burnAfterReading(debugOutput);
+        EXPECT_LT(std::abs(expectedR11 - R11),  0.05 * expectedR11);
+        EXPECT_LT(std::abs(expectedR21 - R21), -0.05 * expectedR21);
+        EXPECT_LT(std::abs(expectedR22 - R22),  0.05 * expectedR22);
+        EXPECT_LT(std::abs(expectedR51 - R51),  0.05 * expectedR51);
+        EXPECT_LT(std::abs(expectedR52 - R52),  0.05 * expectedR52);
+        EXPECT_LT(std::abs(expectedR61 - R61),  0.05 * expectedR61);
+        EXPECT_LT(std::abs(expectedR62 - R62), -0.05 * expectedR62);
     }
 
-    std::cout.rdbuf(defaultCout);;
-    std::cerr.rdbuf(defaultCerr);;
 
     OpalData::deleteInstance();
     delete parser;
@@ -160,13 +148,7 @@ TEST(GaussTest, FullSigmaTest1) {
 }
 
 TEST(GaussTest, FullSigmaTest2) {
-
-    std::streambuf *defaultCout;
-    std::streambuf *defaultCerr;
-    std::ostringstream debugOutput;
-    defaultCout = std::cout.rdbuf(debugOutput.rdbuf());
-    defaultCerr = std::cerr.rdbuf(debugOutput.rdbuf());
-
+    OpalTestUtilities::SilenceTest silencer(true);
     char inputFileName[] = "GaussDistributionTest.in";
     std::string input = "OPTION, ECHO=FALSE;\n"
         "OPTION, CZERO=FALSE;\n"
@@ -214,8 +196,6 @@ TEST(GaussTest, FullSigmaTest2) {
         is = new FileStream(inputFileName);
     } catch(...) {
         is = 0;
-        std::cout.rdbuf(defaultCout);;
-        std::cerr.rdbuf(defaultCerr);;
         throw new OpalException("FullSigmaTest", "Could not read string");
     }
 
@@ -224,8 +204,6 @@ TEST(GaussTest, FullSigmaTest2) {
         try {
             parser->run(is);
         } catch (...) {
-            std::cout.rdbuf(defaultCout);;
-            std::cerr.rdbuf(defaultCerr);;
             throw new OpalException("FullSigmaTest", "Could not parse input");
         }
     }
@@ -235,8 +213,6 @@ TEST(GaussTest, FullSigmaTest2) {
         distObj = OPAL->find("DIST1");
     } catch(...) {
         distObj = 0;
-        std::cout.rdbuf(defaultCout);;
-        std::cerr.rdbuf(defaultCerr);;
         throw new OpalException("FullSigmaTest", "Could not find distribution");
     }
 
@@ -267,18 +243,15 @@ TEST(GaussTest, FullSigmaTest2) {
         const double expectedR61 = 1.362;
         const double expectedR62 = -0.2685;
 
-        EXPECT_LT(std::abs(expectedR11 - R11),  0.05 * expectedR11) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR21 - R21), -0.05 * expectedR21) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR22 - R22),  0.05 * expectedR22) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR51 - R51),  0.05 * expectedR51) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR52 - R52),  0.05 * expectedR52) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR61 - R61),  0.05 * expectedR61) << ::burnAfterReading(debugOutput);
-        EXPECT_LT(std::abs(expectedR62 - R62), -0.05 * expectedR62) << ::burnAfterReading(debugOutput);
+        EXPECT_LT(std::abs(expectedR11 - R11),  0.05 * expectedR11);
+        EXPECT_LT(std::abs(expectedR21 - R21), -0.05 * expectedR21);
+        EXPECT_LT(std::abs(expectedR22 - R22),  0.05 * expectedR22);
+        EXPECT_LT(std::abs(expectedR51 - R51),  0.05 * expectedR51);
+        EXPECT_LT(std::abs(expectedR52 - R52),  0.05 * expectedR52);
+        EXPECT_LT(std::abs(expectedR61 - R61),  0.05 * expectedR61);
+        EXPECT_LT(std::abs(expectedR62 - R62), -0.05 * expectedR62);
     }
 
-    std::cout.rdbuf(defaultCout);;
-    std::cerr.rdbuf(defaultCerr);;
-
     OpalData::deleteInstance();
     delete parser;
     delete gmsg;
@@ -286,4 +259,4 @@ TEST(GaussTest, FullSigmaTest2) {
     delete[] arg;
 
     std::remove(inputFileName);
-}
\ No newline at end of file
+}
diff --git a/tests/opal_test_utilities/CMakeLists.txt b/tests/opal_test_utilities/CMakeLists.txt
new file mode 100644
index 000000000..7e11c477c
--- /dev/null
+++ b/tests/opal_test_utilities/CMakeLists.txt
@@ -0,0 +1,8 @@
+set (_SRCS
+)
+
+include_directories (
+  ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+add_sources(${_SRCS})
diff --git a/tests/opal_test_utilities/SilenceTest.h b/tests/opal_test_utilities/SilenceTest.h
new file mode 100644
index 000000000..6d3f5bb2f
--- /dev/null
+++ b/tests/opal_test_utilities/SilenceTest.h
@@ -0,0 +1,56 @@
+/*
+ *  Copyright (c) 2012, Chris Rogers
+ *  All rights reserved.
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions are met:
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *  2. Redistributions in binary form must reproduce the above copyright notice,
+ *     this list of conditions and the following disclaimer in the documentation
+ *     and/or other materials provided with the distribution.
+ *  3. Neither the name of STFC nor the names of its contributors may be used to
+ *     endorse or promote products derived from this software without specific
+ *     prior written permission.
+ *
+ *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ *  POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef OPALTESTUTILITIES_SILENCETEST_H_
+
+namespace OpalTestUtilities {
+/** Shutup test output */
+class SilenceTest {
+  public:
+    SilenceTest(bool willSilence) {
+        if (willSilence) {
+            _defaultCout = std::cout.rdbuf(_debugOutput.rdbuf());
+            _defaultCerr = std::cerr.rdbuf(_debugOutput.rdbuf());
+        }
+    }
+
+    ~SilenceTest() { // return buffer to normal on delete
+        std::cout.rdbuf(_defaultCout);
+        std::cerr.rdbuf(_defaultCerr);
+    }
+
+  private:
+    SilenceTest(); // disable default ctor
+    SilenceTest(const SilenceTest& test); // disable default copy ctor
+
+    std::ostringstream _debugOutput;
+    std::streambuf *_defaultCout;
+    std::streambuf *_defaultCerr;
+};
+}
+
+#endif //OPALTESTUTILITIES_SILENCETEST_H_
-- 
GitLab