diff --git a/docs/source/algorithms/NMoldyn4Interpolation-v1.rst b/docs/source/algorithms/NMoldyn4Interpolation-v1.rst
index b14f416241a9091966cd75c8fde0f9f883fec868..6209957c9c0dd207260d1e56a70b27cb52b3ebee 100644
--- a/docs/source/algorithms/NMoldyn4Interpolation-v1.rst
+++ b/docs/source/algorithms/NMoldyn4Interpolation-v1.rst
@@ -38,9 +38,9 @@ set**
     osiris = Rebin(osiris, [-0.6, 0.02, 0.6])
     #interpolate the two workspaces
     interpolated_ws = NMoldyn4Interpolation(sim_ws, osiris)
-    print 'No. of Q-values in simulation = ' + str(sim_ws.getNumberHistograms())
-    print 'No. of Q-values in reference = ' + str(osiris.getNumberHistograms())
-    print 'No. of Q-values in interpolated set = '+ str(interpolated_ws.getNumberHistograms())
+    print('No. of Q-values in simulation = {}'.format(sim_ws.getNumberHistograms()))
+    print('No. of Q-values in reference = {}'.format(osiris.getNumberHistograms()))
+    print('No. of Q-values in interpolated set = {}'.format(interpolated_ws.getNumberHistograms()))
 
 Output:
 
diff --git a/docs/source/algorithms/NRCalculateSlitResolution-v1.rst b/docs/source/algorithms/NRCalculateSlitResolution-v1.rst
index 5885a547200c0426e61ddc06f37a7ed4fd3bfcdb..bdd830f369c35ee35ef6af028d065c6bb5ed052b 100644
--- a/docs/source/algorithms/NRCalculateSlitResolution-v1.rst
+++ b/docs/source/algorithms/NRCalculateSlitResolution-v1.rst
@@ -50,7 +50,7 @@ Usage
 
   ws = Load('INTER00013460')
   res = NRCalculateSlitResolution(Workspace = ws, TwoTheta = 0.7 * 2)
-  print("Resolution: %.4f" % res)
+  print("Resolution: {:.4f}".format(res))
 
 .. testoutput::
 
diff --git a/docs/source/algorithms/NormaliseByCurrent-v1.rst b/docs/source/algorithms/NormaliseByCurrent-v1.rst
index 45c3e3fb2d065804f9b00307da242bb404686fd9..ad9fbc58e87251e0fa5ff0a85b098a997b5d7833 100644
--- a/docs/source/algorithms/NormaliseByCurrent-v1.rst
+++ b/docs/source/algorithms/NormaliseByCurrent-v1.rst
@@ -55,16 +55,16 @@ Usage
    log_p = run1.getLogData('gd_prtn_chrg')
 
    # Print the log value
-   print "Good Proton Charge =",log_p.value
+   print("Good Proton Charge = {}".format(log_p.value))
 
    #Run the Algorithm
    wsN = NormaliseByCurrent(ws)
    norm_factor = wsN.getRun().getLogData('NormalizationFactor').value
 
    #Print results
-   print "Before normalisation", ws.readY(0);
-   print "After normalisation ", wsN.readY(0);
-   print "Normalisation factor", norm_factor;
+   print("Before normalisation {}".format(ws.readY(0)))
+   print("After normalisation  {}".format(wsN.readY(0)))
+   print("Normalisation factor {}".format(norm_factor))
 
 
 Output:
diff --git a/docs/source/algorithms/NormaliseByDetector-v1.rst b/docs/source/algorithms/NormaliseByDetector-v1.rst
index 92c2c43bde8b9ebd10bb4070e6ff9ec028e114f4..a9d4792005530a0e4bbcdc2cfe64d265fe28fb35 100644
--- a/docs/source/algorithms/NormaliseByDetector-v1.rst
+++ b/docs/source/algorithms/NormaliseByDetector-v1.rst
@@ -200,10 +200,10 @@ Usage
   #Now we are ready to run the correction
   wsCorrected = NormaliseByDetector(ws)
 
-  print ("The correction will divide the data by an increasing linear function.")
-  print ("f(x) = 2x + 1")
+  print("The correction will divide the data by an increasing linear function.")
+  print("f(x) = 2x + 1")
   for i in range(0,wsCorrected.blocksize(),10):
-    print ("The correct value in bin %i is %.2f compared to %.2f" % (i,wsCorrected.readY(0)[i],ws.readY(0)[i]))
+    print("The correct value in bin {} is {:.2f} compared to {:.2f}".format(i,wsCorrected.readY(0)[i],ws.readY(0)[i]))
 
   #clean up the file
   if os.path.exists(param_file_path):
diff --git a/docs/source/algorithms/NormaliseByPeakArea-v1.rst b/docs/source/algorithms/NormaliseByPeakArea-v1.rst
index 5fce5080943e3a116e9e6bd7680a21d69ae7326e..fb3e25279aeddd584f1a0cfb7ae3c89c4766b9f9 100644
--- a/docs/source/algorithms/NormaliseByPeakArea-v1.rst
+++ b/docs/source/algorithms/NormaliseByPeakArea-v1.rst
@@ -55,10 +55,10 @@ Usage
     normalised, yspace, fitted, symmetrised = \
       NormaliseByPeakArea(InputWorkspace=tof_ws, Mass=1.0079,Sum=False)
 
-    print "Number of normalised spectra is: %d" % normalised.getNumberHistograms()
-    print "Number of Y-space spectra is: %d" % yspace.getNumberHistograms()
-    print "Number of fitted spectra is: %d" % fitted.getNumberHistograms()
-    print "Number of symmetrised spectra is: %d" % symmetrised.getNumberHistograms()
+    print("Number of normalised spectra is: {}".format(normalised.getNumberHistograms()))
+    print("Number of Y-space spectra is: {}".format(yspace.getNumberHistograms()))
+    print("Number of fitted spectra is: {}".format(fitted.getNumberHistograms()))
+    print("Number of symmetrised spectra is: {}".format(symmetrised.getNumberHistograms()))
 
 .. testoutput:: NormaliseNoSumOutput
 
@@ -87,10 +87,10 @@ Usage
     normalised, yspace, fitted, symmetrised = \
       NormaliseByPeakArea(InputWorkspace=tof_ws, Mass=1.0079,Sum=True)
 
-    print "Number of normalised spectra is: %d" % normalised.getNumberHistograms()
-    print "Number of Y-space spectra is: %d" % yspace.getNumberHistograms()
-    print "Number of fitted spectra is: %d" % fitted.getNumberHistograms()
-    print "Number of symmetrised spectra is: %d" % symmetrised.getNumberHistograms()
+    print("Number of normalised spectra is: {}".format(normalised.getNumberHistograms()))
+    print("Number of Y-space spectra is: {}".format(yspace.getNumberHistograms()))
+    print("Number of fitted spectra is: {}".format(fitted.getNumberHistograms()))
+    print("Number of symmetrised spectra is: {}".format(symmetrised.getNumberHistograms()))
 
 .. testoutput:: NormaliseWithSummedOutput
 
diff --git a/docs/source/algorithms/NormaliseByThickness-v1.rst b/docs/source/algorithms/NormaliseByThickness-v1.rst
index 8843002d2daddccb7e1dd17274339bb04eb0660d..437f5584ca7c451946eb66d8420a50b666d862cb 100644
--- a/docs/source/algorithms/NormaliseByThickness-v1.rst
+++ b/docs/source/algorithms/NormaliseByThickness-v1.rst
@@ -24,11 +24,11 @@ Usage
     norm=NormaliseByThickness(raw,SampleThickness=10)
 
     #do a quick check
-    print norm[1]
-    print "Min(raw)=",raw.dataY(0).min()
-    print "Min(norm)=",norm[0].dataY(0).min()
-    print "Max(raw)=",raw.dataY(0).max()
-    print "Max(norm)=",norm[0].dataY(0).max()   
+    print(norm[1])
+    print("Min(raw)= {}".format(raw.dataY(0).min()))
+    print("Min(norm)= {}".format(norm[0].dataY(0).min()))
+    print("Max(raw)= {}".format(raw.dataY(0).max()))
+    print("Max(norm)= {}".format(norm[0].dataY(0).max()))
     
     
 .. testcleanup:: NormaliseByThicness
diff --git a/docs/source/algorithms/NormaliseSpectra-v1.rst b/docs/source/algorithms/NormaliseSpectra-v1.rst
index ac5c15344a6162a0e70809562beaa0b503ae8d1c..d7e428c066435768ca84bae1a9401c48a9d23468 100644
--- a/docs/source/algorithms/NormaliseSpectra-v1.rst
+++ b/docs/source/algorithms/NormaliseSpectra-v1.rst
@@ -32,7 +32,7 @@ Usage
   out_ws = NormaliseSpectra(InputWorkspace=ws)
   
   # Print resulting y values
-  print out_ws.readY(0)
+  print(out_ws.readY(0))
 
 Output:  
   
diff --git a/docs/source/algorithms/NormaliseToMonitor-v1.rst b/docs/source/algorithms/NormaliseToMonitor-v1.rst
index 61a340bcd1e99d2bce15d43a8b31a6c137211b24..1580d16a7d168417a5a084a853e279b5e64c475c 100644
--- a/docs/source/algorithms/NormaliseToMonitor-v1.rst
+++ b/docs/source/algorithms/NormaliseToMonitor-v1.rst
@@ -81,13 +81,13 @@ Usage
 
    wsN = NormaliseToMonitor( ws, MonitorID=1 )
 
-   print "Without normalisation"
-   print "Monitor ID=1 %.3f, %.3f" % ( ws.readY(0)[0], ws.readY(0)[1] )
-   print "Selected data %.6f, %.6f" % ( ws.readY(6)[0], ws.readY(3)[1] )
+   print("Without normalisation")
+   print("Monitor ID=1 {:.3f}, {:.3f}".format(ws.readY(0)[0], ws.readY(0)[1]))
+   print("Selected data {:.6f}, {:.6f}".format(ws.readY(6)[0], ws.readY(3)[1]))
 
-   print "With Normalisation"
-   print "Monitor ID=1 %.3f, %.3f" % ( wsN.readY(0)[0], wsN.readY(0)[1] )
-   print "Selected data %.6f, %.6f" % ( wsN.readY(6)[0], wsN.readY(3)[1] )
+   print("With Normalisation")
+   print("Monitor ID=1 {:.3f}, {:.3f}".format(wsN.readY(0)[0], wsN.readY(0)[1]))
+   print("Selected data {:.6f}, {:.6f}".format(wsN.readY(6)[0], wsN.readY(3)[1]))
 
 Output:
 
diff --git a/docs/source/algorithms/NormaliseToUnity-v1.rst b/docs/source/algorithms/NormaliseToUnity-v1.rst
index 6d40f1162ce480b898fa116f17ecc4458123d7e1..059fe586a0999153986feb5a61ff5a4639a896bf 100644
--- a/docs/source/algorithms/NormaliseToUnity-v1.rst
+++ b/docs/source/algorithms/NormaliseToUnity-v1.rst
@@ -29,20 +29,21 @@ Usage
    # Run algorithm
    wsNorm = NormaliseToUnity (ws)
 
-   print "Normalised Workspace"
+   print("Normalised Workspace")
    for i in range(4):
-      print "[ %.4f,%.4f,%.4f, %.4f, %.4f ]" % (wsNorm.readY(i)[0], wsNorm.readY(i)[1], 
-      wsNorm.readY(i)[2], wsNorm.readY(i)[3], wsNorm.readY(i)[4],)
+      print("[ {:.4f}, {:.4f}, {:.4f}, {:.4f}, {:.4f} ]".format(
+            wsNorm.readY(i)[0], wsNorm.readY(i)[1], wsNorm.readY(i)[2],
+            wsNorm.readY(i)[3], wsNorm.readY(i)[4]))
 
 Output:
 
 .. testoutput:: ExNormaliseToUnitySimple
 
    Normalised Workspace
-   [ 0.2239,0.0065,0.0065, 0.0065, 0.0065 ]
-   [ 0.2239,0.0065,0.0065, 0.0065, 0.0065 ]
-   [ 0.2239,0.0065,0.0065, 0.0065, 0.0065 ]
-   [ 0.2239,0.0065,0.0065, 0.0065, 0.0065 ]
+   [ 0.2239, 0.0065, 0.0065, 0.0065, 0.0065 ]
+   [ 0.2239, 0.0065, 0.0065, 0.0065, 0.0065 ]
+   [ 0.2239, 0.0065, 0.0065, 0.0065, 0.0065 ]
+   [ 0.2239, 0.0065, 0.0065, 0.0065, 0.0065 ]
 
 .. categories::
 
diff --git a/docs/source/algorithms/NormaliseVanadium-v1.rst b/docs/source/algorithms/NormaliseVanadium-v1.rst
index 0a0307f0fd22a86e72534e27cb18354e5112c277..f75bac177bb04e2417f2a5756c08011ae0737665 100644
--- a/docs/source/algorithms/NormaliseVanadium-v1.rst
+++ b/docs/source/algorithms/NormaliseVanadium-v1.rst
@@ -23,7 +23,7 @@ Usage
     inst = LoadEmptyInstrument(Filename='IDFs_for_UNIT_TESTING/MINITOPAZ_Definition.xml')
     vanadium = CreateWorkspace(DataX='0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5', DataY='10.574151,10.873,11.07348,11.22903,11.42286,11.47365,11.37375,11.112,10.512181,10.653397', UnitX='wavelength', ParentWorkspace=inst)
     norm_van = NormaliseVanadium(InputWorkspace=vanadium)
-    print "Wavelength = ", norm_van.readX(0)[2], " Y = ", norm_van.readY(0)[2]
+    print("Wavelength =  {}  Y =  {:.11f}".format(norm_van.readX(0)[2], norm_van.readY(0)[2]))
     
 Output:
 
diff --git a/docs/source/algorithms/OSIRISDiffractionReduction-v1.rst b/docs/source/algorithms/OSIRISDiffractionReduction-v1.rst
index dbcaee183c6cd7c6593e0305e815bf2b5a986f6d..87bf5b6dc369d097dea3b1e95c6ac002e6be133e 100644
--- a/docs/source/algorithms/OSIRISDiffractionReduction-v1.rst
+++ b/docs/source/algorithms/OSIRISDiffractionReduction-v1.rst
@@ -67,7 +67,7 @@ Usage
     vanadium = [os.path.join(os.path.expanduser("~"), van  + ".nxs") for van in vanadium]
     ws = OSIRISDiffractionReduction(Sample=','.join(samples), Vanadium=','.join(vanadium), CalFile="osiris_041_RES10.cal")
 
-    print "Number of Spectra: %d, Number of bins: %d" % (ws.getNumberHistograms(), ws.blocksize())
+    print("Number of Spectra: {}, Number of bins: {}".format(ws.getNumberHistograms(), ws.blocksize()))
 
 Output:
 
diff --git a/docs/source/algorithms/OneMinusExponentialCor-v1.rst b/docs/source/algorithms/OneMinusExponentialCor-v1.rst
index 873ca5a627ca7f35874a91083520a7987fe2429e..a35ca4c96320fc16fc5c2d5cc51384ed7f0d25de 100644
--- a/docs/source/algorithms/OneMinusExponentialCor-v1.rst
+++ b/docs/source/algorithms/OneMinusExponentialCor-v1.rst
@@ -31,13 +31,13 @@ Usage
 .. testcode:: ExOneMinusExp
 
     ws=CreateWorkspace([1,2,3],[1,1,1])
-    print "You can divide the data by the factor"
+    print("You can divide the data by the factor")
     wsOut=OneMinusExponentialCor(ws,2,3,"Divide")
-    print wsOut.readY(0)
+    print(wsOut.readY(0))
 
-    print "Or multiply"
+    print("Or multiply")
     wsOut=OneMinusExponentialCor(ws,2,3,"Multiply")
-    print wsOut.readY(0)
+    print(wsOut.readY(0))
 
 
 Output:
diff --git a/docs/source/algorithms/OptimizeCrystalPlacement-v1.rst b/docs/source/algorithms/OptimizeCrystalPlacement-v1.rst
index 540d3c0d13d4fb7a18cab0c786da3f1ca171f79b..24ab34c7652054f289d6d59da43ecf9d1d1507c3 100644
--- a/docs/source/algorithms/OptimizeCrystalPlacement-v1.rst
+++ b/docs/source/algorithms/OptimizeCrystalPlacement-v1.rst
@@ -56,7 +56,7 @@ Usage
     LoadIsawUB(ws,"ls5637.mat")
     wsd = OptimizeCrystalPlacement(ws)
     (wsPeakOut,fitInfoTable,chi2overDoF,nPeaks,nParams,nIndexed,covrianceInfoTable) = OptimizeCrystalPlacement(ws)
-    print "Chi2: %.4f" % chi2overDoF
+    print("Chi2: {:.4f}".format(chi2overDoF))
 
 Output:
 
diff --git a/docs/source/algorithms/OptimizeLatticeForCellType-v1.rst b/docs/source/algorithms/OptimizeLatticeForCellType-v1.rst
index 9dd20d3a95364bedf4b1a027ecc19a6d54eb597f..fb3a18d2d23c9753442f432370756fab1ccdb558 100644
--- a/docs/source/algorithms/OptimizeLatticeForCellType-v1.rst
+++ b/docs/source/algorithms/OptimizeLatticeForCellType-v1.rst
@@ -29,11 +29,11 @@ Usage
 
     ws=LoadIsawPeaks("TOPAZ_3007.peaks")
     FindUBUsingFFT(ws,MinD=8.0,MaxD=13.0)
-    print "Before Optimization:"
-    print ws.sample().getOrientedLattice().getUB()
+    print("Before Optimization:")
+    print(ws.sample().getOrientedLattice().getUB())
     OptimizeLatticeForCellType(ws,CellType="Monoclinic")
-    print "\nAfter Optimization:"
-    print ws.sample().getOrientedLattice().getUB()
+    print("\nAfter Optimization:")
+    print(ws.sample().getOrientedLattice().getUB())
 
 
 Output:
diff --git a/docs/source/algorithms/PDFFourierTransform-v1.rst b/docs/source/algorithms/PDFFourierTransform-v1.rst
index 199eaab9bd8db1d9caf6cdb9f477b2dbaca2f503..09aee1c8d5d6ef0927379848a518d521c0f6be05 100644
--- a/docs/source/algorithms/PDFFourierTransform-v1.rst
+++ b/docs/source/algorithms/PDFFourierTransform-v1.rst
@@ -144,17 +144,17 @@ Usage
 
 .. testcode:: ExPDFFouurierTransform
 
-    # Simulates Load of a workspace with all necessary parameters #################
+    # Simulates Load of a workspace with all necessary parameters
     import numpy as np;
-    xx= np.array(range(0,100))*0.1
-    yy = np.exp(-((xx)/.5)**2)
-    ws=CreateWorkspace(DataX=xx,DataY=yy,UnitX='MomentumTransfer')
-    Rt= PDFFourierTransform(ws,InputSofQType='S(Q)',PDFType='g(r)');   
-    #
+    xx = np.array(range(0,100))*0.1
+    yy = np.exp(-(2.0 * xx)**2)
+    ws = CreateWorkspace(DataX=xx, DataY=yy, UnitX='MomentumTransfer')
+    Rt = PDFFourierTransform(ws, InputSofQType='S(Q)', PDFType='g(r)')   
+
     # Look at sample results:
-    print 'part of S(Q) and its correlation function'
-    for i in xrange(0,10): 
-       print '! {0:4.2f} ! {1:5f} ! {2:f} ! {3:5f} !'.format(xx[i],yy[i],Rt.readX(0)[i],Rt.readY(0)[i])
+    print('part of S(Q) and its correlation function')
+    for i in range(10): 
+       print('! {0:4.2f} ! {1:5f} ! {2:f} ! {3:5f} !'.format(xx[i], yy[i], Rt.readX(0)[i], Rt.readY(0)[i]))