diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
new file mode 100644
index 0000000000000000000000000000000000000000..951eb39e9911ddecb55cded5ec548f190b79c40c
--- /dev/null
+++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/LoadNMoldyn3Ascii.py
@@ -0,0 +1,413 @@
+#pylint: disable=invalid-name,no-init,too-many-locals,too-many-branches
+
+from mantid.simpleapi import *
+from mantid.kernel import *
+from mantid.api import *
+
+import os
+import numpy as np
+import math
+
+
+def _split_line(a):
+    elements = a.split()  # split line on character
+    extracted = []
+    for n in elements:
+        extracted.append(float(n))
+    return extracted  # values as list
+
+
+def _find_starts(data, c, l1):
+    for l in range(l1, len(data)):
+        char = data[l]
+        if char.startswith(c):
+            line = l
+            break
+    return line
+
+
+def _find_tab_starts(data, c, l1):
+    for l in range(l1, len(data)):
+        char = data[l][1:]
+        if char.startswith(c):
+            line = l
+            break
+    return line
+
+
+def _find_ends(data, c, l1):
+    for l in range(l1, len(data)):
+        char = data[l]
+        if char.endswith(c):
+            line = l
+            break
+    return line
+
+
+def _find_char(data, c, l1):
+    for l in range(l1, len(data)):
+        char = data[l]
+        if char.find(c):
+            line = l
+            break
+    return line
+
+
+def _make_list(a, l1, l2):
+    data = ''
+    for m in range(l1, l2 + 1):
+        data += a[m]
+        alist = data.split(',')
+    return alist
+
+
+class LoadNMoldyn3Ascii(PythonAlgorithm):
+
+    _file_name = None
+    _functions = None
+    _out_ws = None
+
+
+    def category(self):
+        return 'PythonAlgorithms;Inelastic;Simulation'
+
+
+    def summary(self):
+        return 'Imports functions from CDL and ASCII files output by nMOLDYN 3.'
+
+
+    def PyInit(self):
+        self.declareProperty(FileProperty('Filename', '',
+                                          action=FileAction.Load,
+                                          extensions=['.cdl', '.dat']),
+                                          doc='File path for data')
+
+        self.declareProperty(StringArrayProperty('Functions'),
+                             doc='Names of functions to attempt to load from file')
+
+        self.declareProperty(WorkspaceProperty('OutputWorkspace', '',
+                                               direction=Direction.Output),
+                             doc='Output workspace name')
+
+
+    def validateInputs(self):
+        issues = dict()
+
+        sample_filename = self.getPropertyValue('Filename')
+        file_type = os.path.splitext(sample_filename)[1]
+        function_list = self.getProperty('Functions').value
+
+        if len(function_list) == 0 and file_type == '.cdl':
+            issues['Functions'] = 'Must specify at least one function when loading a CDL file'
+
+        if len(function_list) > 0 and file_type == '.dat':
+            issues['Functions'] = 'Cannot specify functions when loading an ASCII file'
+
+        return issues
+
+
+    def PyExec(self):
+        # Do setup
+        self._setup()
+
+        # Load data file
+        data, name, ext = self._load_file()
+
+        # Run nMOLDYN import
+        if ext == 'cdl':
+            self._cdl_import(data, name)
+        elif ext == 'dat':
+            self._ascii_import(data, name)
+        else:
+            raise RuntimeError('Unrecognised file format: %s' % ext)
+
+        # Set the output workspace
+        self.setProperty('OutputWorkspace', self._out_ws)
+
+
+    def _setup(self):
+        """
+        Gets algorithm properties.
+        """
+        self._file_name = self.getPropertyValue('Filename')
+        self._out_ws = self.getPropertyValue('OutputWorkspace')
+
+        raw_functions = self.getProperty('Functions').value
+        self._functions = [x.strip() for x in raw_functions]
+
+
+    def _load_file(self):
+        """
+        Attempts to load the sample file.
+
+        @returns A tuple with the ASCII data, sample file name and file extension
+        """
+
+        # Get some data about the file
+        path = self._file_name
+        base = os.path.basename(path)
+        name = os.path.splitext(base)[0]
+        ext = os.path.splitext(path)[1]
+
+        # Remove dot from extension
+        if len(ext) > 1:
+            ext = ext[1:]
+
+        logger.debug('Base filename for %s is %s' % (self._file_name, name))
+        logger.debug('File type of %s is %s' % (self._file_name, ext))
+
+        if not os.path.isfile(path):
+            path = FileFinder.getFullPath(path)
+
+        logger.information('Got file path for %s: %s' % (self._file_name, path))
+
+        # Open file and get data
+        handle = open(path, 'r')
+        data = []
+        for line in handle:
+            line = line.rstrip()
+            data.append(line)
+        handle.close()
+
+        return data, name, ext
+
+
+    def _find_dimensions(self, data):
+        """
+        Gets the number of Q, time and frequency values in given raw data.
+
+        @param data Raw data to search
+        """
+
+        num_q_values = _find_tab_starts(data, 'NQVALUES', 0)
+        num_time_values = _find_tab_starts(data, 'NTIMES', 0)
+        num_freq_values = _find_tab_starts(data, 'NFREQUENCIES', 0)
+
+        q_el = data[num_q_values].split()
+        num_q = int(q_el[2])
+        t_el = data[num_time_values].split()
+        num_t = int(t_el[2])
+        f_el = data[num_freq_values].split()
+        num_f = int(f_el[2])
+
+        logger.debug(data[2][1:-1])
+        logger.debug(data[3][1:-1])
+        logger.debug(data[6][1:-1])
+
+        return num_q, num_t, num_f
+
+
+    def _cdl_import(self, data, name):
+        """
+        Import data from CDL file.
+
+        @param data Raw data
+        @param name Name of data file
+        """
+
+        logger.notice('Loading CDL file: %s' % name)
+
+        len_data = len(data)
+
+        # raw head
+        nQ, nT, nF = self._find_dimensions(data)
+        ldata = _find_starts(data, 'data:', 0)
+        lq1 = _find_starts(data, ' q =', ldata)  # start Q values
+        lq2 = _find_starts(data, ' q =', lq1 - 1)
+        Qlist = _make_list(data, lq1, lq2)
+        if nQ != len(Qlist):
+            raise RUntimeError('Error reading Q values')
+        Qf = Qlist[0].split()
+        Q = [float(Qf[2]) / 10.0]
+        for m in range(1, nQ - 1):
+            Q.append(float(Qlist[m]) / 10.0)
+
+        Q.append(float(Qlist[nQ - 1][:-1]) / 10.0)
+        logger.information('Q values = ' + str(Q))
+
+        lt1 = _find_starts(data, ' time =', lq2)  # start T values
+        lt2 = _find_ends(data, ';', lt1)
+        Tlist = _make_list(data, lt1, lt2)
+        if nT != len(Tlist):
+            raise RuntimeError('Error reading Time values')
+
+        Tf = Tlist[0].split()
+        T = [float(Tf[2])]
+        for m in range(1, nT - 1):
+            T.append(float(Tlist[m]))
+
+        T.append(float(Tlist[nT - 1][:-1]))
+        T.append(2 * T[nT - 1] - T[nT - 2])
+        logger.information('T values = ' + str(T[:2]) + ' to ' + str(T[-3:]))
+
+        lf1 = _find_starts(data, ' frequency =', lq2)  # start F values
+        lf2 = _find_ends(data, ';', lf1)
+        Flist = _make_list(data, lf1, lf2)
+        if nF != len(Flist):
+            raise RuntimeError('Error reading Freq values')
+
+        Ff = Flist[0].split()
+        F = [float(Ff[2])]
+        for m in range(1, nF - 1):
+            F.append(float(Flist[m]))
+
+        F.append(float(Flist[nF - 1][:-1]))
+        F.append(2 * F[nF - 1] - T[nF - 2])
+        logger.information('F values = ' + str(F[:2]) + ' to ' + str(F[-3:]))
+
+        # Function
+        output_ws_list = list()
+        for func in self._functions:
+            start = []
+            lstart = lt2
+            if func[:3] == 'Fqt':
+                nP = nT
+                xEn = np.array(T)
+                eZero = np.zeros(nT)
+                xUnit = 'TOF'
+            elif func[:3] == 'Sqw':
+                nP = nF
+                xEn = np.array(F)
+                eZero = np.zeros(nF)
+                xUnit = 'Energy'
+            else:
+                raise RuntimeError('Failed to parse function string ' + func)
+
+            for n in range(0, nQ):
+                for m in range(lstart, len_data):
+                    char = data[m]
+                    if char.startswith('  // ' + func):
+                        start.append(m)
+                        lstart = m + 1
+            lend = _find_ends(data, ';', lstart)
+            start.append(lend + 1)
+
+            # Throw error if we couldn't find the function
+            if len(start) < 2:
+                raise RuntimeError('Failed to parse function string ' + func)
+
+            Qaxis = ''
+            for n in range(0, nQ):
+                logger.information(str(start))
+                logger.information('Reading : ' + data[start[n]])
+
+                Slist = _make_list(data, start[n] + 1, start[n + 1] - 1)
+                if n == nQ - 1:
+                    Slist[nP - 1] = Slist[nP - 1][:-1]
+                S = []
+                for m in range(0, nP):
+                    S.append(float(Slist[m]))
+                if nP != len(S):
+                    raise RuntimeError('Error reading S values')
+                else:
+                    logger.information('S values = ' + str(S[:2]) + ' to ' + str(S[-2:]))
+                if n == 0:
+                    Qaxis += str(Q[n])
+                    xDat = xEn
+                    yDat = np.array(S)
+                    eDat = eZero
+                else:
+                    Qaxis += ',' + str(Q[n])
+                    xDat = np.append(xDat, xEn)
+                    yDat = np.append(yDat, np.array(S))
+                    eDat = np.append(eDat, eZero)
+            outWS = name + '_' + func
+
+            CreateWorkspace(OutputWorkspace=outWS,
+                            DataX=xDat,
+                            DataY=yDat,
+                            DataE=eDat,
+                            Nspec=nQ,
+                            UnitX=xUnit,
+                            VerticalAxisUnit='MomentumTransfer',
+                            VerticalAxisValues=Qaxis)
+
+            output_ws_list.append(outWS)
+
+        GroupWorkspaces(InputWorkspaces=output_ws_list,
+                        OutputWorkspace=self._out_ws)
+
+
+    def _ascii_import(self, data, name):
+        """
+        Import ASCII data.
+
+        @param data Raw ASCII data
+        @param name Name of data file
+        """
+
+        from IndirectCommon import getEfixed
+        from IndirectNeutron import ChangeAngles, InstrParas
+
+        logger.notice('Loading ASCII data: %s' % name)
+
+        val = _split_line(data[3])
+        Q = []
+        for n in range(1, len(val)):
+            Q.append(val[n])
+
+        nQ = len(Q)
+        x = []
+        y = []
+        for n in range(4, len(data)):
+            val = _split_line(data[n])
+            x.append(val[0])
+            yval = val[1:]
+            y.append(yval)
+
+        nX = len(x)
+        logger.information('nQ = ' + str(nQ))
+        logger.information('nT = ' + str(nX))
+
+        xT = np.array(x)
+        eZero = np.zeros(nX)
+        #Qaxis = ''
+        for m in range(0, nQ):
+            logger.information('Q[' + str(m + 1) + '] : ' + str(Q[m]))
+
+            S = []
+            for n in range(0, nX):
+                S.append(y[n][m])
+
+            if m == 0:
+                #Qaxis += str(Q[m])
+                xDat = xT
+                yDat = np.array(S)
+                eDat = eZero
+            else:
+                #Qaxis += ',' + str(Q[m])
+                xDat = np.append(xDat, xT)
+                yDat = np.append(yDat, np.array(S))
+                eDat = np.append(eDat, eZero)
+
+        CreateWorkspace(OutputWorkspace=self._out_ws,
+                        DataX=xDat,
+                        DataY=yDat,
+                        DataE=eDat,
+                        Nspec=nQ,
+                        UnitX='TOF')
+
+        Qmax = Q[nQ - 1]
+        instr = 'MolDyn'
+        ana = 'simul'
+        if Qmax <= 2.0:
+            refl = '2'
+        else:
+            refl = '4'
+
+        InstrParas(self._out_ws, instr, ana, refl)
+        efixed = getEfixed(self._out_ws)
+        logger.information('Qmax = ' + str(Qmax) + ' ; efixed = ' + str(efixed))
+        pi4 = 4.0 * math.pi
+        wave = 1.8 * math.sqrt(25.2429 / efixed)
+        theta = []
+        for n in range(0, nQ):
+            qw = wave * Q[n] / pi4
+            ang = 2.0 * math.degrees(math.asin(qw))
+            theta.append(ang)
+
+        ChangeAngles(self._out_ws, instr, theta)
+
+
+# Register algorithm with Mantid
+AlgorithmFactory.subscribe(LoadNMoldyn3Ascii)
diff --git a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MolDyn.py b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MolDyn.py
index 8b6db2241b2d7e72f0ceda0a3bc106590064cc67..6706a380b8ff3e93a011e0a2ba84ad4ac9fa95e6 100644
--- a/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MolDyn.py
+++ b/Code/Mantid/Framework/PythonInterface/plugins/algorithms/WorkflowAlgorithms/MolDyn.py
@@ -4,61 +4,8 @@ from mantid.kernel import *
 from mantid.api import *
 
 import os
-import numpy as np
-import math
 
 
-def _split_line(a):
-    elements = a.split()  # split line on character
-    extracted = []
-    for n in elements:
-        extracted.append(float(n))
-    return extracted  # values as list
-
-
-def _find_starts(data, c, l1):
-    for l in range(l1, len(data)):
-        char = data[l]
-        if char.startswith(c):
-            line = l
-            break
-    return line
-
-
-def _find_tab_starts(data, c, l1):
-    for l in range(l1, len(data)):
-        char = data[l][1:]
-        if char.startswith(c):
-            line = l
-            break
-    return line
-
-
-def _find_ends(data, c, l1):
-    for l in range(l1, len(data)):
-        char = data[l]
-        if char.endswith(c):
-            line = l
-            break
-    return line
-
-
-def _find_char(data, c, l1):
-    for l in range(l1, len(data)):
-        char = data[l]
-        if char.find(c):
-            line = l
-            break
-    return line
-
-
-def _make_list(a, l1, l2):
-    data = ''
-    for m in range(l1, l2 + 1):
-        data += a[m]
-        alist = data.split(',')
-    return alist
-
 #pylint: disable=too-many-instance-attributes
 class MolDyn(PythonAlgorithm):
 
@@ -72,16 +19,19 @@ class MolDyn(PythonAlgorithm):
     _out_ws = None
     _mtd_plot = None
 
+
     def category(self):
         return 'Workflow\\Inelastic;PythonAlgorithms;Inelastic;Simulation'
 
+
     def summary(self):
-        return 'Imports nMOLDYN simulations from CDL and ASCII files.'
+        return 'Imports and processes simulated functions from nMOLDYN.'
+
 
     def PyInit(self):
         self.declareProperty(FileProperty('Filename', '',
-                                          action=FileAction.OptionalLoad,
-                                          extensions=['.cdl', '.dat']),\
+                                          action=FileAction.Load,
+                                          extensions=['.cdl', '.dat']),
                                           doc='File path for data')
 
         self.declareProperty(StringArrayProperty('Functions'),
@@ -110,18 +60,10 @@ class MolDyn(PythonAlgorithm):
     def validateInputs(self):
         issues = dict()
 
-        sample_filename = self.getPropertyValue('Filename')
-        function_list = self.getProperty('Functions').value
         symm = self.getProperty('SymmetriseEnergy').value
         res_ws = self.getPropertyValue('Resolution')
         e_max = self.getPropertyValue('MaxEnergy')
 
-        if len(function_list) == 0 and os.path.splitext(sample_filename)[1] == 'cdl':
-            issues['Functions'] = 'Must specify at least one function when loading a CDL file'
-
-        if len(function_list) > 0 and os.path.splitext(sample_filename)[1] == 'dat':
-            issues['Functions'] = 'Cannot specify functions when loading an ASCII file'
-
         if res_ws is not '' and e_max is '':
             issues['MaxEnergy'] = 'MaxEnergy must be set when convolving with an instrument resolution'
 
@@ -130,25 +72,19 @@ class MolDyn(PythonAlgorithm):
 
         return issues
 
+
     #pylint: disable=too-many-branches
     def PyExec(self):
         from IndirectImport import import_mantidplot
-
         self._mtd_plot = import_mantidplot()
 
         # Do setup
         self._setup()
 
-        # Load data file
-        data, name, ext = self._load_file()
-
         # Run nMOLDYN import
-        if ext == 'cdl':
-            self._cdl_import(data, name)
-        elif ext == 'dat':
-            self._ascii_import(data, name)
-        else:
-            raise RuntimeError('Unrecognised file format: %s' % ext)
+        LoadNMoldyn3Ascii(Filename=self.getPropertyValue('Filename'),
+                          OutputWorkspace=self._out_ws,
+                          Functions=self.getPropertyValue('Functions'))
 
         # Do processing specific to workspaces in energy
         if isinstance(mtd[self._out_ws], WorkspaceGroup):
@@ -224,12 +160,8 @@ class MolDyn(PythonAlgorithm):
         self._plot = self.getProperty('Plot').value
         self._save = self.getProperty('Save').value
 
-        self._sam_path = self.getPropertyValue('Filename')
         self._symmetrise = self.getProperty('SymmetriseEnergy').value
 
-        raw_functions = self.getProperty('Functions').value
-        self._functions = [x.strip() for x in raw_functions]
-
         emax_str = self.getPropertyValue('MaxEnergy')
         self._emax = None
         if emax_str != '':
@@ -239,280 +171,6 @@ class MolDyn(PythonAlgorithm):
         self._out_ws = self.getPropertyValue('OutputWorkspace')
 
 
-    def _load_file(self):
-        """
-        Attempts to load the sample file.
-
-        @returns A tuple with the ASCII data, sample file name and file extension
-        """
-
-        # Get some data about the file
-        path = self._sam_path
-        base = os.path.basename(path)
-        name = os.path.splitext(base)[0]
-        ext = os.path.splitext(path)[1]
-
-        # Remove dot from extension
-        if len(ext) > 1:
-            ext = ext[1:]
-
-        logger.debug('Base filename for %s is %s' % (self._sam_path, name))
-        logger.debug('File type of %s is %s' % (self._sam_path, ext))
-
-        if not os.path.isfile(path):
-            path = FileFinder.getFullPath(path)
-
-        logger.information('Got file path for %s: %s' % (self._sam_path, path))
-
-        # Open file and get data
-        handle = open(path, 'r')
-        data = []
-        for line in handle:
-            line = line.rstrip()
-            data.append(line)
-        handle.close()
-
-        return data, name, ext
-
-
-    def _find_dimensions(self, data):
-        """
-        Gets the number of Q, time and frequency values in given raw data.
-
-        @param data Raw data to search
-        """
-
-        num_q_values = _find_tab_starts(data, 'NQVALUES', 0)
-        num_time_values = _find_tab_starts(data, 'NTIMES', 0)
-        num_freq_values = _find_tab_starts(data, 'NFREQUENCIES', 0)
-
-        q_el = data[num_q_values].split()
-        num_q = int(q_el[2])
-        t_el = data[num_time_values].split()
-        num_t = int(t_el[2])
-        f_el = data[num_freq_values].split()
-        num_f = int(f_el[2])
-
-        logger.debug(data[2][1:-1])
-        logger.debug(data[3][1:-1])
-        logger.debug(data[6][1:-1])
-
-        return num_q, num_t, num_f
-
-
-    #pylint: disable=too-many-locals,too-many-branches
-    def _cdl_import(self, data, name):
-        """
-        Import data from CDL file.
-
-        @param data Raw data
-        @param name Name of data file
-        """
-
-        logger.notice('Loading CDL file: %s' % name)
-
-        len_data = len(data)
-
-        # raw head
-        nQ, nT, nF = self._find_dimensions(data)
-        ldata = _find_starts(data, 'data:', 0)
-        lq1 = _find_starts(data, ' q =', ldata)  # start Q values
-        lq2 = _find_starts(data, ' q =', lq1 - 1)
-        Qlist = _make_list(data, lq1, lq2)
-        if nQ != len(Qlist):
-            raise RUntimeError('Error reading Q values')
-        Qf = Qlist[0].split()
-        Q = [float(Qf[2]) / 10.0]
-        for m in range(1, nQ - 1):
-            Q.append(float(Qlist[m]) / 10.0)
-
-        Q.append(float(Qlist[nQ - 1][:-1]) / 10.0)
-        logger.information('Q values = ' + str(Q))
-
-        lt1 = _find_starts(data, ' time =', lq2)  # start T values
-        lt2 = _find_ends(data, ';', lt1)
-        Tlist = _make_list(data, lt1, lt2)
-        if nT != len(Tlist):
-            raise RuntimeError('Error reading Time values')
-
-        Tf = Tlist[0].split()
-        T = [float(Tf[2])]
-        for m in range(1, nT - 1):
-            T.append(float(Tlist[m]))
-
-        T.append(float(Tlist[nT - 1][:-1]))
-        T.append(2 * T[nT - 1] - T[nT - 2])
-        logger.information('T values = ' + str(T[:2]) + ' to ' + str(T[-3:]))
-
-        lf1 = _find_starts(data, ' frequency =', lq2)  # start F values
-        lf2 = _find_ends(data, ';', lf1)
-        Flist = _make_list(data, lf1, lf2)
-        if nF != len(Flist):
-            raise RuntimeError('Error reading Freq values')
-
-        Ff = Flist[0].split()
-        F = [float(Ff[2])]
-        for m in range(1, nF - 1):
-            F.append(float(Flist[m]))
-
-        F.append(float(Flist[nF - 1][:-1]))
-        F.append(2 * F[nF - 1] - T[nF - 2])
-        logger.information('F values = ' + str(F[:2]) + ' to ' + str(F[-3:]))
-
-        # Function
-        output_ws_list = list()
-        for func in self._functions:
-            start = []
-            lstart = lt2
-            if func[:3] == 'Fqt':
-                nP = nT
-                xEn = np.array(T)
-                eZero = np.zeros(nT)
-                xUnit = 'TOF'
-            elif func[:3] == 'Sqw':
-                nP = nF
-                xEn = np.array(F)
-                eZero = np.zeros(nF)
-                xUnit = 'Energy'
-            else:
-                raise RuntimeError('Failed to parse function string ' + func)
-
-            for n in range(0, nQ):
-                for m in range(lstart, len_data):
-                    char = data[m]
-                    if char.startswith('  // ' + func):
-                        start.append(m)
-                        lstart = m + 1
-            lend = _find_ends(data, ';', lstart)
-            start.append(lend + 1)
-
-            # Throw error if we couldn't find the function
-            if len(start) < 2:
-                raise RuntimeError('Failed to parse function string ' + func)
-
-            Qaxis = ''
-            for n in range(0, nQ):
-                logger.information(str(start))
-                logger.information('Reading : ' + data[start[n]])
-
-                Slist = _make_list(data, start[n] + 1, start[n + 1] - 1)
-                if n == nQ - 1:
-                    Slist[nP - 1] = Slist[nP - 1][:-1]
-                S = []
-                for m in range(0, nP):
-                    S.append(float(Slist[m]))
-                if nP != len(S):
-                    raise RuntimeError('Error reading S values')
-                else:
-                    logger.information('S values = ' + str(S[:2]) + ' to ' + str(S[-2:]))
-                if n == 0:
-                    Qaxis += str(Q[n])
-                    xDat = xEn
-                    yDat = np.array(S)
-                    eDat = eZero
-                else:
-                    Qaxis += ',' + str(Q[n])
-                    xDat = np.append(xDat, xEn)
-                    yDat = np.append(yDat, np.array(S))
-                    eDat = np.append(eDat, eZero)
-            outWS = name + '_' + func
-
-            CreateWorkspace(OutputWorkspace=outWS,
-                            DataX=xDat,
-                            DataY=yDat,
-                            DataE=eDat,
-                            Nspec=nQ,
-                            UnitX=xUnit,
-                            VerticalAxisUnit='MomentumTransfer',
-                            VerticalAxisValues=Qaxis)
-
-            output_ws_list.append(outWS)
-
-        GroupWorkspaces(InputWorkspaces=output_ws_list,
-                        OutputWorkspace=self._out_ws)
-
-
-    #pylint: disable=too-many-locals
-    def _ascii_import(self, data, name):
-        """
-        Import ASCII data.
-
-        @param data Raw ASCII data
-        @param name Name of data file
-        """
-
-        from IndirectCommon import getEfixed
-        from IndirectNeutron import ChangeAngles, InstrParas
-
-        logger.notice('Loading ASCII data: %s' % name)
-
-        val = _split_line(data[3])
-        Q = []
-        for n in range(1, len(val)):
-            Q.append(val[n])
-
-        nQ = len(Q)
-        x = []
-        y = []
-        for n in range(4, len(data)):
-            val = _split_line(data[n])
-            x.append(val[0])
-            yval = val[1:]
-            y.append(yval)
-
-        nX = len(x)
-        logger.information('nQ = ' + str(nQ))
-        logger.information('nT = ' + str(nX))
-
-        xT = np.array(x)
-        eZero = np.zeros(nX)
-        #Qaxis = ''
-        for m in range(0, nQ):
-            logger.information('Q[' + str(m + 1) + '] : ' + str(Q[m]))
-
-            S = []
-            for n in range(0, nX):
-                S.append(y[n][m])
-
-            if m == 0:
-                #Qaxis += str(Q[m])
-                xDat = xT
-                yDat = np.array(S)
-                eDat = eZero
-            else:
-                #Qaxis += ',' + str(Q[m])
-                xDat = np.append(xDat, xT)
-                yDat = np.append(yDat, np.array(S))
-                eDat = np.append(eDat, eZero)
-
-        CreateWorkspace(OutputWorkspace=self._out_ws,
-                        DataX=xDat,
-                        DataY=yDat,
-                        DataE=eDat,
-                        Nspec=nQ, UnitX='TOF')
-
-        Qmax = Q[nQ - 1]
-        instr = 'MolDyn'
-        ana = 'simul'
-        if Qmax <= 2.0:
-            refl = '2'
-        else:
-            refl = '4'
-
-        InstrParas(self._out_ws, instr, ana, refl)
-        efixed = getEfixed(self._out_ws)# RunParas(self._out_ws, instr, name, name)
-        logger.information('Qmax = ' + str(Qmax) + ' ; efixed = ' + str(efixed))
-        pi4 = 4.0 * math.pi
-        wave = 1.8 * math.sqrt(25.2429 / efixed)
-        theta = []
-        for n in range(0, nQ):
-            qw = wave * Q[n] / pi4
-            ang = 2.0 * math.degrees(math.asin(qw))
-            theta.append(ang)
-
-        ChangeAngles(self._out_ws, instr, theta)
-
-
     def _create_res_ws(self, num_sample_hist):
         """
         Creates a resolution workspace.
diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
index de20c3b1338580fe65312b5b74d9c29ff1d6a9f4..74c699a4f65e6bdf71457f24c71ef133409075c1 100644
--- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
+++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/CMakeLists.txt
@@ -41,6 +41,7 @@ set ( TEST_PY_FILES
   LoadLiveDataTest.py
   LoadLogPropertyTableTest.py
   LoadMultipleGSSTest.py
+  LoadNMoldyn3AsciiTest.py
   MaskAngleTest.py
   MaskBTPTest.py
   MaskWorkspaceToCalFileTest.py
diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd4868a60cd8526927b7b9bd4c960fd0bc1171f0
--- /dev/null
+++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/LoadNMoldyn3AsciiTest.py
@@ -0,0 +1,92 @@
+#pylint: disable=too-many-public-methods
+
+import unittest
+from mantid.simpleapi import *
+from mantid.api import *
+
+class MolDynTest(unittest.TestCase):
+
+    _cdl_filename = 'NaF_DISF.cdl'
+    _dat_filename = 'WSH_test.dat'
+
+
+    def test_load_fqt_from_cdl(self):
+        """
+        Load an F(Q, t) function from an nMOLDYN 3 .cdl file
+        """
+
+        moldyn_group = LoadNMoldyn3Ascii(Filename=self._cdl_filename,
+                                         Functions=['Fqt-total'],
+                                         OutputWorkspace='__LoadNMoldyn3Ascii_test')
+
+        self.assertTrue(isinstance(moldyn_group, WorkspaceGroup))
+        self.assertEqual(len(moldyn_group), 1)
+        self.assertEqual(moldyn_group[0].name(), 'NaF_DISF_Fqt-total')
+
+        iqt_ws = moldyn_group[0]
+        self.assertTrue(isinstance(iqt_ws, MatrixWorkspace))
+        self.assertTrue(iqt_ws.getNumberHistograms(), 1)
+
+        # X axis should be in TOF for an Fqt function
+        units = iqt_ws.getAxis(0).getUnit().unitID()
+        self.assertEqual(units, 'TOF')
+
+
+    def test_load_sqw_from_cdl(self):
+        """
+        Load an S(Q, w) function from a nMOLDYN 3 .cdl file
+        """
+
+        moldyn_group = LoadNMoldyn3Ascii(Filename=self._cdl_filename,
+                                         Functions=['Sqw-total'],
+                                         OutputWorkspace='__LoadNMoldyn3Ascii_test')
+
+        self.assertTrue(isinstance(moldyn_group, WorkspaceGroup))
+        self.assertEqual(moldyn_group[0].name(), 'NaF_DISF_Sqw-total')
+
+        sqw_ws = moldyn_group[0]
+        self.assertTrue(isinstance(sqw_ws, MatrixWorkspace))
+        self.assertTrue(sqw_ws.getNumberHistograms(), 1)
+
+        # X axis should be in Energy for an Sqw function
+        units = sqw_ws.getAxis(0).getUnit().unitID()
+        self.assertEqual(units, 'Energy')
+
+
+    def test_load_from_dat(self):
+        """
+        Load a function from an nMOLDYN 3 .dat file
+        """
+
+        moldyn_ws = LoadNMoldyn3Ascii(Filename=self._dat_filename,
+                                      OutputWorkspace='__LoadNMoldyn3Ascii_test')
+
+        self.assertTrue(isinstance(moldyn_ws, MatrixWorkspace))
+        self.assertTrue(moldyn_ws.getNumberHistograms(), 12)
+
+
+    def test_function_validation_cdl(self):
+        """
+        Tests that the algorithm cannot be run when no functions are specified
+        when loading a .cdl file.
+        """
+        self.assertRaises(RuntimeError,
+                          LoadNMoldyn3Ascii,
+                          Filename=self._cdl_filename,
+                          OutputWorkspace='__LoadNMoldyn3Ascii_test')
+
+
+    def test_function_validation_dat(self):
+        """
+        Tests that the algorithm cannot be run when functions are specified
+        when loading a .dat file.
+        """
+        self.assertRaises(RuntimeError,
+                          LoadNMoldyn3Ascii,
+                          Filename=self._dat_filename,
+                          Functions=['Sqw-total'],
+                          OutputWorkspace='__LoadNMoldyn3Ascii_test')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/MolDynTest.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/MolDynTest.py
index 2c2cee3d3004128baf9c93574242f5e520bcec38..f9c57950dabcdc7dc1399543c88d8f2d2bb34dd5 100644
--- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/MolDynTest.py
+++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/MolDynTest.py
@@ -1,3 +1,5 @@
+#pylint: disable=too-many-public-methods,invalid-name
+
 import unittest
 from mantid.simpleapi import *
 from mantid.api import *
@@ -5,39 +7,18 @@ from mantid.api import *
 
 class MolDynTest(unittest.TestCase):
 
-    def test_loadFqt(self):
-        # Load an Fwt function from a nMOLDYN file
-        moldyn_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Fqt-total'])
-
-        # A workspace should exist in this name format
-        self.assertTrue(mtd.doesExist('NaF_DISF_Fqt-total'))
-
-        # X axis should be in TOF for an Fqt function
-        units = mtd['NaF_DISF_Fqt-total'].getAxis(0).getUnit().unitID()
-        self.assertEqual(units, 'TOF')
-
-
-    def test_loadSqw(self):
-        # Load an Sqw function from a nMOLDYN file
-        moldyn_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Sqw-total'])
-
-        # A workspace should exist in this name format
-        self.assertTrue(mtd.doesExist('NaF_DISF_Sqw-total'))
-
-        # X axis should be in Energy for an Sqw function
-        units = mtd['NaF_DISF_Sqw-total'].getAxis(0).getUnit().unitID()
-        self.assertEqual(units, 'Energy')
-
-
     def test_loadSqwWithEMax(self):
         # Load an Sqw function from a nMOLDYN file
-        moldyn_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Sqw-total'], MaxEnergy="1.0")
+        moldyn_group = MolDyn(Filename='NaF_DISF.cdl',
+                              Functions=['Sqw-total'],
+                              MaxEnergy="1.0")
 
-        # A workspace should exist in this name format
-        self.assertTrue(mtd.doesExist('NaF_DISF_Sqw-total'))
+        self.assertTrue(isinstance(moldyn_group, WorkspaceGroup))
+        self.assertEqual(len(moldyn_group), 1)
+        self.assertEqual(moldyn_group[0].name(), 'NaF_DISF_Sqw-total')
 
         # Get max enery from result workspace
-        x_data = mtd['NaF_DISF_Sqw-total'].dataX(0)
+        x_data = moldyn_group[0].dataX(0)
         x_max = x_data[len(x_data) - 1]
 
         # Check that it is less that what was passed to algorithm
@@ -46,13 +27,16 @@ class MolDynTest(unittest.TestCase):
 
     def test_loadSqwWithSymm(self):
         # Load an Sqw function from a nMOLDYN file
-        moldyn_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Sqw-total'], SymmetriseEnergy=True)
+        moldyn_group = MolDyn(Filename='NaF_DISF.cdl',
+                              Functions=['Sqw-total'],
+                              SymmetriseEnergy=True)
 
-        # A workspace should exist in this name format
-        self.assertTrue(mtd.doesExist('NaF_DISF_Sqw-total'))
+        self.assertTrue(isinstance(moldyn_group, WorkspaceGroup))
+        self.assertEqual(len(moldyn_group), 1)
+        self.assertEqual(moldyn_group[0].name(), 'NaF_DISF_Sqw-total')
 
         # Get max and min energy from result workspace
-        x_data = mtd['NaF_DISF_Sqw-total'].dataX(0)
+        x_data = moldyn_group[0].dataX(0)
         x_max = x_data[len(x_data) - 1]
         x_min = x_data[0]
 
@@ -62,13 +46,23 @@ class MolDynTest(unittest.TestCase):
 
     def test_loadSqwWithRes(self):
         # Create a sample workspace thet looks like an instrument resolution
-        sample_res = CreateSampleWorkspace(NumBanks=1, BankPixelWidth=1, XUnit='Energy', XMin=-10, XMax=10, BinWidth=0.1)
+        sample_res = CreateSampleWorkspace(NumBanks=1,
+                                           BankPixelWidth=1,
+                                           XUnit='Energy',
+                                           XMin=-10,
+                                           XMax=10,
+                                           BinWidth=0.1)
 
         # Load an Sqw function from a nMOLDYN file
-        moldyn_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Sqw-total'], MaxEnergy="1.0", SymmetriseEnergy=True, Resolution=sample_res)
+        moldyn_group = MolDyn(Filename='NaF_DISF.cdl',
+                              Functions=['Sqw-total'],
+                              MaxEnergy="1.0",
+                              SymmetriseEnergy=True,
+                              Resolution=sample_res)
 
-        # A workspace should exist in this name format
-        self.assertTrue(mtd.doesExist('NaF_DISF_Sqw-total'))
+        self.assertTrue(isinstance(moldyn_group, WorkspaceGroup))
+        self.assertEqual(len(moldyn_group), 1)
+        self.assertEqual(moldyn_group[0].name(), 'NaF_DISF_Sqw-total')
 
 
     def test_loadSqwWithResWithNoEMaxFails(self):
@@ -77,11 +71,19 @@ class MolDynTest(unittest.TestCase):
         """
 
         # Create a sample workspace thet looks like an instrument resolution
-        sample_res = CreateSampleWorkspace(NumBanks=1, BankPixelWidth=1, XUnit='Energy', XMin=-10, XMax=10, BinWidth=0.1)
+        sample_res = CreateSampleWorkspace(NumBanks=1,
+                                           BankPixelWidth=1,
+                                           XUnit='Energy',
+                                           XMin=-10,
+                                           XMax=10,
+                                           BinWidth=0.1)
 
         # Load an Sqw function from a nMOLDYN file
         self.assertRaises(RuntimeError, MolDyn,
-                          Filename='NaF_DISF.cdl', Functions=['Sqw-total'], Resolution=sample_res, OutputWorkspace='moldyn_group')
+                          Filename='NaF_DISF.cdl',
+                          Functions=['Sqw-total'],
+                          Resolution=sample_res,
+                          OutputWorkspace='moldyn_group')
 
 
 if __name__ == '__main__':
diff --git a/Code/Mantid/Testing/Data/UnitTest/WSH_test.dat.md5 b/Code/Mantid/Testing/Data/UnitTest/WSH_test.dat.md5
new file mode 100644
index 0000000000000000000000000000000000000000..fc92d56413d3038a0ef2b00f75682ef30e099578
--- /dev/null
+++ b/Code/Mantid/Testing/Data/UnitTest/WSH_test.dat.md5
@@ -0,0 +1 @@
+5985c297c7b330280d9c64a9d12605d0
\ No newline at end of file
diff --git a/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst b/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e8e00ea6fb3d8e8ad9fad75c664ae3b6bc8fb0b9
--- /dev/null
+++ b/Code/Mantid/docs/source/algorithms/LoadNMoldyn3Ascii-v1.rst
@@ -0,0 +1,40 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+This algorithm is a loader from simulation data saved by version 3 of the
+nMOLDYN package, simulations can be loaded form either a plain ASCII (``.dat``)
+file or ``.cdl`` file.
+
+Usage
+-----
+
+.. include:: ../usagedata-note.txt
+
+**Example - Loading a simulation from a CDL file.**
+
+.. testcode:: ExLoadCDLFile
+
+    out_ws_group = LoadNMoldyn3Ascii(Filename='NaF_DISF.cdl',
+                                     Functions=['Fqt-total', 'Sqw-total'])
+
+    for ws_name in out_ws_group.getNames():
+      print ws_name
+
+Output:
+
+.. testoutput:: ExLoadCDLFile
+
+    NaF_DISF_Fqt-total
+    NaF_DISF_Sqw-total
+
+.. categories::
+
+.. sourcelink::
diff --git a/Code/Mantid/docs/source/algorithms/MolDyn-v1.rst b/Code/Mantid/docs/source/algorithms/MolDyn-v1.rst
index dfd8ce920b023e3169401e457ba881dfad6a6486..22b793f2b1bfcf469b2a12a32dc6ebcd649ef601 100644
--- a/Code/Mantid/docs/source/algorithms/MolDyn-v1.rst
+++ b/Code/Mantid/docs/source/algorithms/MolDyn-v1.rst
@@ -9,27 +9,32 @@
 Description
 -----------
 
-This algorithm is a loader from simulation data from the nMOLDYN software,
-simulations can be loaded form either a plain ASCII file or CDL file.
+This algorithm is used to load and process simualtion data from the nMOLDYN
+package.
 
-When loading from a CDL file one or multiple functions can be loaded, when
-loading a single function an instrument resolution workspace can be provided
-which the loaded function is convoluted with to allow comparison with actual
-instrument data.
+Currently this supports loading the ``.cdl`` and ``.dat`` files created by
+version 3 of nMOLDYN (using the :ref:`LoadNMoldyn3Ascii
+<algm-LoadNMoldyn3Ascii>` algorithm).
 
-When loading from ASCII function selection and convolution with an instrument
-resolution are unavailable.
+When loading from a ``.cdl`` file from nMOLDYN 3 one or multiple functions can
+be loaded, when loading a single function an instrument resolution workspace can
+be provided which the loaded function is convoluted with to allow comparison
+with actual instrument data.
+
+When loading from a ``.dat`` file from nMOLDYN 3 function selection and
+convolution with an instrument resolution are unavailable.
 
 Usage
 -----
 
 .. include:: ../usagedata-note.txt
 
-**Example - Loading a simulation from a CDL file..**
+**Example - Loading a simulation from a CDL file.**
 
 .. testcode:: ExLoadCDLFile
 
-    out_ws_group = MolDyn(Filename='NaF_DISF.cdl', Functions=['Fqt-total', 'Sqw-total'])
+    out_ws_group = MolDyn(Filename='NaF_DISF.cdl',
+                          Functions=['Fqt-total', 'Sqw-total'])
 
     for ws_name in out_ws_group.getNames():
       print ws_name