diff --git a/docs/source/algorithms/PDFFourierTransform-v1.rst b/docs/source/algorithms/PDFFourierTransform-v1.rst
index 09aee1c8d5d6ef0927379848a518d521c0f6be05..e1ab0936bfd4d201b9efc06ce15b4e96e2e533a1 100644
--- a/docs/source/algorithms/PDFFourierTransform-v1.rst
+++ b/docs/source/algorithms/PDFFourierTransform-v1.rst
@@ -157,7 +157,6 @@ Usage
        print('! {0:4.2f} ! {1:5f} ! {2:f} ! {3:5f} !'.format(xx[i], yy[i], Rt.readX(0)[i], Rt.readY(0)[i]))
 
 
-
 .. testcleanup:: ExPDFFouurierTransform
 
    DeleteWorkspace(ws)
diff --git a/docs/source/algorithms/PaddingAndApodization-v1.rst b/docs/source/algorithms/PaddingAndApodization-v1.rst
index a366fffa392f442b37124b408548cb7dfd04ac56..a2a6dd46e16eef89178a4aae65e588424af214c2 100644
--- a/docs/source/algorithms/PaddingAndApodization-v1.rst
+++ b/docs/source/algorithms/PaddingAndApodization-v1.rst
@@ -37,7 +37,7 @@ Usage
    x = [1,2,3,4,5,6]
    input = CreateWorkspace(x,y)
    output=PaddingAndApodization(InputWorkspace=input,ApodizationFunction="Gaussian",DecayConstant=2.44,Padding=0,)
-   print  "output: ",['{0:.2f}'.format(value)  for value in output.readY(0)]
+   print("output:  {}".format(['{0:.2f}'.format(value) for value in output.readY(0)]))
    
 Output:
 
@@ -55,7 +55,7 @@ Output:
    x = [1,2,3,4,5,6]
    input = CreateWorkspace(x,y)
    output=PaddingAndApodization(InputWorkspace=input,Padding=2,)
-   print  "output: ",['{0:.2f}'.format(value)  for value in output.readY(0)]
+   print("output:  {}".format(['{0:.2f}'.format(value)  for value in output.readY(0)]))
    
 Output:
 
@@ -73,7 +73,7 @@ Output:
    x = [1,2,3,4,5,6]
    input = CreateWorkspace(x,y)
    output=PaddingAndApodization(InputWorkspace=input,ApodizationFunction="Gaussian",DecayConstant=2.44,Padding=2,)
-   print  "output: ",['{0:.2f}'.format(value)  for value in output.readY(0)]
+   print("output:  {}".format(['{0:.2f}'.format(value)  for value in output.readY(0)]))
    
 Output:
 
diff --git a/docs/source/algorithms/Pause-v1.rst b/docs/source/algorithms/Pause-v1.rst
index 6c0b7a3ae44e17b6ce6975ea683848c3b7177cbb..3608265c5ab15ceb0f2ef93cbdeb3eb352424094 100644
--- a/docs/source/algorithms/Pause-v1.rst
+++ b/docs/source/algorithms/Pause-v1.rst
@@ -28,7 +28,7 @@ Usage
    start_time = time.clock()
    Pause(0.05)
    end_time = time.clock()
-   print ("The algorithm paused for %.2f seconds." % (end_time-start_time))
+   print("The algorithm paused for {:.2f} seconds.".format(end_time-start_time))
 	
 Output:
 
diff --git a/docs/source/algorithms/PawleyFit-v1.rst b/docs/source/algorithms/PawleyFit-v1.rst
index d41328f21a72f3d5e733b00dfd75750ad1950148..2d595ce71e36b1d5d8f4801b24d7cd0e5b2248c7 100644
--- a/docs/source/algorithms/PawleyFit-v1.rst
+++ b/docs/source/algorithms/PawleyFit-v1.rst
@@ -36,7 +36,7 @@ For the usage example there is a calculated, theoretical diffraction pattern (in
                                   Atoms='Si 0 0 0 1.0 0.05',
                                   a=5.43, LatticeSpacingMin=0.7)
 
-    print "Silicon has", Si.rowCount(), "unique reflections with d > 0.7."
+    print("Silicon has {} unique reflections with d > 0.7.".format(Si.rowCount()))
 
     # Find peaks in the spectrum
     si_peaks = PoldiPeakSearch(si_spectrum)
@@ -47,7 +47,7 @@ For the usage example there is a calculated, theoretical diffraction pattern (in
     si_peaks_indexed = AnalysisDataService.retrieve('si_peaks_indexed_Si')
 
     # 3 peaks have two possibilities for indexing, because their d-values are identical
-    print "The number of peaks that were indexed:", si_peaks_indexed.rowCount()
+    print("The number of peaks that were indexed: {}".format(si_peaks_indexed.rowCount()))
 
     # Run the actual fit with lattice parameters that are slightly off
     si_fitted, si_cell, si_params, chi_square = PawleyFit(si_spectrum,
@@ -61,10 +61,10 @@ For the usage example there is a calculated, theoretical diffraction pattern (in
     a_err = np.round(si_cell.cell(0, 2), 6)
     a_diff = np.round(np.fabs(a - 5.4311946), 6)
 
-    print "The lattice parameter was refined to a =", a, "+/-", a_err
-    print "The deviation from the actual parameter (a=5.4311946) is:", a_diff
-    print "This difference corresponds to", np.round(a_diff / a_err, 2), "standard deviations."
-    print "The reduced chi square of the fit is:", np.round(chi_square, 3)
+    print("The lattice parameter was refined to a = {} +/- {}".format(a, a_err))
+    print("The deviation from the actual parameter (a=5.4311946) is: {}".format(a_diff))
+    print("This difference corresponds to {:.2f} standard deviations.".format(np.round(a_diff / a_err, 2)))
+    print("The reduced chi square of the fit is: {:.2f}".format(np.round(chi_square, 3)))
 
 Running this script will generate a bit of output about the results of the different steps. At the end the lattice parameter differs less than one standard deviation from the actual value.
 
diff --git a/docs/source/algorithms/PeaksInRegion-v1.rst b/docs/source/algorithms/PeaksInRegion-v1.rst
index 9481f3552528e837c8a8b7c2d9cf02a777ab0ef3..d31f554b61c8fcc59fd6d0d4341b0551a2fd25ed 100644
--- a/docs/source/algorithms/PeaksInRegion-v1.rst
+++ b/docs/source/algorithms/PeaksInRegion-v1.rst
@@ -26,15 +26,15 @@ Usage
    
    # Find peaks in region when the Peak sits in the centre of a box
    in_box_table = PeaksInRegion(peaks, CoordinateFrame='HKL', PeakRadius=0.1, Extents=[-1.0,1.0,-1.0,1.0,-1.0,1.0], CheckPeakExtents=True)
-   print in_box_table.row(0)
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**in_box_table.row(0)))
    
    # Find peaks in region when the peak is just outside the box (by radius)
    just_outside_box_table = PeaksInRegion(peaks, CoordinateFrame='HKL', PeakRadius=0.999, Extents=[1.0,2.0,-1.0,1.0,-1.0,1.0], CheckPeakExtents=True)
-   print just_outside_box_table.row(0)
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**just_outside_box_table.row(0)))
    
    # Find peaks in region when the peak is just inside the box (by radius)
-   just_intesecting_box_table = PeaksInRegion(peaks, CoordinateFrame='HKL', PeakRadius=1.00, Extents=[1.0,2.0,-1.0,1.0,-1.0,1.0], CheckPeakExtents=True)
-   print just_intesecting_box_table.row(0)
+   just_intersecting_box_table = PeaksInRegion(peaks, CoordinateFrame='HKL', PeakRadius=1.00, Extents=[1.0,2.0,-1.0,1.0,-1.0,1.0], CheckPeakExtents=True)
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**just_intersecting_box_table.row(0)))
    
 Output:
 
diff --git a/docs/source/algorithms/PeaksOnSurface-v1.rst b/docs/source/algorithms/PeaksOnSurface-v1.rst
index 82685a55679dfb6b1015c709e6d5ccaf01986851..2120e70924cc30974d541d4fd1654970672f421a 100644
--- a/docs/source/algorithms/PeaksOnSurface-v1.rst
+++ b/docs/source/algorithms/PeaksOnSurface-v1.rst
@@ -28,17 +28,23 @@ Usage
    
    # Peak is on the plane
    out_of_plane_offset = 0
-   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=1.0, CoordinateFrame='HKL', Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset], Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
-   print tbl.row(0)
+   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=1.0, CoordinateFrame='HKL',
+                        Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset],
+                        Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**tbl.row(0)))
    
    # Peak is off the plane, and does not intesect it
    out_of_plane_offset = 1.000
-   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=0.999,  CoordinateFrame='HKL', Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset], Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
-   print tbl.row(0)
+   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=0.999,  CoordinateFrame='HKL',
+                        Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset],
+                        Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**tbl.row(0)))
    
    # Peak is off the plane, but does intesect it when radius is made larger
-   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=1.000,  CoordinateFrame='HKL', Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset], Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
-   print tbl.row(0)
+   tbl = PeaksOnSurface(InputWorkspace=peaks, PeakRadius=1.000,  CoordinateFrame='HKL',
+                        Vertex1=[1.0, -1.0, out_of_plane_offset], Vertex2=[-1.0,-1.0,out_of_plane_offset],
+                        Vertex3=[-1.0, 1.0,out_of_plane_offset], Vertex4=[1.0, 1.0,out_of_plane_offset])
+   print("{{'Distance': {Distance}, 'PeakIndex': {PeakIndex}, 'Intersecting': {Intersecting}}}".format(**tbl.row(0)))
    
 Output:
 
diff --git a/docs/source/algorithms/PerformIndexOperations-v1.rst b/docs/source/algorithms/PerformIndexOperations-v1.rst
index 8cf4e0664f65bae59aa2549476257f0b9f280c00..8884f4c92dd4a072e7596e16fdf01b46fe6ae982 100644
--- a/docs/source/algorithms/PerformIndexOperations-v1.rst
+++ b/docs/source/algorithms/PerformIndexOperations-v1.rst
@@ -40,9 +40,9 @@ Usage
    ws2 = PerformIndexOperations( ws, "0+1,3,4")
 
    #print result
-   print ws2.readY(0)
-   print ws2.readY(1)
-   print ws2.readY(2)
+   print(ws2.readY(0))
+   print(ws2.readY(1))
+   print(ws2.readY(2))
    
 Output:
 
diff --git a/docs/source/algorithms/PhaseQuad-v1.rst b/docs/source/algorithms/PhaseQuad-v1.rst
index 0cefabfd3a45a5d6ef630d5a1dc8dff5e97b3591..30f6d6da89ecb6bfaaaed11826ccfe8ae8e9b2c5 100644
--- a/docs/source/algorithms/PhaseQuad-v1.rst
+++ b/docs/source/algorithms/PhaseQuad-v1.rst
@@ -54,7 +54,7 @@ Usage
        phi = 2*pi*i/32.
        tab.addRow([1, 0.2, phi])
    ows = PhaseQuad(InputWorkspace='MUSR00022725', PhaseTable='tab')
-   print "Output workspace has", ows.getNumberHistograms(), "histograms"
+   print("Output workspace has {} histograms".format(ows.getNumberHistograms()))
 
 Output:
 
diff --git a/docs/source/algorithms/PlotAsymmetryByLogValue-v1.rst b/docs/source/algorithms/PlotAsymmetryByLogValue-v1.rst
index b249120dcb5d00bc1f6a7bbe56c38e2342a087d7..bcdd84cedb0f572a9197683c2a27158e2aefd9c9 100644
--- a/docs/source/algorithms/PlotAsymmetryByLogValue-v1.rst
+++ b/docs/source/algorithms/PlotAsymmetryByLogValue-v1.rst
@@ -55,8 +55,8 @@ Usage
                                 TimeMin = 0.55,
                                 TimeMax = 12.0);
 
-   print "Y values (asymmetry):", ws.readY(0)
-   print "X values (sample magn. field):", ws.readX(0)
+   print("Y values (asymmetry): {}".format(ws.readY(0)))
+   print("X values (sample magn. field): {}".format(ws.readX(0)))
 
 Output:
 
@@ -77,11 +77,11 @@ Output:
                                 Red = 1,
                                 Green = 2);
 
-   print "Y values (difference):", ws.readY(0)
-   print "Y values (red):", ws.readY(1)
-   print "Y values (green):", ws.readY(2)
-   print "Y values (sum):", ws.readY(3)
-   print "X values (sample magn. field):", ws.readX(0)
+   print("Y values (difference): {}".format(ws.readY(0)))
+   print("Y values (red): {}".format(ws.readY(1)))
+   print("Y values (green): {}".format(ws.readY(2)))
+   print("Y values (sum): {}".format(ws.readY(3)))
+   print("X values (sample magn. field): {}".format(ws.readX(0)))
 
 Output:
 
@@ -98,8 +98,8 @@ Output:
 .. testcode:: ExCustomGrouping
 
    # Skip spectra 35
-   fwd_spectra = range(33,35) + range(36,65)
-
+   fwd_spectra = [x for x in range(33, 65) if x != 35]
+   
    # Skip spectra 1 and 2
    bwd_spectra = range(3, 33)
 
@@ -111,10 +111,10 @@ Output:
                                 ForwardSpectra = fwd_spectra,
                                 BackwardSpectra = bwd_spectra)
 
-   print "No of forward spectra used:", len(fwd_spectra)
-   print "No of backward spectra used:", len(bwd_spectra)
-   print "Y values (asymmetry):", ws.readY(0)
-   print "X values (sample magn. field):", ws.readX(0)
+   print("No of forward spectra used: {}".format(len(fwd_spectra)))
+   print("No of backward spectra used: {}".format(len(bwd_spectra)))
+   print("Y values (asymmetry): {}".format(ws.readY(0)))
+   print("X values (sample magn. field): {}".format(ws.readX(0)))
 
 Output:
 
@@ -136,8 +136,8 @@ Output:
                                 TimeMax = 12.0,
                                 DeadTimeCorrType = 'FromRunData');
 
-   print "Y values (asymmetry):", ws.readY(0)
-   print "X values (sample magn. field):", ws.readX(0)
+   print("Y values (asymmetry): {}".format(ws.readY(0)))
+   print("X values (sample magn. field): {}".format(ws.readX(0)))
 
 Output:
 
@@ -154,8 +154,8 @@ Output:
                                 LastRun="MUSR00015191",
                                 LogValue="sample_temp",
                                 Function="Mean")
-   print "Y values (asymmetry):", ws.readY(0)
-   print "X values (sample magn. field):", ws.readX(0)
+   print("Y values (asymmetry): {}".format(ws.readY(0)))
+   print("X values (sample magn. field): {}".format(ws.readX(0)))
 
 Output:
 
diff --git a/docs/source/algorithms/PlotPeakByLogValue-v1.rst b/docs/source/algorithms/PlotPeakByLogValue-v1.rst
index 5a8111b601f626daf11ba4fbbbb732a5a4c925d3..2abc311e9edcd62cf5adcbc398e0b0b3d52b5215 100644
--- a/docs/source/algorithms/PlotPeakByLogValue-v1.rst
+++ b/docs/source/algorithms/PlotPeakByLogValue-v1.rst
@@ -118,7 +118,7 @@ Usage
     ref = np.empty(len(peak_centres))
     ref.fill(10098.6)
 
-    print np.allclose(ref, peak_centres, 1e-3)
+    print(np.allclose(ref, peak_centres, 1e-3))
 
 Output:
 
diff --git a/docs/source/algorithms/Plus-v1.rst b/docs/source/algorithms/Plus-v1.rst
index 65d9c7c450e983abf1a23dcc646efdcfeb3042b9..ae8542c8964a4bea6378a872c4a0ceabad1dc68b 100644
--- a/docs/source/algorithms/Plus-v1.rst
+++ b/docs/source/algorithms/Plus-v1.rst
@@ -53,9 +53,9 @@ Usage
     # perform the algorithm
     ws = Plus(ws1, ws2)
 
-    print "The X values are: " + str(ws.readX(0))
-    print "The Y values are: " + str(ws.readY(0))
-    print "The updated Error values are: " + str(ws.readE(0))
+    print("The X values are: " + str(ws.readX(0)))
+    print("The Y values are: " + str(ws.readY(0)))
+    print("The updated Error values are: " + str(ws.readE(0)))
 
 Output:
 
@@ -83,9 +83,9 @@ Output:
     # perform the algorithm
     ws = ws1 + ws2
 
-    print "The X values are: " + str(ws.readX(0))
-    print "The Y values are: " + str(ws.readY(0))
-    print "The updated Error values are: " + str(ws.readE(0))
+    print("The X values are: " + str(ws.readX(0)))
+    print("The Y values are: " + str(ws.readY(0)))
+    print("The updated Error values are: " + str(ws.readE(0)))
 
 Output:
 
@@ -113,9 +113,9 @@ Output:
     # perform the algorithm
     ws += ws1
 
-    print "The X values are: " + str(ws.readX(0))
-    print "The Y values are: " + str(ws.readY(0))
-    print "The updated Error values are: " + str(ws.readE(0))
+    print("The X values are: " + str(ws.readX(0)))
+    print("The Y values are: " + str(ws.readY(0)))
+    print("The updated Error values are: " + str(ws.readE(0)))
 
 Output:
 
@@ -139,9 +139,9 @@ Output:
     # perform the algorithm
     ws = ws1 + 2.5
 
-    print "The X values are: " + str(ws.readX(0))
-    print "The Y values are: " + str(ws.readY(0))
-    print "The updated Error values are: " + str(ws.readE(0))
+    print("The X values are: " + str(ws.readX(0)))
+    print("The Y values are: " + str(ws.readY(0)))
+    print("The updated Error values are: " + str(ws.readE(0)))
 
 Output:
 
diff --git a/docs/source/algorithms/PoissonErrors-v1.rst b/docs/source/algorithms/PoissonErrors-v1.rst
index a7fca58bc4f3a66897b2d33f86509205dba8831e..ff71601af0216fb910ace03932b3b515157f2bda 100644
--- a/docs/source/algorithms/PoissonErrors-v1.rst
+++ b/docs/source/algorithms/PoissonErrors-v1.rst
@@ -32,9 +32,9 @@ Usage
     ws = PoissonErrors(ws1, ws2)
 
     #X-values aren't touched at all, Y-Values are used but not altered, E-Values are 0 if rhsY is 0 or (rhsE/rhsY)*lshY if they are non-zero
-    print "The X values are: " + str(ws.readX(0))
-    print "The Y values are: " + str(ws.readY(0))
-    print "The updated Error values are: " + str(ws.readE(0))
+    print("The X values are: " + str(ws.readX(0)))
+    print("The Y values are: " + str(ws.readY(0)))
+    print("The updated Error values are: " + str(ws.readE(0)))
 
 Output:
 
diff --git a/docs/source/algorithms/PoldiAnalyseResiduals-v1.rst b/docs/source/algorithms/PoldiAnalyseResiduals-v1.rst
index 4869b7df1f8891715a9c567c418b7327bf9d13d7..db4d53d023643f1509e39640ba992acc1263319e 100644
--- a/docs/source/algorithms/PoldiAnalyseResiduals-v1.rst
+++ b/docs/source/algorithms/PoldiAnalyseResiduals-v1.rst
@@ -51,7 +51,7 @@ The following example shows how to calculate the residuals following a fit perfo
 
     residual_data = residuals_Si.readY(0)
 
-    print "Residuals are in the range: [", round(min(residual_data), 2), ", ", round(max(residual_data), 2), "]"
+    print("Residuals are in the range: [ {:.2f} ,  {:.2f} ]".format(round(min(residual_data), 2), round(max(residual_data), 2)))
 
 The output contains the range in which residuals are found:
     
diff --git a/docs/source/algorithms/PoldiAutoCorrelation-v5.rst b/docs/source/algorithms/PoldiAutoCorrelation-v5.rst
index 40a0dd6a00a580e14fd07d1e453472d5697fa669..06d1a4746d6ce6343dde224b860e3c26879be2d2 100644
--- a/docs/source/algorithms/PoldiAutoCorrelation-v5.rst
+++ b/docs/source/algorithms/PoldiAutoCorrelation-v5.rst
@@ -44,7 +44,7 @@ PoldiAutoCorrelation operates on a MatrixWorkspace with a valid POLDI instrument
     raw_6904 = LoadSINQFile(Filename = "poldi2013n006904.hdf", Instrument = "POLDI")
     
     # Print the number of spectra in the workspace. It should be 400, one for each detector wire
-    print "The workspace contains", raw_6904.getNumberHistograms(), "spectra."
+    print("The workspace contains {} spectra.".format(raw_6904.getNumberHistograms()))
     
     # For most calculations, an instrument definition is needed, so it's loaded as well
     LoadInstrument(raw_6904, InstrumentName = "POLDI", RewriteSpectraMap=True)
@@ -54,7 +54,7 @@ PoldiAutoCorrelation operates on a MatrixWorkspace with a valid POLDI instrument
     
     # The first spectrum contains the correlation data. In this case there should be 5531 bins.
     # On other data or different wavelength limits, this number will be different.
-    print "The correlation spectrum has", len(correlated_6904.readY(0)), "data points."
+    print("The correlation spectrum has {} data points.".format(len(correlated_6904.readY(0))))
 
 Output:    
     
diff --git a/docs/source/algorithms/PoldiCreatePeaksFromCell-v1.rst b/docs/source/algorithms/PoldiCreatePeaksFromCell-v1.rst
index 2f67cb786d7470bb2c434214317ea6e23c11bdbb..474852eb1f01b2eeb3ef66087403bf74180b5fba 100644
--- a/docs/source/algorithms/PoldiCreatePeaksFromCell-v1.rst
+++ b/docs/source/algorithms/PoldiCreatePeaksFromCell-v1.rst
@@ -40,7 +40,7 @@ The following usage example illustrates how the algorithm can be used to generat
                         a=4.126,
                         LatticeSpacingMin=0.55, LatticeSpacingMax=4.0)
 
-    print "CsCl has", csClReflections.rowCount(), "unique reflections in the range between 0.55 and 4.0 Angstrom."
+    print("CsCl has {} unique reflections in the range between 0.55 and 4.0 Angstrom.".format(csClReflections.rowCount()))
 
 Output: