diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index a8f4806709628440274a87670fff8406a8a71a14..bc7b88c1c146670cafd7799a5366e901ec8e4a79 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,6 +2,18 @@ This document contains the ``fslpy`` release history in reverse chronological
 order.
 
 
+3.20.0 (Wednesday 10th July 2024)
+---------------------------------
+
+
+Added
+^^^^^
+
+
+* New functions/methods in the :mod:`.featanalysis` module and
+  :class:`.FEATImage` class for accessing F-test results (!460).
+
+
 3.19.1 (Wednesday 26th June 2024)
 ---------------------------------
 
diff --git a/fsl/data/featanalysis.py b/fsl/data/featanalysis.py
index 3e528ba233d750af70bb52be68fab514df67121d..94c776894c134aa53380aeb482d5815ce4e38002 100644
--- a/fsl/data/featanalysis.py
+++ b/fsl/data/featanalysis.py
@@ -22,10 +22,12 @@ following functions are provided:
    isFirstLevelAnalysis
    loadDesign
    loadContrasts
+   loadFTests
    loadFsf
    loadSettings
    getThresholds
    loadClusterResults
+   loadFEATDesignFile
 
 
 The following functions return the names of various files of interest:
@@ -39,11 +41,14 @@ The following functions return the names of various files of interest:
    getPEFile
    getCOPEFile
    getZStatFile
+   getZFStatFile
    getClusterMaskFile
+   getFClusterMaskFile
 """
 
 
 import                   collections
+import                   io
 import                   logging
 import os.path        as op
 import numpy          as np
@@ -166,55 +171,69 @@ def loadContrasts(featdir):
     :arg featdir: A FEAT directory.
     """
 
-    matrix       = None
-    numContrasts = 0
-    names        = {}
-    designcon    = op.join(featdir, 'design.con')
+    filename = op.join(featdir, 'design.con')
 
-    log.debug('Loading FEAT contrasts from {}'.format(designcon))
+    log.debug('Loading FEAT contrasts from %s', filename)
 
-    with open(designcon, 'rt') as f:
+    try:
+        designcon    = loadFEATDesignFile(filename)
+        contrasts    = np.genfromtxt(io.StringIO(designcon['Matrix']), ndmin=2)
+        numContrasts = int(designcon['NumContrasts'])
+        names        = []
 
-        while True:
-            line = f.readline().strip()
+        if numContrasts != contrasts.shape[0]:
+            raise RuntimeError(f'Matrix shape {contrasts.shape} does not '
+                               f'match number of contrasts {numContrasts}')
 
-            if line.startswith('/ContrastName'):
-                tkns       = line.split(None, 1)
-                num        = [c for c in tkns[0] if c.isdigit()]
-                num        = int(''.join(num))
+        contrasts = [list(row) for row in contrasts]
 
-                # The /ContrastName field may not
-                # actually have a name specified
-                if len(tkns) > 1:
-                    name       = tkns[1].strip()
-                    names[num] = name
+        for i in range(numContrasts):
+            cname = designcon.get(f'ContrastName{i + 1}', '')
+            if cname == '':
+                cname = f'{i + 1}'
+            names.append(cname)
 
-            elif line.startswith('/NumContrasts'):
-                numContrasts = int(line.split()[1])
+    except Exception as e:
+        log.debug('Error reading %s: %s', filename, e, exc_info=True)
+        raise RuntimeError(f'{filename} does not appear '
+                           'to be a valid design.con file') from e
 
-            elif line == '/Matrix':
-                break
+    return names, contrasts
 
-        matrix = np.loadtxt(f, ndmin=2)
 
-    if matrix       is None             or \
-       numContrasts != matrix.shape[0]:
-        raise RuntimeError('{} does not appear to be a '
-                           'valid design.con file'.format(designcon))
+def loadFTests(featdir):
+    """Loads F-tests from a FEAT directory. Returns a list of f-test vectors
+    (each of which is a list itself), where each vector contains a 1 or a 0
+    denoting the contrasts included in the F-test.
 
-    # Fill in any missing contrast names
-    if len(names) != numContrasts:
-        for i in range(numContrasts):
-            if i + 1 not in names:
-                names[i + 1] = str(i + 1)
+    :arg featdir: A FEAT directory.
+    """
 
-    names     = [names[c + 1] for c in range(numContrasts)]
-    contrasts = []
+    filename = op.join(featdir, 'design.fts')
 
-    for row in matrix:
-        contrasts.append(list(row))
+    if not op.exists(filename):
+        return []
 
-    return names, contrasts
+    log.debug('Loading FEAT F-tests from %s', filename)
+
+    try:
+        desfts = loadFEATDesignFile(filename)
+        ftests = np.genfromtxt(io.StringIO(desfts['Matrix']), ndmin=2)
+        ncols  = int(desfts['NumWaves'])
+        nrows  = int(desfts['NumContrasts'])
+
+        if ftests.shape != (nrows, ncols):
+            raise RuntimeError(f'Matrix shape {ftests.shape} does not match '
+                               f'number of EVs/FTests ({ncols}, {nrows})')
+
+        ftests = [list(row) for row in ftests]
+
+    except Exception as e:
+        log.debug('Error reading %s: %s', filename, e, exc_info=True)
+        raise RuntimeError(f'{filename} does not appear '
+                           'to be a valid design.fts file') from e
+
+    return ftests
 
 
 def loadFsf(designfsf):
@@ -228,7 +247,7 @@ def loadFsf(designfsf):
 
     settings  = collections.OrderedDict()
 
-    log.debug('Loading FEAT settings from {}'.format(designfsf))
+    log.debug('Loading FEAT settings from %s', designfsf)
 
     with open(designfsf, 'rt') as f:
 
@@ -310,19 +329,22 @@ def isFirstLevelAnalysis(settings):
     return settings['level'] == '1'
 
 
-def loadClusterResults(featdir, settings, contrast):
+def loadClusterResults(featdir, settings, contrast, ftest=False):
     """If cluster thresholding was used in the FEAT analysis, this function
     will load and return the cluster results for the specified (0-indexed)
-    contrast number.
+    contrast or f-test.
 
-    If there are no cluster results for the given contrast, ``None`` is
-    returned.
+    If there are no cluster results for the given contrast/f-test, ``None``
+    is returned.
 
     An error will be raised if the cluster file cannot be parsed.
 
     :arg featdir:  A FEAT directory.
     :arg settings: A FEAT settings dictionary.
-    :arg contrast: 0-indexed contrast number.
+    :arg contrast: 0-indexed contrast or f-test number.
+    :arg ftest:    If ``False`` (default), return cluster results for
+                   the contrast numbered ``contrast``. Otherwise, return
+                   cluster results for the f-test numbered ``contrast``.
 
     :returns:      A list of ``Cluster`` instances, each of which contains
                    information about one cluster. A ``Cluster`` object has the
@@ -343,11 +365,16 @@ def loadClusterResults(featdir, settings, contrast):
                                   gravity.
                      ``zcogz``    Z voxel coordinate of cluster centre of
                                   gravity.
-                     ``copemax``  Maximum COPE value in cluster.
-                     ``copemaxx`` X voxel coordinate of maximum COPE value.
+                     ``copemax``  Maximum COPE value in cluster (not
+                                  present for f-tests).
+                     ``copemaxx`` X voxel coordinate of maximum COPE value
+                                  (not present for f-tests).
                      ``copemaxy`` Y voxel coordinate of maximum COPE value.
+                                  (not present for f-tests).
                      ``copemaxz`` Z voxel coordinate of maximum COPE value.
+                                  (not present for f-tests).
                      ``copemean`` Mean COPE of all voxels in the cluster.
+                                  (not present for f-tests).
                      ============ =========================================
     """
 
@@ -357,8 +384,11 @@ def loadClusterResults(featdir, settings, contrast):
     # the ZMax/COG etc coordinates
     # are usually in voxel coordinates
     coordXform  = np.eye(4)
-    clusterFile = op.join(
-        featdir, 'cluster_zstat{}.txt'.format(contrast + 1))
+
+    if ftest: prefix = 'cluster_zfstat'
+    else:     prefix = 'cluster_zstat'
+
+    clusterFile = op.join(featdir, f'{prefix}{contrast + 1}.txt')
 
     if not op.exists(clusterFile):
 
@@ -367,8 +397,7 @@ def loadClusterResults(featdir, settings, contrast):
         # the cluster file will instead be called
         # 'cluster_zstatX_std.txt', so we'd better
         # check for that too.
-        clusterFile = op.join(
-            featdir, 'cluster_zstat{}_std.txt'.format(contrast + 1))
+        clusterFile = op.join(featdir, f'{prefix}{contrast + 1}_std.txt')
 
         if not op.exists(clusterFile):
             return None
@@ -379,9 +408,6 @@ def loadClusterResults(featdir, settings, contrast):
         # later on.
         coordXform = fslimage.Image(getDataFile(featdir)).worldToVoxMat
 
-    log.debug('Loading cluster results for contrast {} from {}'.format(
-        contrast, clusterFile))
-
     # The cluster.txt file is converted
     # into a list of Cluster objects,
     # each of which encapsulates
@@ -398,10 +424,18 @@ def loadClusterResults(featdir, settings, contrast):
 
             # if cluster thresholding was not used,
             # the cluster table will not contain
-            # P valuse.
+            # P values.
             if not hasattr(self, 'p'):    self.p    = 1.0
             if not hasattr(self, 'logp'): self.logp = 0.0
 
+            # F-test cluster results will not have
+            # COPE-* results
+            if not hasattr(self, 'copemax'):  self.copemax  = np.nan
+            if not hasattr(self, 'copemaxx'): self.copemaxx = np.nan
+            if not hasattr(self, 'copemaxy'): self.copemaxy = np.nan
+            if not hasattr(self, 'copemaxz'): self.copemaxz = np.nan
+            if not hasattr(self, 'copemean'): self.copemean = np.nan
+
     # This dict provides a mapping between
     # Cluster object attribute names, and
     # the corresponding column name in the
@@ -433,10 +467,9 @@ def loadClusterResults(featdir, settings, contrast):
         'COPE-MAX Z (mm)'  : 'copemaxz',
         'COPE-MEAN'        : 'copemean'}
 
-    # An error will be raised if the
-    # cluster file does not exist (e.g.
-    # if the specified contrast index
-    # is invalid)
+    log.debug('Loading cluster results for contrast %s from %s',
+              contrast, clusterFile)
+
     with open(clusterFile, 'rt') as f:
 
         # Get every line in the file,
@@ -458,12 +491,11 @@ def loadClusterResults(featdir, settings, contrast):
     colNames     =  colNames.split('\t')
     clusterLines = [cl      .split('\t') for cl in clusterLines]
 
-    # Turn each cluster line into a
-    # Cluster instance. An error will
-    # be raised if the columm names
-    # are unrecognised (i.e. not in
-    # the colmap above), or if the
-    # file is poorly formed.
+    # Turn each cluster line into a Cluster
+    # instance. An error will be raised if the
+    # columm names are unrecognised (i.e. not
+    # in the colmap above), or if the file is
+    # poorly formed.
     clusters = [Cluster(**dict(zip(colNames, cl))) for cl in clusterLines]
 
     # Make sure all coordinates are in voxels -
@@ -489,6 +521,35 @@ def loadClusterResults(featdir, settings, contrast):
     return clusters
 
 
+def loadFEATDesignFile(filename):
+    """Load a FEAT design file, e.g. ``design.mat``, ``design.con``, ``design.fts``.
+
+    These files contain key-value pairs, and are formatted according to an
+    undocumented structure where each key is of the form "/KeyName", and is
+    followed immediately by a whitespace character, and then the value.
+
+    :arg filename: File to load
+    :returns:      A dictionary of key-value pairs. The values are all left
+                   as strings.
+    """
+
+    fields = {}
+
+    with open(filename, 'rt') as f:
+        content = f.read()
+
+    content = content.split('/')
+    for line in content:
+        line = line.strip()
+        if line == '':
+            continue
+
+        name, value  = line.split(maxsplit=1)
+        fields[name] = value
+
+    return fields
+
+
 def getDataFile(featdir):
     """Returns the name of the file in the FEAT directory which contains
     the model input data (typically called ``filtered_func_data.nii.gz``).
@@ -532,7 +593,7 @@ def getPEFile(featdir, ev):
     :arg featdir: A FEAT directory.
     :arg ev:      The EV number (0-indexed).
     """
-    pefile = op.join(featdir, 'stats', 'pe{}'.format(ev + 1))
+    pefile = op.join(featdir, 'stats', f'pe{ev + 1}')
     return fslimage.addExt(pefile, mustExist=True)
 
 
@@ -544,7 +605,7 @@ def getCOPEFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    copefile = op.join(featdir, 'stats', 'cope{}'.format(contrast + 1))
+    copefile = op.join(featdir, 'stats', f'cope{contrast + 1}')
     return fslimage.addExt(copefile, mustExist=True)
 
 
@@ -556,10 +617,22 @@ def getZStatFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    zfile = op.join(featdir, 'stats', 'zstat{}'.format(contrast + 1))
+    zfile = op.join(featdir, 'stats', f'zstat{contrast + 1}')
     return fslimage.addExt(zfile, mustExist=True)
 
 
+def getZFStatFile(featdir, ftest):
+    """Returns the path of the Z-statistic file for the specified F-test.
+
+    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
+
+    :arg featdir: A FEAT directory.
+    :arg ftest:   The F-test number (0-indexed).
+    """
+    zffile = op.join(featdir, 'stats', f'zfstat{ftest + 1}')
+    return fslimage.addExt(zffile, mustExist=True)
+
+
 def getClusterMaskFile(featdir, contrast):
     """Returns the path of the cluster mask file for the specified contrast.
 
@@ -568,5 +641,17 @@ def getClusterMaskFile(featdir, contrast):
     :arg featdir:  A FEAT directory.
     :arg contrast: The contrast number (0-indexed).
     """
-    mfile = op.join(featdir, 'cluster_mask_zstat{}'.format(contrast + 1))
+    mfile = op.join(featdir, f'cluster_mask_zstat{contrast + 1}')
+    return fslimage.addExt(mfile, mustExist=True)
+
+
+def getFClusterMaskFile(featdir, ftest):
+    """Returns the path of the cluster mask file for the specified f-test.
+
+    Raises a :exc:`~fsl.utils.path.PathError` if the file does not exist.
+
+    :arg featdir:  A FEAT directory.
+    :arg contrast: The f-test number (0-indexed).
+    """
+    mfile = op.join(featdir, f'cluster_mask_zfstat{ftest + 1}')
     return fslimage.addExt(mfile, mustExist=True)
diff --git a/fsl/data/featdesign.py b/fsl/data/featdesign.py
index ce5b41f592609b5dbd6bc8b987e04016567a761c..c60b31871d154f3e0ef5a2e7312bdd0618b4d7f2 100644
--- a/fsl/data/featdesign.py
+++ b/fsl/data/featdesign.py
@@ -160,7 +160,7 @@ class FEATFSFDesign(object):
         # Print a warning if we're
         # using an old version of FEAT
         if version < 6:
-            log.warning('Unsupported FEAT version: {}'.format(version))
+            log.warning('Unsupported FEAT version: %s', version)
 
         # We need to parse the EVS a bit
         # differently depending on whether
@@ -210,8 +210,7 @@ class FEATFSFDesign(object):
                 continue
 
             if (not self.__loadVoxEVs) or (ev.filename is None):
-                log.warning('Voxel EV image missing '
-                            'for ev {}'.format(ev.index))
+                log.warning('Voxel EV image missing for ev %s', ev.index)
                 continue
 
             design[:, ev.index] = ev.getData(x, y, z)
@@ -250,8 +249,7 @@ class VoxelwiseEVMixin(object):
         if op.exists(filename):
             self.__filename = filename
         else:
-            log.warning('Voxelwise EV file does not '
-                        'exist: {}'.format(filename))
+            log.warning('Voxelwise EV file does not exist: %s', filename)
             self.__filename = None
 
         self.__image = None
@@ -502,11 +500,11 @@ def getFirstLevelEVs(featDir, settings, designMat):
     #   - voxelwise EVs
     for origIdx in range(origEVs):
 
-        title    = settings[        'evtitle{}'  .format(origIdx + 1)]
-        shape    = int(settings[    'shape{}'    .format(origIdx + 1)])
-        convolve = int(settings[    'convolve{}' .format(origIdx + 1)])
-        deriv    = int(settings[    'deriv_yn{}' .format(origIdx + 1)])
-        basis    = int(settings.get('basisfnum{}'.format(origIdx + 1), -1))
+        title    = settings[        f'evtitle{origIdx + 1}']
+        shape    = int(settings[    f'shape{origIdx + 1}'])
+        convolve = int(settings[    f'convolve{origIdx + 1}'])
+        deriv    = int(settings[    f'deriv_yn{origIdx + 1}'])
+        basis    = int(settings.get(f'basisfnum{origIdx + 1}', -1))
 
         # Normal EV. This is just a column
         # in the design matrix, defined by
@@ -525,8 +523,7 @@ def getFirstLevelEVs(featDir, settings, designMat):
 
             # The addExt function will raise an
             # error if the file does not exist.
-            filename = op.join(
-                featDir, 'designVoxelwiseEV{}'.format(origIdx + 1))
+            filename = op.join(featDir, f'designVoxelwiseEV{origIdx + 1}')
             filename = fslimage.addExt(filename, True)
 
             evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
@@ -607,7 +604,7 @@ def getFirstLevelEVs(featDir, settings, designMat):
         startIdx = len(evs) + 1
         if voxConfLocs != list(range(startIdx, startIdx + len(voxConfFiles))):
             raise FSFError('Unsupported voxelwise confound ordering '
-                           '({} -> {})'.format(startIdx, voxConfLocs))
+                           f'({startIdx} -> {voxConfLocs})')
 
         # Create the voxelwise confound EVs.
         # We make a name for the EV from the
@@ -680,7 +677,7 @@ def getHigherLevelEVs(featDir, settings, designMat):
     for origIdx in range(origEVs):
 
         # All we need is the title
-        title = settings['evtitle{}'.format(origIdx + 1)]
+        title = settings[f'evtitle{origIdx + 1}']
         evs.append(NormalEV(len(evs), origIdx, title))
 
     # Only the input file is specified for
@@ -689,7 +686,7 @@ def getHigherLevelEVs(featDir, settings, designMat):
     # name.
     for origIdx in range(voxEVs):
 
-        filename = settings['evs_vox_{}'.format(origIdx + 1)]
+        filename = settings[f'evs_vox_{origIdx + 1}']
         title    = op.basename(fslimage.removeExt(filename))
         evs.append(VoxelwiseEV(len(evs), origIdx, title, filename))
 
@@ -705,12 +702,12 @@ def loadDesignMat(designmat):
     :arg designmat: Path to the ``design.mat`` file.
     """
 
-    log.debug('Loading FEAT design matrix from {}'.format(designmat))
+    log.debug('Loading FEAT design matrix from %s', designmat)
 
     matrix = np.loadtxt(designmat, comments='/', ndmin=2)
 
     if matrix is None or matrix.size == 0 or len(matrix.shape) != 2:
-        raise FSFError('{} does not appear to be a '
-                       'valid design.mat file'.format(designmat))
+        raise FSFError(f'{designmat} does not appear '
+                       'to be a valid design.mat file')
 
     return matrix
diff --git a/fsl/data/featimage.py b/fsl/data/featimage.py
index 80aad19893ae4c504d2a0828e003c67134588fb0..430c71829501e5ea606bf6546a0cb6dffe89b294 100644
--- a/fsl/data/featimage.py
+++ b/fsl/data/featimage.py
@@ -63,8 +63,8 @@ class FEATImage(fslimage.Image):
             path = op.join(path, 'filtered_func_data')
 
         if not featanalysis.isFEATImage(path):
-            raise ValueError('{} does not appear to be data '
-                             'from a FEAT analysis'.format(path))
+            raise ValueError(f'{path} does not appear to be '
+                             'data from a FEAT analysis')
 
         featDir     = op.dirname(path)
         settings    = featanalysis.loadSettings( featDir)
@@ -72,9 +72,11 @@ class FEATImage(fslimage.Image):
         if featanalysis.hasStats(featDir):
             design      = featanalysis.loadDesign(   featDir, settings)
             names, cons = featanalysis.loadContrasts(featDir)
+            ftests      = featanalysis.loadFTests(   featDir)
         else:
             design      = None
             names, cons = [], []
+            ftests      = []
 
         fslimage.Image.__init__(self, path, **kwargs)
 
@@ -83,26 +85,31 @@ class FEATImage(fslimage.Image):
         self.__design        = design
         self.__contrastNames = names
         self.__contrasts     = cons
+        self.__ftests        = ftests
         self.__settings      = settings
 
         self.__residuals     =  None
         self.__pes           = [None] * self.numEVs()
         self.__copes         = [None] * self.numContrasts()
         self.__zstats        = [None] * self.numContrasts()
+        self.__zfstats       = [None] * self.numFTests()
         self.__clustMasks    = [None] * self.numContrasts()
+        self.__fclustMasks   = [None] * self.numFTests()
 
         if 'name' not in kwargs:
-            self.name = '{}: {}'.format(self.__analysisName, self.name)
+            self.name = f'{self.__analysisName}: {self.name}'
 
 
     def __del__(self):
         """Clears references to any loaded images."""
-        self.__design     = None
-        self.__residuals  = None
-        self.__pes        = None
-        self.__copes      = None
-        self.__zstats     = None
-        self.__clustMasks = None
+        self.__design      = None
+        self.__residuals   = None
+        self.__pes         = None
+        self.__copes       = None
+        self.__zstats      = None
+        self.__zfstats     = None
+        self.__clustMasks  = None
+        self.__fclustMasks = None
 
 
     def getFEATDir(self):
@@ -191,6 +198,11 @@ class FEATImage(fslimage.Image):
         return len(self.__contrasts)
 
 
+    def numFTests(self):
+        """Returns the number of f-tests in the analysis."""
+        return len(self.__ftests)
+
+
     def contrastNames(self):
         """Returns a list containing the name of each contrast in the analysis.
         """
@@ -206,6 +218,15 @@ class FEATImage(fslimage.Image):
         return [list(c) for c in self.__contrasts]
 
 
+    def ftests(self):
+        """Returns a list containing the analysis f-test vectors.
+
+        See :func:`.featanalysis.loadFTests`
+
+        """
+        return [list(f) for f in self.__ftests]
+
+
     def thresholds(self):
         """Returns the statistical thresholds used in the analysis.
 
@@ -214,14 +235,16 @@ class FEATImage(fslimage.Image):
         return featanalysis.getThresholds(self.__settings)
 
 
-    def clusterResults(self, contrast):
-        """Returns the clusters found in the analysis.
+    def clusterResults(self, contrast, ftest=False):
+        """Returns the clusters found in the analysis for the specified
+        contrast or f-test.
 
         See :func:.featanalysis.loadClusterResults`
         """
         return featanalysis.loadClusterResults(self.__featDir,
                                                self.__settings,
-                                               contrast)
+                                               contrast,
+                                               ftest)
 
 
     def getPE(self, ev):
@@ -229,12 +252,10 @@ class FEATImage(fslimage.Image):
 
         if self.__pes[ev] is None:
             pefile = featanalysis.getPEFile(self.__featDir, ev)
+            evname = self.evNames()[ev]
             self.__pes[ev] = fslimage.Image(
                 pefile,
-                name='{}: PE{} ({})'.format(
-                    self.__analysisName,
-                    ev + 1,
-                    self.evNames()[ev]))
+                name=f'{self.__analysisName}: PE{ev + 1} ({evname})')
 
         return self.__pes[ev]
 
@@ -246,7 +267,7 @@ class FEATImage(fslimage.Image):
             resfile = featanalysis.getResidualFile(self.__featDir)
             self.__residuals = fslimage.Image(
                 resfile,
-                name='{}: residuals'.format(self.__analysisName))
+                name=f'{self.__analysisName}: residuals')
 
         return self.__residuals
 
@@ -256,12 +277,10 @@ class FEATImage(fslimage.Image):
 
         if self.__copes[con] is None:
             copefile = featanalysis.getCOPEFile(self.__featDir, con)
+            conname  = self.contrastNames()[con]
             self.__copes[con] = fslimage.Image(
                 copefile,
-                name='{}: COPE{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
+                name=f'{self.__analysisName}: COPE{con + 1} ({conname})')
         return self.__copes[con]
 
 
@@ -270,35 +289,54 @@ class FEATImage(fslimage.Image):
         """
 
         if self.__zstats[con] is None:
-            zfile = featanalysis.getZStatFile(self.__featDir, con)
-
+            zfile   = featanalysis.getZStatFile(self.__featDir, con)
+            conname = self.contrastNames()[con]
             self.__zstats[con] = fslimage.Image(
                 zfile,
-                name='{}: zstat{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
-
+                name=f'{self.__analysisName}: zstat{con + 1} ({conname})')
         return self.__zstats[con]
 
 
+    def getZFStats(self, ftest):
+        """Returns the Z statistic image for the given f-test (0-indexed). """
+
+        if self.__zfstats[ftest] is None:
+            zfile = featanalysis.getZFStatFile(self.__featDir, ftest)
+            self.__zfstats[ftest] = fslimage.Image(
+                zfile,
+                name=f'{self.__analysisName}: zfstat{ftest + 1}')
+        return self.__zfstats[ftest]
+
+
     def getClusterMask(self, con):
         """Returns the cluster mask image for the given contrast (0-indexed).
         """
 
         if self.__clustMasks[con] is None:
-            mfile = featanalysis.getClusterMaskFile(self.__featDir, con)
-
+            mfile   = featanalysis.getClusterMaskFile(self.__featDir, con)
+            conname = self.contrastNames()[con]
             self.__clustMasks[con] = fslimage.Image(
                 mfile,
-                name='{}: cluster mask for zstat{} ({})'.format(
-                    self.__analysisName,
-                    con + 1,
-                    self.contrastNames()[con]))
+                name=f'{self.__analysisName}: cluster mask '
+                     f'for zstat{con + 1} ({conname})')
 
         return self.__clustMasks[con]
 
 
+    def getFClusterMask(self, ftest):
+        """Returns the cluster mask image for the given f-test (0-indexed).
+        """
+
+        if self.__fclustMasks[ftest] is None:
+            mfile = featanalysis.getFClusterMaskFile(self.__featDir, ftest)
+            self.__fclustMasks[ftest] = fslimage.Image(
+                mfile,
+                name=f'{self.__analysisName}: cluster mask '
+                     f'for zfstat{ftest + 1}')
+
+        return self.__fclustMasks[ftest]
+
+
     def fit(self, contrast, xyz):
         """Calculates the model fit for the given contrast vector
         at the given voxel. See the :func:`modelFit` function.
diff --git a/fsl/tests/__init__.py b/fsl/tests/__init__.py
index c1b105dcd45052efeed7ebe6ba5afd0d1f188e93..a758db12e28746a44d9324a0b8a190590d87239c 100644
--- a/fsl/tests/__init__.py
+++ b/fsl/tests/__init__.py
@@ -288,7 +288,8 @@ def make_mock_feat_analysis(featdir,
                             copes=True,
                             zstats=True,
                             residuals=True,
-                            clustMasks=True):
+                            clustMasks=True,
+                            zfstats=True):
 
     if xform is None:
         xform = np.eye(4)
@@ -346,6 +347,11 @@ def make_mock_feat_analysis(featdir,
         otherFiles .extend(files)
         otherShapes.extend([shape] * len(files))
 
+    if zfstats:
+        files = glob.glob(op.join(featdir, 'stats', 'zfstat*nii.gz'))
+        otherFiles .extend(files)
+        otherShapes.extend([shape] * len(files))
+
     if residuals:
         files = glob.glob(op.join(featdir, 'stats', 'res4d.nii.gz'))
         otherFiles .extend(files)
diff --git a/fsl/tests/test_featanalysis.py b/fsl/tests/test_featanalysis.py
index 6f07109ed18b5c20f6cd810988a20deb97f26f87..b35b48f8c846d98b22988d75d582614ef5e793db 100644
--- a/fsl/tests/test_featanalysis.py
+++ b/fsl/tests/test_featanalysis.py
@@ -210,6 +210,77 @@ def test_loadContrasts():
                 featanalysis.loadContrasts(featdir)
 
 
+def test_loadFTests():
+
+    goodtests = [
+        ("""
+         /NumWaves 4
+         /NumContrasts 3
+         /Matrix
+         0 1 0 1
+         0 0 1 1
+         1 1 1 1
+         """,
+         [[0, 1, 0, 1],
+          [0, 0, 1, 1],
+          [1, 1, 1, 1]]),
+        ("""
+         /NumWaves 10
+         /NumContrasts 2
+         /Matrix
+         0 1 0 1 0 1 1 0 0 1
+         0 0 1 1 1 0 0 1 0 0
+         """,
+         [[0, 1, 0, 1, 0, 1, 1, 0, 0, 1],
+          [0, 0, 1, 1, 1, 0, 0, 1, 0, 0]]),
+    ]
+    badtests  = [
+        """
+        /NumWaves 10
+        /NumContrasts 2
+        """,
+        """
+        /NumContrasts 2
+        /Matrix
+        1 0
+        0 1
+        """,
+        """
+        /NumWaves Burgers
+        /NumContrasts 2
+        /Matrix
+        1 0
+        0 1
+        """,
+        """
+        /Matrix
+        1 0
+        0 1
+        """,
+        """
+        /NumWaves  4
+        /NumContrasts 3
+        /Matrix
+        1 0 0 0 1 0 0
+        0 1 0 0 1 0 0
+        """,
+    ]
+
+    with tests.testdir() as testdir:
+        featdir = op.join(testdir, 'analysis.feat')
+        for contents, expect in goodtests:
+            designcon = op.join(featdir, 'design.fts')
+            tests.make_dummy_file(designcon, textwrap.dedent(contents).strip())
+            assert featanalysis.loadFTests(featdir) == expect
+
+        for contents in badtests:
+            designcon = op.join(featdir, 'design.fts')
+            tests.make_dummy_file(designcon, textwrap.dedent(contents).strip())
+            with pytest.raises(Exception):
+                featanalysis.loadFTests(featdir)
+
+
+
 def test_loadSettings():
 
     contents = """
@@ -289,7 +360,10 @@ def test_loadClusterResults():
                 '2ndlevel_1.gfeat/cope1.feat', '2ndlevel_1.gfeat/cope2.feat',
                 '2ndlevel_2.gfeat/cope1.feat', '2ndlevel_2.gfeat/cope2.feat']
     ncontrasts = [2, 2, 2, 1, 1, 1, 1]
-    nclusters  = [[1, 5], [2, 2], [3, 5], [7], [1], [10], [27]]
+    nftests    = [0, 0, 1, 0, 0, 1, 1]
+
+    # nclusters = [contrastclusters] + [ftestclusters]
+    nclusters  = [[1, 5], [2, 2], [3, 5, 3], [7], [1], [10, 8], [27, 21]]
 
     with pytest.raises(Exception):
         featanalysis.loadClusterResults('notafeatdir')
@@ -320,14 +394,20 @@ def test_loadClusterResults():
                 fslimage.Image(data, xform=xform).save(datafile)
 
             settings = featanalysis.loadSettings(featdir)
+            # contrasts
             for c in range(ncontrasts[i]):
                 clusters = featanalysis.loadClusterResults(
                     featdir, settings, c)
-
                 assert len(clusters) == nclusters[i][c]
+            # f-tests
+            for c in range(nftests[i]):
+                clusters = featanalysis.loadClusterResults(
+                    featdir, settings, c, ftest=True)
+                assert len(clusters) == nclusters[i][c + ncontrasts[i]]
 
             # Test calling the function on a feat dir
             # which doesn't have any cluster results
+            # (2ndlevel_2.gfeat)
             if i == len(featdirs) - 1:
                 for clustfile in glob.glob(op.join(featdir, 'cluster*txt')):
                     os.remove(clustfile)
@@ -498,6 +578,53 @@ def test_getZStatFile():
                         featanalysis.getZStatFile(featdir, zi)
 
 
+def test_getZStatFile():
+    testcases = [
+        (['analysis.feat/stats/zstat1.nii.gz',
+          'analysis.feat/stats/zstat2.nii.gz'], True),
+        (['analysis.feat/stats/zstat1.nii.gz'], True),
+        (['analysis.feat/stats/zstat0.nii.gz'], False),
+        (['analysis.feat/stats/zstat1.txt'],    False),
+    ]
+
+    for paths, shouldPass in testcases:
+        with tests.testdir(paths) as testdir:
+            featdir = op.join(testdir, 'analysis.feat')
+
+            for zi in range(len(paths)):
+                expect = op.join(
+                    featdir, 'stats', 'zstat{}.nii.gz'.format(zi + 1))
+
+                if shouldPass:
+                    assert featanalysis.getZStatFile(featdir, zi) == expect
+                else:
+                    with pytest.raises(fslpath.PathError):
+                        featanalysis.getZStatFile(featdir, zi)
+
+
+def test_getZFStatFile():
+    testcases = [
+        (['analysis.feat/stats/zfstat1.nii.gz',
+          'analysis.feat/stats/zfstat2.nii.gz'], True),
+        (['analysis.feat/stats/zfstat1.nii.gz'], True),
+        (['analysis.feat/stats/zfstat0.nii.gz'], False),
+        (['analysis.feat/stats/zfstat1.txt'],    False),
+    ]
+    for paths, shouldPass in testcases:
+        with tests.testdir(paths) as testdir:
+            featdir = op.join(testdir, 'analysis.feat')
+
+            for zi in range(len(paths)):
+                expect = op.join(
+                    featdir, 'stats', 'zfstat{}.nii.gz'.format(zi + 1))
+
+                if shouldPass:
+                    assert featanalysis.getZFStatFile(featdir, zi) == expect
+                else:
+                    with pytest.raises(fslpath.PathError):
+                        featanalysis.getZFStatFile(featdir, zi)
+
+
 def test_getClusterMaskFile():
     testcases = [
         (['analysis.feat/cluster_mask_zstat1.nii.gz',
@@ -520,3 +647,27 @@ def test_getClusterMaskFile():
                 else:
                     with pytest.raises(fslpath.PathError):
                         featanalysis.getClusterMaskFile(featdir, ci)
+
+
+def test_getFClusterMaskFile():
+    testcases = [
+        (['analysis.feat/cluster_mask_zfstat1.nii.gz',
+          'analysis.feat/cluster_mask_zfstat2.nii.gz'], True),
+        (['analysis.feat/cluster_mask_zfstat1.nii.gz'], True),
+        (['analysis.feat/cluster_mask_zfstat0.nii.gz'], False),
+        (['analysis.feat/cluster_mask_zfstat1.txt'],    False),
+    ]
+
+    for paths, shouldPass in testcases:
+        with tests.testdir(paths) as testdir:
+            featdir = op.join(testdir, 'analysis.feat')
+
+            for ci in range(len(paths)):
+                expect = op.join(
+                    featdir, 'cluster_mask_zfstat{}.nii.gz'.format(ci + 1))
+
+                if shouldPass:
+                    assert featanalysis.getFClusterMaskFile(featdir, ci) == expect
+                else:
+                    with pytest.raises(fslpath.PathError):
+                        featanalysis.getFClusterMaskFile(featdir, ci)
diff --git a/fsl/tests/test_featimage.py b/fsl/tests/test_featimage.py
index 370e797366597d44afaf39414e81cdf7f185bd5e..05f5cf5e3f8e1517f28e0459723f0bc6a361cc01 100644
--- a/fsl/tests/test_featimage.py
+++ b/fsl/tests/test_featimage.py
@@ -88,7 +88,8 @@ def test_FEATImage_attributes():
                     copes=False,
                     zstats=False,
                     residuals=False,
-                    clustMasks=False)
+                    clustMasks=False,
+                    zfstats=False)
             else:
                 featdir = op.join(datadir, featdir)
 
@@ -100,6 +101,7 @@ def test_FEATImage_attributes():
             design   = featdesign.FEATFSFDesign(featdir, settings)
             desmat   = design.getDesign()
             evnames  = [ev.title for ev in design.getEVs()]
+            ftests   = featanalysis.loadFTests(featdir)
             contrastnames, contrasts = featanalysis.loadContrasts(featdir)
 
             assert np.all(np.isclose(fi.shape,         shape))
@@ -115,8 +117,10 @@ def test_FEATImage_attributes():
             assert fi.numEVs()                 == desmat.shape[1]
             assert fi.evNames()                == evnames
             assert fi.numContrasts()           == len(contrasts)
+            assert fi.numFTests()              == len(ftests)
             assert fi.contrastNames()          == contrastnames
             assert fi.contrasts()              == contrasts
+            assert fi.ftests()                 == ftests
             assert np.all(np.isclose(fi.getDesign(), desmat))
 
             assert fi.thresholds() == featanalysis.getThresholds(settings)
@@ -153,9 +157,10 @@ def test_FEATImage_imageAccessors():
             shape4D = shape
             shape   = shape4D[:3]
 
-            fi    = featimage.FEATImage(featdir)
-            nevs  = fi.numEVs()
-            ncons = fi.numContrasts()
+            fi      = featimage.FEATImage(featdir)
+            nevs    = fi.numEVs()
+            ncons   = fi.numContrasts()
+            nftests = fi.numFTests()
 
             # Testing the FEATImage internal cache
             for i in range(2):
@@ -166,6 +171,9 @@ def test_FEATImage_imageAccessors():
                     assert fi.getCOPE(       con).shape == shape
                     assert fi.getZStats(     con).shape == shape
                     assert fi.getClusterMask(con).shape == shape
+                for ft in range(nftests):
+                    assert fi.getZFStats(     ft).shape == shape
+                    assert fi.getFClusterMask(ft).shape == shape
             del fi
             fi = None
 
diff --git a/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_mask_zfstat1.nii.gz b/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_mask_zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f5bed2405c8ed254de724086c5c2ec03466297b9
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_mask_zfstat1.nii.gz
@@ -0,0 +1 @@
+./cluster_mask_zfstat1.nii.gz
diff --git a/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_zfstat1.txt b/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_zfstat1.txt
new file mode 100644
index 0000000000000000000000000000000000000000..26a4ad56d943492258ee25b80060bbf3f6967b91
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/1stlevel_3.feat/cluster_zfstat1.txt
@@ -0,0 +1,4 @@
+Cluster Index	Voxels	P	-log10(P)	Z-MAX	Z-MAX X (vox)	Z-MAX Y (vox)	Z-MAX Z (vox)	Z-COG X (vox)	Z-COG Y (vox)	Z-COG Z (vox)
+3	17	0.00497	2.3	3.47	34	9	2	34.7	9.59	1.96
+2	14	0.0163	1.79	3.44	39	11	3	39.7	11	2.27
+1	13	0.0247	1.61	3.7	29	9	3	29.2	9.01	2.58
diff --git a/fsl/tests/testdata/test_feat/1stlevel_3.feat/design.fts b/fsl/tests/testdata/test_feat/1stlevel_3.feat/design.fts
new file mode 100644
index 0000000000000000000000000000000000000000..b243edb2198719d2935327395415915245f49746
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/1stlevel_3.feat/design.fts
@@ -0,0 +1,5 @@
+/NumWaves	2
+/NumContrasts	1
+
+/Matrix
+1 1
diff --git a/fsl/tests/testdata/test_feat/1stlevel_3.feat/stats/zfstat1.nii.gz b/fsl/tests/testdata/test_feat/1stlevel_3.feat/stats/zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..c396817eca92baed3a21b62b3b78ee694e7d9157
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/1stlevel_3.feat/stats/zfstat1.nii.gz
@@ -0,0 +1 @@
+./stats/zfstat1.nii.gz
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_mask_zfstat1.nii.gz b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_mask_zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..e98358baa21c5a2565eba4051b472930ed854796
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_mask_zfstat1.nii.gz
@@ -0,0 +1 @@
+./cluster_mask_zstat1.nii.gz
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_zfstat1_std.txt b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_zfstat1_std.txt
new file mode 100644
index 0000000000000000000000000000000000000000..11679b5b741f745dc64eea54a2cec31d26d8c7a4
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/cluster_zfstat1_std.txt
@@ -0,0 +1,9 @@
+Cluster Index	Voxels	P	-log10(P)	Z-MAX	Z-MAX X (mm)	Z-MAX Y (mm)	Z-MAX Z (mm)	Z-COG X (mm)	Z-COG Y (mm)	Z-COG Z (mm)
+8	43	9.08e-09	8.04	3.06	-46	28	10	-47.7	22.6	5.11
+7	32	9.54e-07	6.02	3.5	54	32	0	53.4	29.7	1.21
+6	21	0.00017	3.77	3.93	16	-8	8	16.4	-5.41	2.97
+5	20	0.000285	3.55	4	54	24	8	50.3	21	7.62
+4	14	0.00763	2.12	3.37	-20	10	8	-21	15.3	9.1
+3	13	0.0137	1.86	2.85	-38	12	14	-40.8	15	11.8
+2	12	0.0251	1.6	3.63	42	-48	16	42.5	-47.3	12.4
+1	11	0.0462	1.34	3.11	-28	46	16	-25.7	43.3	14.9
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/design.fts b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/design.fts
new file mode 100644
index 0000000000000000000000000000000000000000..129044f64ed86575c86dd60038bc7661b71a706d
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/design.fts
@@ -0,0 +1,5 @@
+/NumWaves	1
+/NumContrasts	1
+
+/Matrix
+1
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/stats/zfstat1.nii.gz b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/stats/zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..c396817eca92baed3a21b62b3b78ee694e7d9157
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope1.feat/stats/zfstat1.nii.gz
@@ -0,0 +1 @@
+./stats/zfstat1.nii.gz
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_mask_zfstat1.nii.gz b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_mask_zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..f5bed2405c8ed254de724086c5c2ec03466297b9
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_mask_zfstat1.nii.gz
@@ -0,0 +1 @@
+./cluster_mask_zfstat1.nii.gz
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_zfstat1_std.txt b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_zfstat1_std.txt
new file mode 100644
index 0000000000000000000000000000000000000000..acc340e1893b0779fcd89a6134330fa0a83f59a5
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/cluster_zfstat1_std.txt
@@ -0,0 +1,22 @@
+Cluster Index	Voxels	P	-log10(P)	Z-MAX	Z-MAX X (mm)	Z-MAX Y (mm)	Z-MAX Z (mm)	Z-COG X (mm)	Z-COG Y (mm)	Z-COG Z (mm)
+21	41	9.48e-09	8.02	3.46	42	-58	0	41.2	-63.6	-0.463
+20	33	2.98e-07	6.53	3.5	32	38	-10	35.4	37.9	-7.22
+19	33	2.98e-07	6.53	3.65	-30	-74	0	-30.4	-73.3	0.808
+18	33	2.98e-07	6.53	3.59	18	-14	2	17.3	-16.4	1.79
+17	24	2.26e-05	4.64	2.93	60	-58	-8	61.1	-56.6	-7.76
+16	22	6.3e-05	4.2	3.06	40	-38	8	40.9	-41.4	6.31
+15	21	0.000106	3.97	3.28	22	-54	4	20.7	-50.5	1.9
+14	19	0.00031	3.51	3.5	54	-70	6	50.5	-73.3	1.96
+13	18	0.000537	3.27	3.92	30	-16	2	27.7	-20.7	3.76
+12	15	0.00297	2.53	4.27	-58	26	14	-60.6	26.8	15.3
+11	15	0.00297	2.53	3.01	-30	-44	4	-27.4	-45.3	5.42
+10	15	0.00297	2.53	3.05	-10	-42	-8	-7.07	-42.9	-5.51
+9	14	0.00539	2.27	3.2	34	-16	-4	37.7	-17.7	-3.36
+8	14	0.00539	2.27	3.4	4	-50	2	2.8	-48.1	0.972
+7	14	0.00539	2.27	3.44	20	26	-10	16.5	24.1	-9.28
+6	13	0.00989	2	3.88	30	-94	-2	26.3	-97.6	-2.54
+5	13	0.00989	2	3.26	-18	-38	-2	-17.3	-40.6	-1.04
+4	12	0.0184	1.73	3.31	-46	24	-6	-47.6	24.3	-7.44
+3	11	0.0348	1.46	4.37	-32	-46	-6	-31.7	-48.6	-6.41
+2	11	0.0348	1.46	3.03	2	-28	-4	-0.26	-28	-3.62
+1	11	0.0348	1.46	3.06	26	-56	-6	29.7	-53.9	-7.4
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/design.fts b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/design.fts
new file mode 100644
index 0000000000000000000000000000000000000000..129044f64ed86575c86dd60038bc7661b71a706d
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/design.fts
@@ -0,0 +1,5 @@
+/NumWaves	1
+/NumContrasts	1
+
+/Matrix
+1
diff --git a/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/stats/zfstat1.nii.gz b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/stats/zfstat1.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..b013ebacb3cd988aef40a9c41cda9aab568d58b5
--- /dev/null
+++ b/fsl/tests/testdata/test_feat/2ndlevel_2.gfeat/cope2.feat/stats/zfstat1.nii.gz
@@ -0,0 +1 @@
+./cope2.feat/stats/zfstat1.nii.gz
diff --git a/fsl/version.py b/fsl/version.py
index 4cad9f591ca222583b3f3587807ce95fbaf5798b..d40259f541448c9137101a2669c6adfd55cf14f9 100644
--- a/fsl/version.py
+++ b/fsl/version.py
@@ -47,7 +47,7 @@ import            re
 import            string
 
 
-__version__ = '3.20.0.dev0'
+__version__ = '3.21.0.dev0'
 """Current version number, as a string. """