diff --git a/Framework/Crystal/inc/MantidCrystal/PeakAlgorithmHelpers.h b/Framework/Crystal/inc/MantidCrystal/PeakAlgorithmHelpers.h
index f6460157dfcbcacda8fe5147477821f4e9830103..fe6cd38fb9c45251f571b7de98a2545f664e842d 100644
--- a/Framework/Crystal/inc/MantidCrystal/PeakAlgorithmHelpers.h
+++ b/Framework/Crystal/inc/MantidCrystal/PeakAlgorithmHelpers.h
@@ -14,6 +14,12 @@ class V3D;
 
 namespace Mantid::Crystal {
 
+/// return -1 if convention is "Crystallography" and 1 otherwise.
+double qConventionFactor(const std::string &convention);
+
+/// convenience overload to pull the convention from the config service
+double qConventionFactor();
+
 /// Tie together a modulated peak number with its offset
 using MNPOffset = std::tuple<double, double, double, Kernel::V3D>;
 
diff --git a/Framework/Crystal/src/IndexPeaks.cpp b/Framework/Crystal/src/IndexPeaks.cpp
index 716cd1a5053703cb9a64dae92f76170063016262..ed20b30a09a403b06887fcbc1e912b041999b0b4 100644
--- a/Framework/Crystal/src/IndexPeaks.cpp
+++ b/Framework/Crystal/src/IndexPeaks.cpp
@@ -432,9 +432,9 @@ void IndexPeaks::exec() {
     auto &lattice = args.workspace->mutableSample().getOrientedLattice();
     lattice.setMaxOrder(args.satellites.maxOrder);
     lattice.setCrossTerm(args.satellites.crossTerms);
-    for (auto i = 0u; i < 3; ++i) {
-      lattice.setModVec1(args.satellites.modVectors[i]);
-    }
+    lattice.setModVec1(args.satellites.modVectors[0]);
+    lattice.setModVec2(args.satellites.modVectors[1]);
+    lattice.setModVec3(args.satellites.modVectors[2]);
   }
 
   CombinedIndexingStats indexingInfo;
diff --git a/Framework/Crystal/src/PeakAlgorithmHelpers.cpp b/Framework/Crystal/src/PeakAlgorithmHelpers.cpp
index 3081069b74f07996487b9297a5fe171369867118..88cd01b92302da27e44c3de5a3b3855af93d2dbc 100644
--- a/Framework/Crystal/src/PeakAlgorithmHelpers.cpp
+++ b/Framework/Crystal/src/PeakAlgorithmHelpers.cpp
@@ -8,16 +8,35 @@
 #include "MantidKernel/ArrayLengthValidator.h"
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/ConfigService.h"
 #include "MantidKernel/V3D.h"
 
 using Mantid::Kernel::ArrayLengthValidator;
 using Mantid::Kernel::ArrayProperty;
 using Mantid::Kernel::BoundedValidator;
+using Mantid::Kernel::ConfigService;
 using Mantid::Kernel::Direction;
 using Mantid::Kernel::V3D;
 
 namespace Mantid::Crystal {
 
+/**
+ * @return -1 if Q.convention==Crystallography, else return 1.0
+ */
+double qConventionFactor(const std::string &convention) {
+  if (convention == "Crystallography") {
+    return -1.0;
+  }
+  return 1.0;
+}
+
+/**
+ * @return qConventionFactor for Q.convention setting in ConfigService
+ */
+double qConventionFactor() {
+  return qConventionFactor(ConfigService::Instance().getString("Q.convention"));
+}
+
 /**
  * Append the common set of properties that relate to modulation vectors
  * to the given algorithm
diff --git a/Framework/Crystal/src/PredictFractionalPeaks.cpp b/Framework/Crystal/src/PredictFractionalPeaks.cpp
index 0b7600055792f25418ee8e641f188d926edbd00a..3d3827fe7457242e6f99b6a758d0ad72606f7f0d 100644
--- a/Framework/Crystal/src/PredictFractionalPeaks.cpp
+++ b/Framework/Crystal/src/PredictFractionalPeaks.cpp
@@ -229,6 +229,7 @@ IPeaksWorkspace_sptr predictFractionalPeaks(
   using PeakHash = std::array<int, 4>;
   std::vector<PeakHash> alreadyDonePeaks;
 
+  const auto qConvention{Mantid::Crystal::qConventionFactor()};
   V3D currentHKL;
   DblMatrix gonioMatrix;
   int runNumber{0};
@@ -236,8 +237,12 @@ IPeaksWorkspace_sptr predictFractionalPeaks(
   auto progressReporter = searchStrategy.createProgressReporter(alg);
   while (true) {
     for (const auto &mnpOffset : offsets) {
-      const V3D candidateHKL{currentHKL + std::get<3>(mnpOffset)};
-      const V3D qLab = (gonioMatrix * UB * candidateHKL) * 2 * M_PI;
+      const V3D currentIntHKL{std::round(currentHKL[0]),
+                              std::round(currentHKL[1]),
+                              std::round(currentHKL[2])};
+      const V3D candidateHKL{currentIntHKL + std::get<3>(mnpOffset)};
+      const V3D qLab =
+          (gonioMatrix * UB * candidateHKL) * 2 * M_PI * qConvention;
       if (qLab[2] <= 0)
         continue;
 
@@ -267,7 +272,8 @@ IPeaksWorkspace_sptr predictFractionalPeaks(
       else
         continue;
 
-      peak->setHKL(candidateHKL);
+      peak->setHKL(candidateHKL * qConvention);
+      peak->setIntHKL(currentIntHKL * qConvention);
       const double m{std::get<0>(mnpOffset)}, n{std::get<1>(mnpOffset)},
           p{std::get<2>(mnpOffset)};
       if (fabs(m) > 0. || fabs(n) > 0. || fabs(p) > 0.)
diff --git a/Framework/Crystal/src/PredictPeaks.cpp b/Framework/Crystal/src/PredictPeaks.cpp
index 17558ae9e9c06fbaf5d3f78aac2adc780cd37b60..1b721283210b2536072dca806fd0cdfb14916b95 100644
--- a/Framework/Crystal/src/PredictPeaks.cpp
+++ b/Framework/Crystal/src/PredictPeaks.cpp
@@ -9,6 +9,7 @@
 #include "MantidAPI/MatrixWorkspace.h"
 #include "MantidAPI/Run.h"
 #include "MantidAPI/Sample.h"
+#include "MantidCrystal/PeakAlgorithmHelpers.h"
 #include "MantidGeometry/Crystal/BasicHKLFilters.h"
 #include "MantidGeometry/Crystal/EdgePixel.h"
 #include "MantidGeometry/Crystal/HKLFilterWavelength.h"
@@ -38,24 +39,11 @@ using namespace Mantid::DataObjects;
 using namespace Mantid::Geometry;
 using namespace Mantid::Kernel;
 
-namespace {
-/// Small helper function that return -1 if convention
-/// is "Crystallography" and 1 otherwise.
-double get_factor_for_q_convention(const std::string &convention) {
-  if (convention == "Crystallography") {
-    return -1.0;
-  }
-
-  return 1.0;
-}
-} // namespace
-
 /** Constructor
  */
 PredictPeaks::PredictPeaks()
     : m_runNumber(-1), m_inst(), m_pw(), m_sfCalculator(),
-      m_qConventionFactor(get_factor_for_q_convention(
-          ConfigService::Instance().getString("Q.convention"))) {
+      m_qConventionFactor(qConventionFactor()) {
   m_refConds = getAllReflectionConditions();
 }
 
@@ -510,7 +498,7 @@ void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(
    * for the convention stored in the workspace.
    */
   double peaks_q_convention_factor =
-      get_factor_for_q_convention(peaksWorkspace->getConvention());
+      qConventionFactor(peaksWorkspace->getConvention());
 
   for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
     IPeak &p = peaksWorkspace->getPeak(i);
diff --git a/Framework/Crystal/src/PredictSatellitePeaks.cpp b/Framework/Crystal/src/PredictSatellitePeaks.cpp
index 2ab7870ba45aa5fc5d9505f7d2bbc4a6551e6e17..47f15b8df01bbd2380c202b61ebd8eb5082e903d 100644
--- a/Framework/Crystal/src/PredictSatellitePeaks.cpp
+++ b/Framework/Crystal/src/PredictSatellitePeaks.cpp
@@ -9,6 +9,7 @@
 #include "MantidAPI/Run.h"
 #include "MantidAPI/Sample.h"
 #include "MantidAPI/WorkspaceFactory.h"
+#include "MantidCrystal/PeakAlgorithmHelpers.h"
 #include "MantidDataObjects/PeaksWorkspace.h"
 #include "MantidGeometry/Crystal/BasicHKLFilters.h"
 #include "MantidGeometry/Crystal/HKLGenerator.h"
@@ -30,23 +31,11 @@ namespace Crystal {
 
 DECLARE_ALGORITHM(PredictSatellitePeaks)
 
-namespace {
-/// Small helper function that return -1 if convention
-/// is "Crystallography" and 1 otherwise.
-double get_factor_for_q_convention(const std::string &convention) {
-  if (convention == "Crystallography") {
-    return -1.0;
-  }
-  return 1.0;
-}
-} // namespace
-
 /** Constructor
  */
 
 PredictSatellitePeaks::PredictSatellitePeaks()
-    : m_qConventionFactor(get_factor_for_q_convention(
-          ConfigService::Instance().getString("Q.convention"))) {}
+    : m_qConventionFactor(qConventionFactor()) {}
 
 /// Initialise the properties
 void PredictSatellitePeaks::init() {
diff --git a/Framework/Crystal/test/IndexPeaksTest.h b/Framework/Crystal/test/IndexPeaksTest.h
index d27ea00df0cbc2173cab035bc41dbe6f13d9ee3e..a516b870883bcfe6489810918ea82334c1b8dffa 100644
--- a/Framework/Crystal/test/IndexPeaksTest.h
+++ b/Framework/Crystal/test/IndexPeaksTest.h
@@ -448,9 +448,9 @@ public:
     TS_ASSERT_DELTA(V3D(0.333, 0.667, -0.333).norm(),
                     lattice.getModVec(0).norm(), 1e-8)
     TS_ASSERT_DELTA(V3D(-0.333, 0.667, 0.333).norm(),
-                    lattice.getModVec(0).norm(), 1e-8)
+                    lattice.getModVec(1).norm(), 1e-8)
     TS_ASSERT_DELTA(V3D(0.333, -0.667, 0.333).norm(),
-                    lattice.getModVec(0).norm(), 1e-8)
+                    lattice.getModVec(2).norm(), 1e-8)
     TS_ASSERT_EQUALS(true, lattice.getCrossTerm())
   }
 
diff --git a/Framework/Crystal/test/PredictFractionalPeaksTest.h b/Framework/Crystal/test/PredictFractionalPeaksTest.h
index f0a179e4c1246e54450d1744c3dc9ef4e4d2b798..c975ef339a19d9b59af062aa99f494b780127e1a 100644
--- a/Framework/Crystal/test/PredictFractionalPeaksTest.h
+++ b/Framework/Crystal/test/PredictFractionalPeaksTest.h
@@ -263,6 +263,10 @@ public:
     TS_ASSERT_DELTA(peak0.getH(), -4.5, .0001)
     TS_ASSERT_DELTA(peak0.getK(), 7.0, .0001)
     TS_ASSERT_DELTA(peak0.getL(), -4.5, .0001)
+    const auto mainHKL0 = peak0.getIntHKL();
+    TS_ASSERT_DELTA(mainHKL0[0], -5, .0001)
+    TS_ASSERT_DELTA(mainHKL0[1], 7.0, .0001)
+    TS_ASSERT_DELTA(mainHKL0[2], -4, .0001)
     TS_ASSERT_EQUALS(peak0.getDetectorID(), 1129591)
     const auto mnp0 = peak0.getIntMNP();
     TS_ASSERT_DELTA(mnp0.X(), -1., 1e-08)
@@ -273,6 +277,10 @@ public:
     TS_ASSERT_DELTA(peak34.getH(), -7, .0001)
     TS_ASSERT_DELTA(peak34.getK(), 7.5, .0001)
     TS_ASSERT_DELTA(peak34.getL(), -2.5, .0001)
+    const auto mainHKL34 = peak34.getIntHKL();
+    TS_ASSERT_DELTA(mainHKL34[0], -7, .0001)
+    TS_ASSERT_DELTA(mainHKL34[1], 7.0, .0001)
+    TS_ASSERT_DELTA(mainHKL34[2], -3, .0001)
     TS_ASSERT_EQUALS(peak34.getDetectorID(), 1812163)
     const auto mnp34 = peak34.getIntMNP();
     TS_ASSERT_DELTA(mnp34.X(), 0., 1e-08)
diff --git a/Framework/PythonInterface/mantid/api/src/Exports/IPeaksWorkspace.cpp b/Framework/PythonInterface/mantid/api/src/Exports/IPeaksWorkspace.cpp
index c3b0634a4014772eb103de74bddae197b3abe3db..4b75069ad94043176c2be4616df327b7026b1644 100644
--- a/Framework/PythonInterface/mantid/api/src/Exports/IPeaksWorkspace.cpp
+++ b/Framework/PythonInterface/mantid/api/src/Exports/IPeaksWorkspace.cpp
@@ -13,6 +13,7 @@
 #include <boost/none.hpp>
 #include <boost/optional.hpp>
 #include <boost/python/class.hpp>
+#include <boost/python/iterator.hpp>
 #include <boost/python/manage_new_object.hpp>
 #include <boost/python/return_internal_reference.hpp>
 
@@ -206,8 +207,49 @@ void setCell(IPeaksWorkspace &self, const object &col_or_row,
   PeakWorkspaceTableAdaptor tableMap{self};
   tableMap.setProperty(columnName, rowIndex, value);
 }
+
+/// Internal helper to support iteration on PeaksWorkspace in Python
+struct IPeaksWorkspaceIterator {
+  explicit IPeaksWorkspaceIterator(IPeaksWorkspace *const workspace)
+      : m_workspace{workspace}, m_numPeaks{workspace->getNumberPeaks()},
+        m_rowIndex{-1} {
+    assert(workspace);
+  }
+  IPeak *next() {
+    ++m_rowIndex;
+    if (m_rowIndex >= m_numPeaks) {
+      objects::stop_iteration_error();
+    }
+    return m_workspace->getPeakPtr(m_rowIndex);
+  }
+
+private:
+  IPeaksWorkspace *const m_workspace;
+  const int m_numPeaks;
+  int m_rowIndex;
+};
+
+// Create an iterator from the given workspace
+IPeaksWorkspaceIterator makePyIterator(IPeaksWorkspace &self) {
+  return IPeaksWorkspaceIterator(&self);
+}
+
 } // namespace
 
+void export_IPeaksWorkspaceIterator() {
+  class_<IPeaksWorkspaceIterator>("IPeaksWorkspaceIterator", no_init)
+      .def(
+#if PY_VERSION_HEX >= 0x03000000
+          "__next__"
+#else
+          "next"
+#endif
+          ,
+          &IPeaksWorkspaceIterator::next,
+          return_value_policy<reference_existing_object>())
+      .def("__iter__", objects::identity_function());
+}
+
 void export_IPeaksWorkspace() {
   // IPeaksWorkspace class
   class_<IPeaksWorkspace, bases<ITableWorkspace, ExperimentInfo>,
@@ -245,8 +287,8 @@ void export_IPeaksWorkspace() {
             arg("value")),
            "Sets the value of a given cell. If the row_or_column argument is a "
            "number then it is interpreted as a row otherwise it "
-           "is interpreted as a column name.");
-
+           "is interpreted as a column name.")
+      .def("__iter__", makePyIterator);
   //-------------------------------------------------------------------------------------------------
 
   RegisterWorkspacePtrToPython<IPeaksWorkspace>();
diff --git a/Framework/PythonInterface/mantid/api/src/Exports/Sample.cpp b/Framework/PythonInterface/mantid/api/src/Exports/Sample.cpp
index e0e1c66690b1d8dfa30c6aa9bab76b062d6ec910..594233d551259ef099b840fa2d3d14ae8d30896a 100644
--- a/Framework/PythonInterface/mantid/api/src/Exports/Sample.cpp
+++ b/Framework/PythonInterface/mantid/api/src/Exports/Sample.cpp
@@ -41,6 +41,8 @@ void export_Sample() {
       .def("hasOrientedLattice", &Sample::hasOrientedLattice, arg("self"),
            "Returns True if this sample has an oriented lattice, false "
            "otherwise")
+      .def("clearOrientedLattice", &Sample::clearOrientedLattice, arg("self"),
+           "Clears any attached lattice information")
       .def("getCrystalStructure", &Sample::getCrystalStructure, arg("self"),
            return_value_policy<reference_existing_object>(),
            "Get the crystal structure for this sample")
diff --git a/Framework/PythonInterface/mantid/geometry/src/Exports/UnitCell.cpp b/Framework/PythonInterface/mantid/geometry/src/Exports/UnitCell.cpp
index 77d7cdd9fc84530ef8a73193710ba48f5232e844..c56ad2f6df791f7d1ee56c71314da41cf9126986 100644
--- a/Framework/PythonInterface/mantid/geometry/src/Exports/UnitCell.cpp
+++ b/Framework/PythonInterface/mantid/geometry/src/Exports/UnitCell.cpp
@@ -325,6 +325,20 @@ void export_UnitCell() {
            "Set the error in the :math:`\\gamma` angle of the unit cell using "
            "the "
            "``Unit`` parameter.")
+      .def("setModVec1",
+           (void (UnitCell::*)(const V3D &)) & UnitCell::setModVec1,
+           (arg("self"), arg("vec")),
+           "Set the first modulated structure vector")
+      .def("setModVec2",
+           (void (UnitCell::*)(const V3D &)) & UnitCell::setModVec2,
+           (arg("self"), arg("vec")),
+           "Set the second modulated structure vector")
+      .def("setModVec3",
+           (void (UnitCell::*)(const V3D &)) & UnitCell::setModVec3,
+           (arg("self"), arg("vec")),
+           "Set the third modulated structure vector")
+      .def("setMaxOrder", &UnitCell::setMaxOrder,
+           "Set the maximum order of modulated vectors searched")
       .def("volume", (double (UnitCell::*)() const) & UnitCell::volume,
            arg("self"),
            "Return the volume of the unit cell (in :math:`\\rm{\\AA}{^3}`)")
@@ -354,6 +368,8 @@ void export_UnitCell() {
       .def("getMaxOrder", &UnitCell::getMaxOrder, arg("self"),
            "Returns the number of modulation vectors. This will return an "
            "int.")
+      .def("getModVec", &UnitCell::getModVec, (arg("self"), arg("i")),
+           "Returns the ith modulation vector")
       .def("recalculateFromGstar", &recalculateFromGstar,
            (arg("self"), arg("NewGstar")),
            "Recalculate the unit cell parameters from a metric tensor. This "
diff --git a/Framework/PythonInterface/plugins/algorithms/IndexSatellitePeaks.py b/Framework/PythonInterface/plugins/algorithms/IndexSatellitePeaks.py
index eb7ca8c707be91dad7de81303749faa8be046f85..a99451e8d9693b738199c8c84aa2213d9f0fded8 100644
--- a/Framework/PythonInterface/plugins/algorithms/IndexSatellitePeaks.py
+++ b/Framework/PythonInterface/plugins/algorithms/IndexSatellitePeaks.py
@@ -8,9 +8,8 @@ import numpy as np
 from scipy.spatial import KDTree
 import fractional_indexing as indexing
 
-from mantid.kernel import Direction
-from mantid.api import (IPeaksWorkspaceProperty,
-                        ITableWorkspaceProperty, PythonAlgorithm, AlgorithmFactory)
+from mantid.kernel import Direction, V3D
+from mantid.api import (IPeaksWorkspaceProperty, PythonAlgorithm, AlgorithmFactory)
 import mantid.simpleapi as api
 
 
@@ -40,7 +39,7 @@ class IndexSatellitePeaks(PythonAlgorithm):
                              doc="Positions of satellite peaks with fractional \
                              HKL coordinates")
 
-        self.declareProperty(ITableWorkspaceProperty(name="OutputWorkspace",
+        self.declareProperty(IPeaksWorkspaceProperty(name="OutputWorkspace",
                                                      defaultValue="",
                                                      direction=Direction.Output),
                              doc="The indexed satellite peaks. This will be a  \
@@ -100,56 +99,34 @@ class IndexSatellitePeaks(PythonAlgorithm):
             distance, index = peak_map.query(q, k=1)
             hklm[i, 3:] = indices[index]
 
-        indexed = self.create_indexed_workspace(satellites, ndim, hklm)
+        indexed = self.create_indexed_peaksworkspace(satellites, qs, hklm)
         self.setProperty("OutputWorkspace", indexed)
 
-    def create_indexed_workspace(self, fractional_peaks, ndim, hklm):
-        """Create a TableWorkepace that contains indexed peak data.
-
-        This produces a TableWorkepace that looks like a PeaksWorkspace but
-        with the additional index columns included. In future releases support
-        for indexing should be added to the PeaksWorkspace data type itself.
+    def create_indexed_peaksworkspace(self, fractional_peaks, qs, hklm):
+        """Create a PeaksWorkspace that contains indexed peak data.
 
         :param fractional_peaks: the peaks workspace containing peaks with
             fractional HKL values.
-        :param ndim: the number of additional indexing columns to add.
-        :param hklm: the new higher dimensional miller indicies to add.
-        :returns: a table workspace with the indexed peak data
+        :param qs: The set of modulation vectors determined
+        :param hklm: the new higher dimensional miller indices to add.
+        :returns: a peaks workspace with the indexed peak data
         """
-        # Create table with the number of columns we need
-        types = ['int', 'long64', 'double', 'double', 'double', 'double',  'double', 'double',
-                 'double', 'double', 'double', 'float', 'str', 'float', 'float', 'V3D', 'V3D']
-        name = self.getPropertyValue("OutputWorkspace")
-        indexed = api.CreateEmptyTableWorkspace(name)
-        names = fractional_peaks.getColumnNames()
-
-        # Insert the extra columns for the addtional indicies
-        for i in range(ndim - 3):
-            names.insert(5 + i, 'm{}'.format(i + 1))
-            types.insert(5 + i, 'double')
-
-        names = np.array(names)
-        types = np.array(types)
-
-        # Create columns in the table workspace
-        for name, column_type in zip(names, types):
-            indexed.addColumn(column_type, name)
-
-        # Copy all columns from original workspace, ignoring HKLs
-        column_data = []
-        idx = np.arange(0, names.size)
-        hkl_mask = (idx < 2) | (idx > 4 + (ndim - 3))
-        for name in names[hkl_mask]:
-            column_data.append(fractional_peaks.column(name))
-
-        # Insert the addtional HKL columns into the data
-        for i, col in enumerate(hklm.T.tolist()):
-            column_data.insert(i + 2, col)
-
-        # Insert the columns into the table workspace
-        for i in range(fractional_peaks.rowCount()):
-            row = [column_data[j][i] for j in range(indexed.columnCount())]
-            indexed.addRow(row)
+        # pad to 6 columns so we can assume a (hkl) (mnp) layout
+        hklm = np.pad(hklm, pad_width=(0, 6 - hklm.shape[1]), mode='constant',
+                      constant_values=0)
+        indexed = api.CloneWorkspace(fractional_peaks, StoreInADS=False)
+        # save modulation vectors. ensure qs has 3 rows
+        qs = np.pad(qs, pad_width=((0, 3 - qs.shape[0]), (0, 0)), mode='constant',
+                    constant_values=0)
+        lattice = fractional_peaks.sample().getOrientedLattice()
+        lattice.setModVec1(V3D(*qs[0]))
+        lattice.setModVec2(V3D(*qs[1]))
+        lattice.setModVec3(V3D(*qs[2]))
+        # save indices
+        for row, peak in enumerate(indexed):
+            row_indices = hklm[row]
+            peak.setHKL(*row_indices[:3])
+            peak.setIntMNP(V3D(*row_indices[3:]))
 
         return indexed
 
diff --git a/Framework/PythonInterface/plugins/algorithms/LinkedUBs.py b/Framework/PythonInterface/plugins/algorithms/LinkedUBs.py
index f4065bd26aaebd349836fc696feb49f445503956..8af87f90fd7cf60ee18b74d57c19f2d402f4d231 100644
--- a/Framework/PythonInterface/plugins/algorithms/LinkedUBs.py
+++ b/Framework/PythonInterface/plugins/algorithms/LinkedUBs.py
@@ -310,7 +310,7 @@ class LinkedUBs(DataProcessorAlgorithm):
 
             # get the indexing list that sorts dspacing from largest to
             # smallest
-            hkls = np.array([[p['h'], p['k'], p['l']] for p in predictor])
+            hkls = np.array([[p.getH(), p.getK(), p.getL()] for p in predictor])
             idx = dspacings_predicted.argsort()[::-1]
             HKL_predicted = hkls[idx, :]
 
diff --git a/Framework/PythonInterface/plugins/algorithms/SaveReflections.py b/Framework/PythonInterface/plugins/algorithms/SaveReflections.py
index 265743f523cd280445a818967269f4e57c642a11..73c0c3f6227df1b5ad027f7239a22c507843fbdd 100644
--- a/Framework/PythonInterface/plugins/algorithms/SaveReflections.py
+++ b/Framework/PythonInterface/plugins/algorithms/SaveReflections.py
@@ -4,47 +4,79 @@
 #     NScD Oak Ridge National Laboratory, European Spallation Source
 #     & Institut Laue - Langevin
 # SPDX - License - Identifier: GPL - 3.0 +
-import re
+from __future__ import (absolute_import, division, unicode_literals)
+
+import os.path as osp
 import numpy as np
-from mantid.api import AlgorithmFactory, FileProperty, FileAction, PythonAlgorithm, ITableWorkspaceProperty, IPeaksWorkspace
+from mantid.api import (AlgorithmFactory, FileProperty, FileAction, IPeaksWorkspaceProperty, PythonAlgorithm)
 from mantid.kernel import StringListValidator, Direction
-from mantid.simpleapi import SaveHKL
+from mantid.py3compat.enum import Enum
+
+
+def num_modulation_vectors(workspace):
+    """Check if the workspace has any modulation vectors stored and
+    return the number if any
 
-# List of file format names supported by this algorithm
-SUPPORTED_FORMATS = ["Fullprof", "GSAS", "Jana", "SHELX"]
+    :params: workspace :: the workspace to check
+    :returns: The number of stored modulation vectors
+    """
+    sample = workspace.sample()
+    if sample.hasOrientedLattice():
+        lattice = workspace.sample().getOrientedLattice()
+        count = 0
+        for i in range(3):
+            vec = lattice.getModVec(i)
+            if abs(vec.X()) > 0 or abs(vec.Y()) or abs(vec.Z()):
+                count += 1
+        return count
+    else:
+        return 0
+
+
+def get_two_theta(dspacing, wavelength):
+    """Get the two theta value for this peak.
+
+    This is just Bragg's law relating wavelength to scattering angle.
+    :param dspacing: DSpacing of a peak
+    :param wavelength: Wavelength of a peak
+    :returns: the scattering angle for the peak.
+    """
+    theta = 2. * np.arcsin(0.5 * (wavelength / dspacing))
+    return np.rad2deg(theta)
 
 
 def has_modulated_indexing(workspace):
-    """Check if this workspace has more than 3 indicies
+    """Check if this workspace has more than 3 indices
 
     :params: workspace :: the workspace to check
-    :returns: True if the workspace > 3 indicies else False
+    :returns: True if the workspace > 3 indices else False
     """
-    return num_additional_indicies(workspace) > 0
+    return num_modulation_vectors(workspace) > 0
 
 
-def num_additional_indicies(workspace):
-    """Check if this workspace has more than 3 indicies
+def modulation_indices(peak, num_mod_vec):
+    """
+    Gather non-zero modulated structure indices from a peak
 
-    :params: workspace :: the workspace count indicies in
-    :returns: the number of additional indicies present in the workspace
+    :param workspace: A single Peak
+    :param num_mod_vec: The number of modulation vectors set on the workspace
+    :return: A list of the modulation indices
     """
-    return len(get_additional_index_names(workspace))
+    mnp = peak.getIntMNP()
+    return [mnp[i] for i in range(num_mod_vec)]
 
 
 def get_additional_index_names(workspace):
-    """Get the names of the additional indicies to export
+    """Get the names of the additional indices to export
 
     :params: workspace :: the workspace to get column names from
     :returns: the names of any additional columns in the workspace
     """
-    pattern = re.compile("m[1-9]+")
-    names = workspace.getColumnNames()
-    return list(filter(pattern.match, names))
+    num_mod_vec = num_modulation_vectors(workspace)
+    return ["m{}".format(i + 1) for i in range(num_mod_vec)]
 
 
 class SaveReflections(PythonAlgorithm):
-
     def category(self):
         return "DataHandling\\Text;Crystal\\DataHandling"
 
@@ -54,10 +86,11 @@ class SaveReflections(PythonAlgorithm):
     def PyInit(self):
         """Initilize the algorithms properties"""
 
-        self.declareProperty(ITableWorkspaceProperty("InputWorkspace", '', Direction.Input),
+        self.declareProperty(IPeaksWorkspaceProperty("InputWorkspace", '', Direction.Input),
                              doc="The name of the peaks worksapce to save")
 
-        self.declareProperty(FileProperty("Filename", "",
+        self.declareProperty(FileProperty("Filename",
+                                          "",
                                           action=FileAction.Save,
                                           direction=Direction.Input),
                              doc="File with the data from a phonon calculation.")
@@ -65,36 +98,25 @@ class SaveReflections(PythonAlgorithm):
         self.declareProperty(name="Format",
                              direction=Direction.Input,
                              defaultValue="Fullprof",
-                             validator=StringListValidator(SUPPORTED_FORMATS),
+                             validator=StringListValidator(dir(ReflectionFormat)),
                              doc="The output format to export reflections to")
 
+        self.declareProperty(name="SplitFiles",
+                             defaultValue=False,
+                             direction=Direction.Input,
+                             doc="If True save separate files with only the peaks associated"
+                             "with a single modulation vector in a single file. Only "
+                             "applies to JANA format.")
+
     def PyExec(self):
         """Execute the algorithm"""
         workspace = self.getProperty("InputWorkspace").value
-        output_format = self.getPropertyValue("Format")
+        output_format = ReflectionFormat[self.getPropertyValue("Format")]
         file_name = self.getPropertyValue("Filename")
+        split_files = self.getProperty("SplitFiles").value
 
-        file_writer  = self.choose_format(output_format)
-        file_writer(file_name, workspace)
-
-    def choose_format(self, output_format):
-        """Choose the function to use to write out data for this format
-
-        :param output_format: the format to use to output refelctions as. Options are
-            "Fullprof", "GSAS", "Jana", and "SHELX".
+        FORMAT_MAP[output_format]()(file_name, workspace, split_files)
 
-        :returns: file format to use for saving reflections to an ASCII file.
-        """
-
-        if output_format == "Fullprof":
-            return FullprofFormat()
-        elif output_format == "Jana":
-            return JanaFormat()
-        elif output_format == "GSAS" or output_format == "SHELX":
-            return SaveHKLFormat()
-        else:
-            raise ValueError("Unexpected file format {}. Format should be one of {}."
-                             .format(output_format, ", ".join(SUPPORTED_FORMATS)))
 
 # ------------------------------------------------------------------------------------------------------
 
@@ -106,12 +128,12 @@ class FullprofFormat(object):
     This is a 7 columns file format consisting of H, K, L, instensity,
     sigma, crystal domain, and wavelength.
     """
-
-    def __call__(self, file_name, workspace):
+    def __call__(self, file_name, workspace, split_files):
         """Write a PeaksWorkspace to an ASCII file using this formatter.
 
         :param file_name: the file name to output data to.
         :param workspace: the PeaksWorkspace to write to file.
+        :param _: Ignored parameter for compatability with other savers
         """
         with open(file_name, 'w') as f_handle:
             self.write_header(f_handle, workspace)
@@ -123,11 +145,11 @@ class FullprofFormat(object):
         :param f_handle: handle to the file to write to.
         :param workspace: the PeaksWorkspace to save to file.
         """
-        num_hkl = 3+num_additional_indicies(workspace)
+        num_hkl = 3 + num_modulation_vectors(workspace)
         f_handle.write(workspace.getTitle())
         f_handle.write("({}i4,2f12.2,i5,4f10.4)\n".format(num_hkl))
         f_handle.write("  0 0 0\n")
-        names =  "".join(["  {}".format(name) for name in get_additional_index_names(workspace)])
+        names = "".join(["  {}".format(name) for name in get_additional_index_names(workspace)])
         f_handle.write("#  h   k   l{}      Fsqr       s(Fsqr)   Cod   Lambda\n".format(names))
 
     def write_peaks(self, f_handle, workspace):
@@ -136,18 +158,19 @@ class FullprofFormat(object):
         :param f_handle: handle to the file to write to.
         :param workspace: the PeaksWorkspace to save to file.
         """
+        num_mod_vec = num_modulation_vectors(workspace)
         for i, peak in enumerate(workspace):
-            data = [peak['h'],peak['k'],peak['l']]
-            data.extend([peak[name] for name in get_additional_index_names(workspace)])
-
+            data = [peak.getH(), peak.getK(), peak.getL()]
+            data.extend(modulation_indices(peak, num_mod_vec))
             hkls = "".join(["{:>4.0f}".format(item) for item in data])
 
-            data = (peak['Intens'],peak['SigInt'],i+1,peak['Wavelength'])
+            data = (peak.getIntensity(), peak.getSigmaIntensity(), i + 1, peak.getWavelength())
             line = "{:>12.2f}{:>12.2f}{:>5.0f}{:>10.4f}\n".format(*data)
             line = "".join([hkls, line])
 
             f_handle.write(line)
 
+
 # ------------------------------------------------------------------------------------------------------
 
 
@@ -159,80 +182,170 @@ class JanaFormat(object):
     crystal domain, wavelength, 2*theta, transmission, absorption weighted path length (Tbar),
     and thermal diffuse scattering correction (TDS).
 
-    Currently the last three columns are shard coded to 1.0, 0.0, and 0.0 respectively.
+    Currently the last three columns are hard coded to 1.0, 0.0, and 0.0 respectively.
     """
-
-    def __call__(self, file_name, workspace):
-        """Write a PeaksWorkspace to an ASCII file using this formatter.
+    class FileBuilder(object):
+        """Encapsulate information to build a single Jana file"""
+        def __init__(self, filepath, workspace, num_mod_vec, modulation_col_num=None):
+            self._filepath = filepath
+            self._workspace = workspace
+            self._num_mod_vec = num_mod_vec
+            self._modulation_col_num = modulation_col_num
+
+            self._headers = []
+            self._peaks = []
+            self._num_cols = None
+
+        def build_headers(self):
+            headers = self._headers
+            sample = self._workspace.sample()
+            lattice = sample.getOrientedLattice() if sample.hasOrientedLattice() else None
+
+            def append_mod_vector(index, col_num):
+                vec = lattice.getModVec(index)
+                x, y, z = vec.X(), vec.Y(), vec.Z()
+                if abs(x) > 0 or abs(y) > 0 or abs(z) > 0:
+                    headers.append("       {}{: >13.6f}{: >13.6f}{: >13.6f}\n".format(
+                        col_num, x, y, z))
+
+            # propagation vector information if required
+            if lattice is not None:
+                if self._num_mod_vec > 0:
+                    headers.append("# Structural propagation vectors used\n")
+                    if self._modulation_col_num is not None:
+                        headers.append("           1\n")
+                        append_mod_vector(self._modulation_col_num - 1, 1)
+                        headers.append("# Linear combination table used\n")
+                        headers.append("           2\n")
+                        headers.append("       1       1\n")
+                        headers.append("       2      -1\n")
+                    else:
+                        headers.append("           {}\n".format(self._num_mod_vec))
+                        for mod_vec_index in range(3):
+                            append_mod_vector(mod_vec_index, mod_vec_index + 1)
+                # lattice parameters
+                lattice_params = [
+                    lattice.a(),
+                    lattice.b(),
+                    lattice.c(),
+                    lattice.alpha(),
+                    lattice.beta(),
+                    lattice.gamma()
+                ]
+                lattice_params = "".join(["{: >10.4f}".format(value) for value in lattice_params])
+                headers.append("# Lattice parameters   {}\n".format(lattice_params))
+            headers.append("(3i5,2f12.2,i5,4f10.4)\n")
+            # column headers
+            column_names = ["h", "k", "l"]
+            if self._modulation_col_num is not None:
+                modulated_cols = ["m1"]
+            else:
+                modulated_cols = ["m{}".format(i + 1) for i in range(self._num_mod_vec)]
+            column_names.extend(modulated_cols)
+            column_names.extend(
+                ["Fsqr", "s(Fsqr)", "Cod", "Lambda", "Twotheta", "Transm.", "Tbar", "TDS"])
+
+            column_format = "#{:>4}{:>4}{:>4}"
+            column_format += "".join(["{:>4}" for _ in range(len(modulated_cols))])
+            column_format += "{:>12}{:>12}{:>5}{:>10}{:>10}{:>10}{:>10}{:>10}\n"
+            headers.append(column_format.format(*column_names))
+            self._num_cols = len(column_names)
+
+        def build_peaks_info(self):
+            for peak in self._workspace:
+                if self._num_mod_vec > 0:
+                    # if this is a main peak write it out. if not decide if it should be in this file
+                    hkl = peak.getIntHKL()
+                    mnp = peak.getIntMNP()
+                    if self._modulation_col_num is None:
+                        # write all modulation indices
+                        modulation_indices = [mnp[i] for i in range(self._num_mod_vec)]
+                    else:
+                        # is this a main peak or one with the modulation vector matching this file
+                        mnp_index = -1
+                        for i in range(3):
+                            if abs(mnp[i]) > 0.0:
+                                mnp_index = i
+                                break
+                        if mnp_index < 0:
+                            # all(mnp) == 0: main peak
+                            modulation_indices = [0]
+                        elif mnp_index == self._modulation_col_num - 1:
+                            # mnp index matches this file
+                            modulation_indices = [mnp[mnp_index]]
+                        else:
+                            # mnp doesn't match this file
+                            continue
+                else:
+                    # no modulated structure information
+                    hkl = peak.getHKL()
+                    modulation_indices = []
+                self._peaks.append(
+                    self.create_peak_line(hkl, modulation_indices,
+                                          peak.getIntensity(), peak.getSigmaIntensity(),
+                                          peak.getWavelength(),
+                                          get_two_theta(peak.getDSpacing(), peak.getWavelength())))
+
+        def create_peak_line(self, hkl, mnp, intensity, sig_int, wavelength, two_theta):
+            """
+            Write the raw peak data to a file.
+
+            :param f_handle: handle to the file to write to.
+            :param hkl: Integer HKL indices
+            :param mnp: List of mnp values
+            :param intensity: Intensity of the peak
+            :param sig_int: Signal value
+            :param wavelength: Wavelength in angstroms
+            :param two_theta: Two theta of detector
+            """
+            template = "{: >5.0f}{: >5.0f}{: >5.0f}{}{: >12.2f}{: >12.2f}{: >5.0f}{: >10.4f}{: >10.4f}{: >10.4f}{: >10.4f}{: >10.4f}\n"
+            mod_indices = "".join(["{: >5.0f}".format(value) for value in mnp])
+            return template.format(hkl[0], hkl[1], hkl[2], mod_indices, intensity, sig_int, 1, wavelength,
+                                   two_theta, 1.0, 0.0, 0.0)
+
+        def write(self):
+            with open(self._filepath, 'w') as handle:
+                handle.write("".join(self._headers))
+                handle.write("".join(self._peaks))
+
+    def __call__(self, file_name, workspace, split_files):
+        """Write a PeaksWorkspace or TableWorkspace with the appropriate columns
+        to an ASCII file using this formatter.
 
         :param file_name: the file name to output data to.
         :param workspace: the PeaksWorkspace to write to file.
+        :param split_files: if True the peaks associated with each
+                            modulation vector are saved to separate
+                            files. The suffix -mi where identifier=1,2...n
+                            is appended to each file
         """
-        with open(file_name, 'w') as f_handle:
-            self.write_header(f_handle, workspace)
-            self.write_peaks(f_handle, workspace)
+        builders = self._create_file_builders(file_name, workspace, split_files)
+        for builder in builders:
+            builder.build_headers()
+            builder.build_peaks_info()
+            builder.write()
 
-    def write_header(self, f_handle, workspace):
-        """Write the header of the Fullprof file format
-
-        :param f_handle: handle to the file to write to.
-        :param workspace: the PeaksWorkspace to save to file.
-        """
-        if isinstance(workspace, IPeaksWorkspace):
-            sample = workspace.sample()
-            lattice = sample.getOrientedLattice()
-            lattice_params = [lattice.a(), lattice.b(), lattice.c(), lattice.alpha(), lattice.beta(), lattice.gamma()]
-            lattice_params = "".join(["{: >10.4f}".format(value) for value in lattice_params])
-            f_handle.write("# Lattice parameters   {}\n".format(lattice_params))
-
-        f_handle.write("(3i5,2f12.2,i5,4f10.4)\n")
+    def _create_file_builders(self, file_name, workspace, split_files):
+        """Create a sequence of JanaFileBuilder to contain the information.
 
-    def write_peaks(self, f_handle, workspace):
-        """Write all the peaks in the workspace to file.
-
-        :param f_handle: handle to the file to write to.
-        :param workspace: the PeaksWorkspace to save to file.
+        :param file_name: Filename given by user
+        :param workspace: the PeaksWorkspace to write to file.
+        :param split_files: If true create a file for each modulation vector
+                            if there are more than 1
         """
-        column_names = ["h", "k", "l"]
-        column_names.extend(get_additional_index_names(workspace))
-        column_names.extend(["Fsqr", "s(Fsqr)", "Cod", "Lambda", "Twotheta", "Transm.", "Tbar", "TDS"])
-
-        column_format = "#{:>4}{:>4}{:>4}"
-        column_format += "".join(["{:>4}" for _ in range(num_additional_indicies(workspace))])
-        column_format += "{:>12}{:>12}{:>5}{:>10}{:>10}{:>10}{:>10}{:>10}\n"
-
-        f_handle.write(column_format.format(*column_names))
-        for row in workspace:
-            self.write_peak(f_handle, row, workspace)
+        num_mod_vec = num_modulation_vectors(workspace)
+        if split_files and num_mod_vec > 1:
+            name, ext = osp.splitext(file_name)
+            builders = [
+                JanaFormat.FileBuilder('{}-m{}{}'.format(name, col_num, ext), workspace,
+                                       num_mod_vec, col_num)
+                for col_num in range(1, num_mod_vec + 1)
+            ]
+        else:
+            builders = [JanaFormat.FileBuilder(file_name, workspace, num_mod_vec, None)]
 
-    def write_peak(self, f_handle, peak, workspace):
-        """Write a single Peak from the peaks workspace to file.
+        return builders
 
-        :param f_handle: handle to the file to write to.
-        :param peak: the peak object to write to the file.
-        """
-        ms = ["{: >5.0f}".format(peak[name]) for name in get_additional_index_names(workspace)]
-        f_handle.write("{h: >5.0f}{k: >5.0f}{l: >5.0f}".format(**peak))
-        f_handle.write("".join(ms))
-        f_handle.write("{Intens: >12.2f}".format(**peak))
-        f_handle.write("{SigInt: >12.2f}".format(**peak))
-        f_handle.write("{: >5.0f}".format(1))
-        f_handle.write("{Wavelength: >10.4f}".format(**peak))
-        f_handle.write("{: >10.4f}".format(self._get_two_theta(peak)))
-        f_handle.write("{: >10.4f}{: >10.4f}{: >10.4f}".format(1.0, 0.0, 0.0))
-        f_handle.write("\n")
-
-    def _get_two_theta(self, peak):
-        """Get the two theta value for this peak.
-
-        This is just Bragg's law relating wavelength to scattering angle.
-        :param peak: peak object to get the scattering angle for.
-        :returns: the scattering angle for the peak.
-        """
-        d = peak['DSpacing']
-        wavelength = peak['Wavelength']
-        theta = 2.*np.arcsin(0.5*(wavelength / d))
-        return np.rad2deg(theta)
 
 # ------------------------------------------------------------------------------------------------------
 
@@ -244,16 +357,33 @@ class SaveHKLFormat(object):
     The SaveHKL algorithm currently supports both the GSAS and SHELX formats. For
     more information see the SaveHKL algorithm documentation.
     """
-
-    def __call__(self, file_name, workspace):
+    def __call__(self, file_name, workspace, _):
         """Write a PeaksWorkspace to an ASCII file using this formatter.
 
         :param file_name: the file name to output data to.
         :param workspace: the PeaksWorkspace to write to file.
+        :param _: Ignored parameter for compatability with other savers
         """
         if has_modulated_indexing(workspace):
-            raise NotImplementedError("Cannot currently save modulated structures to GSAS or SHELX formats")
+            raise RuntimeError(
+                "Cannot currently save modulated structures to GSAS or SHELX formats")
+
+        from mantid.simpleapi import SaveHKL
         SaveHKL(Filename=file_name, InputWorkspace=workspace, OutputWorkspace=workspace.name())
 
 
+class ReflectionFormat(Enum):
+    Fullprof = 1
+    GSAS = 2
+    Jana = 3
+    SHELX = 4
+
+
+FORMAT_MAP = {
+    ReflectionFormat.Fullprof: FullprofFormat,
+    ReflectionFormat.GSAS: SaveHKLFormat,
+    ReflectionFormat.Jana: JanaFormat,
+    ReflectionFormat.SHELX: SaveHKLFormat
+}
+
 AlgorithmFactory.subscribe(SaveReflections)
diff --git a/Framework/PythonInterface/plugins/algorithms/fractional_indexing.py b/Framework/PythonInterface/plugins/algorithms/fractional_indexing.py
index eb51bf40fcddbcb0309a5ce34b5cbc642f842a6c..e9748e004e741e0ada5f2f92327ae1e5bc19c25d 100644
--- a/Framework/PythonInterface/plugins/algorithms/fractional_indexing.py
+++ b/Framework/PythonInterface/plugins/algorithms/fractional_indexing.py
@@ -209,7 +209,7 @@ def get_hkls(peaks_workspace):
     :param peaks_workspace: the peaks workspace to extract HKL values from
     :return: 2D numpy array of HKL values.
     """
-    return np.array([np.array([peak['h'], peak['k'], peak['l']])
+    return np.array([np.array([peak.getH(), peak.getK(), peak.getL()])
                      for peak in peaks_workspace])
 
 
diff --git a/Framework/PythonInterface/test/python/mantid/api/IPeaksWorkspaceTest.py b/Framework/PythonInterface/test/python/mantid/api/IPeaksWorkspaceTest.py
index 403c747fd5a0b4706c79ef0792d21c0da9010034..fa2174686c2e220bfb67efa2453029efdf78b060 100644
--- a/Framework/PythonInterface/test/python/mantid/api/IPeaksWorkspaceTest.py
+++ b/Framework/PythonInterface/test/python/mantid/api/IPeaksWorkspaceTest.py
@@ -64,11 +64,10 @@ class IPeaksWorkspaceTest(unittest.TestCase):
 
     def test_createPeakHKL(self):
         pws = WorkspaceCreationHelper.createPeaksWorkspace(0, True)
-        lattice = pws.mutableSample().getOrientedLattice()
 
         # Simple test that the creational method is exposed
         p = pws.createPeakHKL([1,1,1])
-        self.assertNotEqual(IPeak,  None)
+        self.assertFalse(IPeak is None)
 
     def test_peak_setQLabFrame(self):
         pws = WorkspaceCreationHelper.createPeaksWorkspace(1, True)
@@ -125,6 +124,20 @@ class IPeaksWorkspaceTest(unittest.TestCase):
         self.assertEqual(pws.cell("QLab", 0), V3D(1,1,1))
         self.assertEqual(pws.cell("QSample", 0), V3D(1,1,1))
 
+    def test_iteration_support(self):
+        pws = WorkspaceCreationHelper.createPeaksWorkspace(0, True)
+        hkls = ([1, 1, 1], [2, 1, 1], [1, 2, 1])
+        for hkl in hkls:
+            pws.addPeak(pws.createPeakHKL(hkl))
+
+        count = 0
+        for index, peak in enumerate(pws):
+            count += 1
+            self.assertTrue(isinstance(peak, IPeak))
+            self.assertAlmostEqual(V3D(*hkls[index]), peak.getHKL())
+
+        self.assertEquals(len(hkls), count)
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/Framework/PythonInterface/test/python/mantid/api/SampleTest.py b/Framework/PythonInterface/test/python/mantid/api/SampleTest.py
index 28d918d26083d3ea4059802cc57fa4db8e385199..c3f40c01be49a5960ec27c23e5ac106ec8cad7ca 100644
--- a/Framework/PythonInterface/test/python/mantid/api/SampleTest.py
+++ b/Framework/PythonInterface/test/python/mantid/api/SampleTest.py
@@ -7,15 +7,24 @@
 from __future__ import (absolute_import, division, print_function)
 
 import unittest
-from mantid.simpleapi import CreateWorkspace
-from mantid.simpleapi import SetSampleMaterial
-from mantid.geometry import CrystalStructure
-from mantid.geometry import CSGObject
+from mantid.simpleapi import CreateSampleWorkspace, CreatePeaksWorkspace, CreateWorkspace, SetSampleMaterial, SetUB
+from mantid.geometry import CrystalStructure, CSGObject, OrientedLattice
 from mantid.api import Sample
 from numpy import pi
 import copy
 
+
 class SampleTest(unittest.TestCase):
+    def test_lattice_accessors(self):
+        instrument_ws = CreateSampleWorkspace()
+        peaks = CreatePeaksWorkspace(instrument_ws, 0)
+        SetUB(peaks, 1, 1, 1, 90, 90, 90)
+        sample = peaks.sample()
+
+        self.assertTrue(sample.hasOrientedLattice())
+        self.assertTrue(isinstance(sample.getOrientedLattice(), OrientedLattice))
+        sample.clearOrientedLattice()
+        self.assertFalse(sample.hasOrientedLattice())
 
     def test_geometry_getters_and_setters(self):
         sample = Sample()
diff --git a/Framework/PythonInterface/test/python/mantid/geometry/UnitCellTest.py b/Framework/PythonInterface/test/python/mantid/geometry/UnitCellTest.py
index 405062d14d89771ee9b8cacd9d03e5c85bc5fde6..eb6083f46a170c3f8b60a7c406e2b4e75793adfc 100644
--- a/Framework/PythonInterface/test/python/mantid/geometry/UnitCellTest.py
+++ b/Framework/PythonInterface/test/python/mantid/geometry/UnitCellTest.py
@@ -59,6 +59,36 @@ class UnitCellTest(unittest.TestCase):
         self.assertEqual(unit.beta(), newUnit.beta())
         self.assertEqual(unit.gamma(), newUnit.gamma())
 
+    def test_getModVec_returns_v3d(self):
+        unit = UnitCell(3, 3, 3)
+
+        self.assertTrue(isinstance(unit.getModVec(0), V3D))
+        self.assertTrue(isinstance(unit.getModVec(1), V3D))
+        self.assertTrue(isinstance(unit.getModVec(2), V3D))
+
+    def test_setModVec_accepts_v3d(self):
+        unit = UnitCell(3, 3, 3)
+
+        vectors = (V3D(1,2,3), V3D(2, 1, 3), V3D(3, 2, 1))
+        unit.setModVec1(vectors[0])
+        unit.setModVec2(vectors[1])
+        unit.setModVec3(vectors[2])
+
+        for index, expected in enumerate(vectors):
+            self.assertEquals(vectors[index], unit.getModVec(index))
+
+    def test_getMaxOrder_returns_int(self):
+        unit = UnitCell(3, 3, 3)
+
+        self.assertTrue(isinstance(unit.getMaxOrder(), int))
+
+    def test_setMaxOrder_accepts_int(self):
+        unit = UnitCell(3, 3, 3)
+
+        unit.setMaxOrder(2)
+
+        self.assertEquals(2, unit.getMaxOrder())
+
     def _check_cell(self, cell):
         self.assertAlmostEqual(cell.a(),2.5,10)
         self.assertAlmostEqual(cell.b(),6,10)
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/IndexSatellitePeaksTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/IndexSatellitePeaksTest.py
index ce0a2e622cf04ce30796922eb92aa4ac7632430a..6093094c1dc25df3c1254166b1161949ac45a6fd 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/IndexSatellitePeaksTest.py
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/IndexSatellitePeaksTest.py
@@ -29,19 +29,27 @@ class IndexSatellitePeaksTest(unittest.TestCase):
         indexed_peaks = IndexSatellitePeaks(self._nuclear_peaks, self._peaks,
                                             Tolerance=0.1,
                                             ClusterThreshold=1.5, NumOfQs=-1)
-        index_values = np.array(indexed_peaks.column("m1"))
 
+        index_values = []
+        for peak in indexed_peaks:
+            mnp = peak.getIntMNP()
+            self.assertAlmostEqual(0.0, mnp[1])
+            self.assertAlmostEqual(0.0, mnp[2])
+            index_values.append(mnp[0])
         npt.assert_array_equal(index_values, expected_values)
-        self.assertRaises(RuntimeError, indexed_peaks.column, "m2")
 
     def test_exec_with_number_of_qs(self):
         expected_values = np.array([-1, -1, 1, 1])
         indexed_peaks = IndexSatellitePeaks(self._nuclear_peaks, self._peaks,
                                             Tolerance=0.1, NumOfQs=2)
-        index_values = np.array(indexed_peaks.column("m1"))
 
+        index_values = []
+        for peak in indexed_peaks:
+            mnp = peak.getIntMNP()
+            self.assertAlmostEqual(0.0, mnp[1])
+            self.assertAlmostEqual(0.0, mnp[2])
+            index_values.append(mnp[0])
         npt.assert_array_equal(index_values, expected_values)
-        self.assertRaises(RuntimeError, indexed_peaks.column, "m2")
 
 
 if __name__ == "__main__":
diff --git a/Framework/PythonInterface/test/python/plugins/algorithms/SaveReflectionsTest.py b/Framework/PythonInterface/test/python/plugins/algorithms/SaveReflectionsTest.py
index 74d18d93757a7a82025527b7babff30a4599be51..6f8900f74f67a4a1719abe99fa458fd94266a6c4 100644
--- a/Framework/PythonInterface/test/python/plugins/algorithms/SaveReflectionsTest.py
+++ b/Framework/PythonInterface/test/python/plugins/algorithms/SaveReflectionsTest.py
@@ -11,16 +11,24 @@ import shutil
 import numpy as np
 
 from mantid.api import FileFinder
-from mantid.simpleapi import SaveReflections, DeleteWorkspace, LoadEmptyInstrument, CreatePeaksWorkspace, SetUB, CreateEmptyTableWorkspace
+from mantid.kernel import V3D
+from mantid.simpleapi import (CloneWorkspace, CreatePeaksWorkspace, DeleteWorkspace,
+                              LoadEmptyInstrument, SetUB, SaveReflections)
 
 
 class SaveReflectionsTest(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls._test_dir = tempfile.mkdtemp()
+
+    @classmethod
+    def tearDownClass(cls):
+        shutil.rmtree(cls._test_dir)
+
     def setUp(self):
         self._workspace = self._create_peaks_workspace()
-        self._test_dir = tempfile.mkdtemp()
 
     def tearDown(self):
-        shutil.rmtree(self._test_dir)
         DeleteWorkspace(self._workspace)
 
     def _create_peaks_workspace(self):
@@ -54,44 +62,25 @@ class SaveReflectionsTest(unittest.TestCase):
         return ws
 
     def _create_modulated_peak_table(self):
-        hklm = np.ones((self._workspace.rowCount(), 2))
-        return self._create_indexed_workspace(self._workspace, 5, hklm)
-
-    def _create_indexed_workspace(self, fractional_peaks, ndim, hklm):
+        # alternating pattern of (1,0),(0,-1) and a (0,0)
+        nrows = self._workspace.rowCount()
+        m1m2 = np.zeros((nrows, 2))
+        m1m2[::2, 0] = 1
+        m1m2[1::2, 1] = -1
+        m1m2[int(nrows / 2)] = [0, 0]
+        return self._create_indexed_workspace(self._workspace, m1m2)
+
+    def _create_indexed_workspace(self, fractional_peaks, m1m2):
         # Create table with the number of columns we need
-        indexed = CreateEmptyTableWorkspace()
-        names = fractional_peaks.getColumnNames()
-        types = fractional_peaks.columnTypes()
-
-        # Insert the extra columns for the addtional indicies
-        for i in range(ndim - 3):
-            names.insert(5 + i, 'm{}'.format(i + 1))
-            types.insert(5 + i, 'double')
-
-        names = np.array(names)
-        types = np.array(types)
-
-        # Create columns in the table workspace
-        for name, column_type in zip(names, types):
-            indexed.addColumn(column_type, name)
+        modulated = CloneWorkspace(fractional_peaks)
+        lattice = modulated.sample().getOrientedLattice()
+        lattice.setModVec1(V3D(0.5, 0.0, 0.5))
+        lattice.setModVec2(V3D(0.333, 0.333, 0.))
+        for row, peak in enumerate(modulated):
+            row_indices = m1m2[row]
+            peak.setIntMNP(V3D(row_indices[0], row_indices[1], 0))
 
-        # Copy all columns from original workspace, ignoring HKLs
-        column_data = []
-        idx = np.arange(0, names.size)
-        hkl_mask = (idx < 5) | (idx > 4 + (ndim - 3))
-        for name in names[hkl_mask]:
-            column_data.append(fractional_peaks.column(name))
-
-        # Insert the addtional HKL columns into the data
-        for i, col in enumerate(hklm.T.tolist()):
-            column_data.insert(i + 2, col)
-
-        # Insert the columns into the table workspace
-        for i in range(fractional_peaks.rowCount()):
-            row = [column_data[j][i] for j in range(indexed.columnCount())]
-            indexed.addRow(row)
-
-        return indexed
+        return modulated
 
     def _get_reference_result(self, name):
         path = FileFinder.getFullPath(name)
@@ -136,7 +125,7 @@ class SaveReflectionsTest(unittest.TestCase):
         # Assert
         self._assert_file_content_equal(reference_result, file_name)
 
-    def test_save_jana_format_modulated(self):
+    def test_save_jana_format_modulated_single_file(self):
         # Arrange
         workspace = self._create_modulated_peak_table()
         reference_result = self._get_reference_result("jana_format_modulated.hkl")
@@ -149,6 +138,38 @@ class SaveReflectionsTest(unittest.TestCase):
         # Assert
         self._assert_file_content_equal(reference_result, file_name)
 
+    def test_save_jana_format_modulated_separate_files(self):
+        # Arrange
+        workspace = self._create_modulated_peak_table()
+        reference_results = [
+            self._get_reference_result("jana_format_modulated-m{}.hkl".format(i))
+            for i in range(1, 3)
+        ]
+        file_name = os.path.join(self._test_dir, "test_jana_modulated.hkl")
+        output_format = "Jana"
+        split_files = True
+
+        # Act
+        SaveReflections(InputWorkspace=workspace,
+                        Filename=file_name,
+                        Format=output_format,
+                        SplitFiles=split_files)
+
+        # Assert
+        self._assert_file_content_equal(reference_results[0],
+                                        os.path.join(self._test_dir, "test_jana_modulated-m1.hkl"))
+        self._assert_file_content_equal(reference_results[1],
+                                        os.path.join(self._test_dir, "test_jana_modulated-m2.hkl"))
+
+    def test_save_jana_with_no_lattice_information(self):
+        peaks = CloneWorkspace(self._workspace)
+        peaks.sample().clearOrientedLattice()
+        file_name = os.path.join(self._test_dir, "test_jana_no_lattice.hkl")
+
+        # Act
+        SaveReflections(InputWorkspace=peaks, Filename=file_name,
+                        Format="Jana", SplitFiles=False)
+
     def test_save_GSAS_format(self):
         # Arrange
         reference_result = self._get_reference_result("gsas_format.hkl")
@@ -159,7 +180,6 @@ class SaveReflectionsTest(unittest.TestCase):
         SaveReflections(InputWorkspace=self._workspace, Filename=file_name, Format=output_format)
 
         # Assert
-        #self._assert_file_content_equal(reference_result, file_name))
         self._assert_file_content_equal(reference_result, file_name)
 
     def test_save_GSAS_format_modulated(self):
@@ -215,13 +235,11 @@ class SaveReflectionsTest(unittest.TestCase):
     # Private api
     def _assert_file_content_equal(self, reference_result, file_name):
         with open(reference_result, 'r') as ref_file:
-            ref_lines = list(map(lambda x: x.strip(), ref_file.readlines()))
+            reference = ref_file.read()
         with open(file_name, 'r') as actual_file:
-            actual_lines = list(map(lambda x: x.strip(), actual_file.readlines()))
-
-        self.assertEqual(len(ref_lines), len(actual_lines))
-        for ref_line, actual_line in zip(ref_lines, actual_lines):
-            self.assertEqual(ref_line, actual_line)
+            actual = actual_file.read()
+            self.maxDiff = 10000
+        self.assertMultiLineEqual(reference, actual)
 
 
 if __name__ == '__main__':
diff --git a/Testing/Data/UnitTest/fullprof_format_modulated.hkl.md5 b/Testing/Data/UnitTest/fullprof_format_modulated.hkl.md5
index fedaea0aaa20334bac27e5fe844c918fcea228d2..5cd769a6db4d6aaf3eec2584c2fff89ba85d6d10 100644
--- a/Testing/Data/UnitTest/fullprof_format_modulated.hkl.md5
+++ b/Testing/Data/UnitTest/fullprof_format_modulated.hkl.md5
@@ -1 +1 @@
-5db1e4a64d6cbac5ae5cbf23dbc00897
+ac3eb47ad5f07fa38d7280c5d70f760f
diff --git a/Testing/Data/UnitTest/jana_format_modulated-m1.hkl.md5 b/Testing/Data/UnitTest/jana_format_modulated-m1.hkl.md5
new file mode 100644
index 0000000000000000000000000000000000000000..d750faf984746cfebd21fc2a1cb5097d6cdf30d2
--- /dev/null
+++ b/Testing/Data/UnitTest/jana_format_modulated-m1.hkl.md5
@@ -0,0 +1 @@
+cb6736cca510fe581bf928c7d5aaf5df
diff --git a/Testing/Data/UnitTest/jana_format_modulated-m2.hkl.md5 b/Testing/Data/UnitTest/jana_format_modulated-m2.hkl.md5
new file mode 100644
index 0000000000000000000000000000000000000000..9cdf3cb0f7dba190dcb3a4389da321657ba4cf00
--- /dev/null
+++ b/Testing/Data/UnitTest/jana_format_modulated-m2.hkl.md5
@@ -0,0 +1 @@
+612072cce6aa723ecc4fb6a58863e4a7
diff --git a/Testing/Data/UnitTest/jana_format_modulated.hkl.md5 b/Testing/Data/UnitTest/jana_format_modulated.hkl.md5
index 148d77dbf3face53dbe97926d2ec6c09a56131c0..28cfe89e0c3af5a467439b344ece64b2c416966f 100644
--- a/Testing/Data/UnitTest/jana_format_modulated.hkl.md5
+++ b/Testing/Data/UnitTest/jana_format_modulated.hkl.md5
@@ -1 +1 @@
-96c8aa0e0a91f2c537ec06f0170d7d25
+48a94ce1d7b1afe5ec8db4384e1582f1
diff --git a/Testing/SystemTests/lib/systemtests/systemtesting.py b/Testing/SystemTests/lib/systemtests/systemtesting.py
index 8f5a5fe5e2e115798df92de6524c69aa430a292a..d89f80919d5e12595ef4ca3d8e7b633acf157404 100644
--- a/Testing/SystemTests/lib/systemtests/systemtesting.py
+++ b/Testing/SystemTests/lib/systemtests/systemtesting.py
@@ -908,7 +908,7 @@ class TestManager(object):
         # A regex check is used to iterate back from the position of '.nxs' and
         # check that the current character is still part of a variable name. This
         # is needed to find the start of the string, hence the total filename.
-        check = re.compile("[A-Za-z0-9_-]")
+        check = re.compile(r"[A-Za-z0-9_-]")
         # In the case of looking for digits inside strings, the strings can start
         # with either " or '
         string_quotation_mark = ["'", '"']
diff --git a/Testing/SystemTests/tests/analysis/ISIS_WISHSingleCrystalReduction.py b/Testing/SystemTests/tests/analysis/ISIS_WISHSingleCrystalReduction.py
index 586c46972b213d8e8b7b548390f05371f17c5c86..c7a794c97003d051a1384765c7b8acae1f8c5cd5 100644
--- a/Testing/SystemTests/tests/analysis/ISIS_WISHSingleCrystalReduction.py
+++ b/Testing/SystemTests/tests/analysis/ISIS_WISHSingleCrystalReduction.py
@@ -61,11 +61,12 @@ class WISHSingleCrystalPeakPredictionTest(MantidSystemTest):
         # The peak at [-5 -1 -7] is known to fall between the gaps of WISH's tubes
         # Specifically check this one is predicted to exist because past bugs have
         # been found in the ray tracing.
-        Peak = namedtuple('Peak', ('DetID', 'BankName', 'h', 'k', 'l'))
-        expected = Peak(DetID=9202086, BankName='WISHpanel09', h=-5.0, k=-1.0, l=-7.0)
+        BasicPeak = namedtuple('Peak', ('DetID', 'BankName', 'h', 'k', 'l'))
+        expected = BasicPeak(DetID=9202086, BankName='WISHpanel09', h=-5.0, k=-1.0, l=-7.0)
         expected_peak_found = False
-        for row in self._filtered:
-            peak = Peak(DetID=row['DetID'], BankName=row['BankName'], h=row['h'], k=row['k'], l=row['l'])
+        for full_peak in self._filtered:
+            peak = BasicPeak(DetID=full_peak.getDetectorID(), BankName=full_peak.getBankName(),
+                             h=full_peak.getH(), k=full_peak.getK(), l=full_peak.getL())
             if peak == expected:
                 expected_peak_found = True
                 break
diff --git a/docs/source/algorithms/IndexSatellitePeaks-v1.rst b/docs/source/algorithms/IndexSatellitePeaks-v1.rst
index 4d7a644b49983da502e88a75b6cc22dd18f99b25..d6e20bb537dadd118080b2282c8da58ce68b93ca 100644
--- a/docs/source/algorithms/IndexSatellitePeaks-v1.rst
+++ b/docs/source/algorithms/IndexSatellitePeaks-v1.rst
@@ -73,10 +73,9 @@ Usage
     nuclear_peaks = Load('WISH_peak_hkl_small.nxs')
     satellite_peaks = Load("refine_satellites_fixed_q_test.nxs")
     indexed_peaks = IndexSatellitePeaks(nuclear_peaks, satellite_peaks, tolerance=0.1, NumOfQs=2)
-    index_values = np.array(indexed_peaks.column("m1"))
 
     for peak in indexed_peaks:
-        print("H: {h:>5} K: {k:>5} L: {l:>5} M1: {m1:>5}".format(**peak))
+        print("H: {:>5} K: {:>5} L: {:>5} M1: {:>5}".format(peak.getH(),peak.getK(),peak.getL(),peak.getIntMNP()[0]))
 
 Output:
 
diff --git a/docs/source/release/v4.3.0/diffraction.rst b/docs/source/release/v4.3.0/diffraction.rst
index ce5cd543f2728a2b3e2eb5060b166211d74a9182..59570624dbd1d81adf98b55cd1e6a079d0e39aa9 100644
--- a/docs/source/release/v4.3.0/diffraction.rst
+++ b/docs/source/release/v4.3.0/diffraction.rst
@@ -38,6 +38,8 @@ Single Crystal Diffraction
 - New algorithm :ref:`ConvertHFIRSCDtoMDE <algm-ConvertHFIRSCDtoMDE-v1>` for converting HFIR single crystal data (from WAND and DEMAND) into MDEventWorkspace in units Q_sample.
 - ``IndexPeaksWithsatellites`` has been deleted as it had been deprecated and superseded by :ref:`IndexPeaks <algm-IndexPeaks-v1>`.
 - The output peak workspace from :ref:`PredictFractionalPeaks<algm-PredictFractionalPeaks-v1>` now keeps the same lattice parameters as the input workspace. 
+- :ref:`SaveReflections <algm-SaveReflections>` now has the option to save peaks to separate files based on their associated modulation vectors
+  when using the Jana format.
 
 Imaging
 -------