diff --git a/xbout/__init__.py b/xbout/__init__.py
index 8cffd586..5eb7f21d 100644
--- a/xbout/__init__.py
+++ b/xbout/__init__.py
@@ -1,5 +1,8 @@
 from .load import open_boutdataset, collect
 
+from .options import OptionsFile
+evaluate = OptionsFile.evaluate
+
 from . import geometries
 from .geometries import register_geometry, REGISTERED_GEOMETRIES
 
diff --git a/xbout/boutdataset.py b/xbout/boutdataset.py
index bf267d00..b76daa41 100644
--- a/xbout/boutdataset.py
+++ b/xbout/boutdataset.py
@@ -14,6 +14,8 @@
 
 from .plotting.animate import animate_poloidal, animate_pcolormesh, animate_line
 from .plotting.utils import _create_norm
+
+
 @register_dataset_accessor('bout')
 class BoutDatasetAccessor:
     """
diff --git a/xbout/options.py b/xbout/options.py
new file mode 100644
index 00000000..bf80255f
--- /dev/null
+++ b/xbout/options.py
@@ -0,0 +1,364 @@
+from collections import UserDict
+from pathlib import Path
+
+
+# These are imported to be used by 'evaluate=True' in Section.get().
+# Change the names to match those used by C++/BOUT++
+from numpy import (pi, sin, cos, tan, arccos as acos, arcsin as asin,
+                   arctan as atan, arctan2 as atan2, sinh, cosh, tanh,
+                   arcsinh as asinh, arccosh as acosh, arctanh as atanh,
+                   exp, log, log10, power as pow, sqrt, ceil, floor,
+                   round, abs)
+
+# TODO file reader
+# TODO substitution of keys within strings
+# TODO Sanitise some of the weirder things BOUT's option parser does:
+#  - Detect escaped arithmetic symbols (+-*/^), dot (.), or brackets ((){}[]) in option names
+#  - Detect escaping through backquotes in option names (`2ndvalue`)
+#  - Evaluate numletter as num*letter (and numbracket as num*bracket)
+#  - Substitute unicode representation of pi
+#  - Ensure option names don't contain ':' or '='
+#  - Check if all expressions are rounded to nearest integer?!
+# TODO ability to read from/write to nested dictionary
+
+
+SECTION_DELIM = ':'
+COMMENT_DELIM = ['#', ';']
+INDENT_STR = '|-- '
+BOOLEAN_STATES = {'true': True, 'on': True,
+                  'false': False, 'off': False}
+
+
+class Section(UserDict):
+    """
+    Section of an options file.
+
+    This section can have multiple sub-sections, and each section stores the
+    options as multiple key-value pairs.
+
+    Parameters
+    ----------
+    name : str
+        Name of the section
+    data : dict, optional
+    parent : str, Section or None
+        The parent Section
+    """
+
+    # Inheriting from UserDict gives keys, items, values, iter, del, & len
+    # methods already, and the data attribute for storing the contents
+
+    def __init__(self, name, parent=None, data=None):
+        if data is None: data = {}
+        super().__init__(data)
+
+        self.name = name
+        # TODO check if parent has a section with the same name?
+        self._parent = parent
+        if isinstance(parent, Section):
+            self.parent = parent.name
+        else:
+            self.parent = parent
+
+    def __getitem__(self, key):
+        return self.get(key, evaluate=False, keep_comments=False)
+
+    def get(self, key, substitute=False, evaluate=False, keep_comments=False):
+        """
+        Fetch the value stored under a certain key.
+
+        Parameters
+        ----------
+        key : str
+        substitute : bool, optional (default: False)
+            If value contains other values, referenced as `section:key`, will
+            substitute them in to the value string.
+        evaluate : bool, optional (default: False)
+            If true, attempts to evaluate the value as an expression by
+            substituting in other values from the options file. Other values
+            are specified by key, or by section:key. Also sets substitute=True.
+            If false, will return value as an unmodified string.
+        keep_comments : bool, optional (default: False)
+            If false, will strip off any inline comments, delimited by '#', and
+            any whitespace before returning.
+            If true, will return the whole line.
+        """
+
+        if evaluate and keep_comments:
+            # TODO relax this?
+            raise ValueError("Cannot keep comments and evaluate the contents.")
+
+        line = self.data[key]
+        if keep_comments or isinstance(line, Section):
+            return line
+        else:
+            # any comment will be discarded, along with any trailing whitespace
+            for delim in COMMENT_DELIM:
+                line, *comments = line.split(delim)
+            value = line.rstrip()
+
+            if substitute or evaluate:
+                value = self._substitute_keys_within(value)
+
+            if evaluate:
+                return self.evaluate(value)
+            else:
+                return value
+
+    def __setitem__(self, key, value):
+        self.set(key, value)
+
+    def set(self, key, value):
+        """
+        Store a value under a certain key.
+
+        Will cast the value to a string (unless value is a new section).
+        If a dictionary is passed, will create a new Section.
+
+        Parameters
+        ----------
+        key : str
+        value : str
+        """
+
+        if isinstance(value, dict):
+            value = Section(name=key, parent=self, data=value)
+        if not isinstance(value, Section):
+            value = str(value)
+        self.data[key] = value
+
+    def lineage(self):
+        """
+        Returns the full route to this section by working back up the tree.
+
+        Returns
+        -------
+        str
+        """
+        if self._parent is None or self._parent == 'root':
+            return self.name
+        else:
+            return self._parent.lineage() + SECTION_DELIM + self.name
+
+    def _find_sections(self, sections):
+        """
+        Recursively find a list of all the section objects below this one, and
+        append them to the list passed before returning them all.
+        """
+        for key in self.keys():
+            val = self[key]
+            if isinstance(val, Section):
+                sections.append(val)
+                sections = val._find_sections(sections)
+        return sections
+
+    def sections(self):
+        """
+        Returns a list of all sections contained, including nested ones.
+
+        Returns
+        -------
+        list of Section objects
+        """
+        return self._find_sections([self])
+
+    def __str__(self):
+        depth = self.lineage().count(SECTION_DELIM)
+        text = INDENT_STR * depth + f"[{self.name}]\n"
+        for key, val in self.items():
+            if isinstance(val, Section):
+                text += str(val)
+            else:
+                text += INDENT_STR * (depth+1) + f"{key} = {val}\n"
+        return text
+
+    def _write(self, file, substitute=False, evaluate=False,
+               keep_comments=True):
+        """
+        Writes out a single section to an open file object as
+        [...:parent:section]
+        key = value     # comment
+        ...
+
+        If it encounters a nested section it will recursively call this method
+        on that nested section.
+        """
+
+        if self.name is not None and self.name != 'root':
+            section_header = f"[{self.lineage()}]\n"
+            file.write(section_header)
+
+        for key in self.keys():
+            entry = self.get(key, substitute, evaluate, keep_comments)
+
+            if isinstance(entry, Section):
+                # Recursively write sections
+                file.write("\n")
+                entry._write(file, substitute, evaluate, keep_comments)
+            else:
+                file.write(f"{key} = {entry}\n")
+
+    def __repr__(self):
+        return f"Section(name='{self.name}', parent='{self.parent}', " \
+               f"data={self.data})"
+
+    @staticmethod
+    def evaluate(value):
+        """
+        Evaluates the string using eval, following the conventions of BOUT++.
+
+        Parameters
+        ----------
+        value : str
+        """
+
+        if value.lower() in BOOLEAN_STATES:
+            # Treat booleans separately to cover lowercase 'true' etc.
+            return BOOLEAN_STATES[value.lower()]
+        else:
+            if '^' in value:
+                value = value.replace('^', '**')
+            return eval(value)
+
+    def _substitute_keys_within(self, value):
+        # Detect if any colons
+
+        # Parse candidate keys contained
+
+        # Loop over candidate keys
+            # Search for matches in whole tree
+            # If found, replace, if not, throw error
+
+        # If any keys left?
+            # Check none of them point to this value (i.e. no cycles)
+            # Call _substitute_keys again
+        # Otherwise return ready for evaluation
+        return value
+
+
+class OptionsTree(Section):
+    """
+    Tree of options, with the same structure as a complete BOUT++ input file.
+
+    This class represents a tree structure. Each section (Section object) can
+    have multiple sub-sections, and each section stores the options as multiple
+    key-value pairs.
+
+    Examples
+    --------
+
+    >>> opts = OptionsTree()  # Create a root
+
+    Specify value of a key in a section "test"
+
+    >>> opts["test"]["key"] = 4
+
+    Get the value of a key in a section "test"
+    If the section does not exist then a KeyError is raised
+
+    >>> print(opts["test"]["key"])
+    4
+
+    To pretty print the options
+
+    >>> print(opts)
+    root
+     |- test
+     |   |- key = 4
+
+    """
+
+    def __init__(self, data=None):
+        super().__init__(name='root', data=data, parent=None)
+
+    # TODO .as_dict() ?
+
+    def write_to(self, file, substitute=False, evaluate=False,
+                 keep_comments=True, lower=False):
+        """
+        Writes out contents to a file, following the format of a BOUT.inp file.
+
+        Parameters
+        ----------
+        file
+        evaluate : bool, optional (default: False)
+            If true, attempts to evaluate the value as an expression by
+            substituting in other values from the options file, before writing.
+            Other values are specified by key, or by section:key.
+            If false, will write value as an unmodified string.
+        keep_comments : bool, optional (default: False)
+            If false, will strip off any inline comments, delimited by '#', and
+            any whitespace before writing.
+            If true, will write out the whole line.
+
+        Returns
+        -------
+        filepath : str
+            Full path of output file written to
+        """
+
+        with open(Path(file), 'w') as f:
+            self._write(file=f, substitute=False, evaluate=False,
+                        keep_comments=True)
+        return str(Path(file).resolve())
+
+    def _read_from(self, filepath, lower):
+        with open(filepath, 'r') as f:
+            # TODO add in first section header?
+            #for l in f:
+            #    line = Line(l)
+            #    if line.is_section():
+            #
+            #        name
+            #        self.data[name] = Section(name, parent, data)
+            #    elif line.is_comment() or line.is_empty()
+            #        continue
+            #    else:
+            data = None
+        return data, str(filepath.resolve())
+
+
+class OptionsFile(OptionsTree):
+    """
+    The full contents of a particular BOUT++ input file.
+
+    This class represents a tree structure of options, all loaded from `file`.
+    Each section (Section object) can have multiple sub-sections, and each
+    section stores the options as multiple key-value pairs.
+
+    Parameters
+    ----------
+    file : str or path-like object
+        Path to file from which to read input options.
+    lower : bool, optional (default: False)
+        If true, converts all strings to lowercase on reading.
+    """
+
+    # TODO gridfile stuff?
+
+    def __init__(self, file, lower=False):
+        contents, self.file = self._read_from(Path(file), lower)
+        super().__init__(data=contents)
+
+    def write(self, substitute=False, evaluate=False, keep_comments=True):
+        """
+        Writes out contents to the file, following the BOUT.inp format.
+
+        Parameters
+        ----------
+        evaluate : bool, optional (default: False)
+            If true, attempts to evaluate the value as an expression by
+            substituting in other values from the options file, before writing.
+            Other values are specified by key, or by section:key.
+            If false, will write value as an unmodified string.
+        keep_comments : bool, optional (default: False)
+            If false, will strip off any inline comments, delimited by '#', and
+            any whitespace before writing.
+            If true, will write out the whole line.
+        """
+
+        self.write_to(self.file, substitute, evaluate, keep_comments)
+
+    def __repr__(self):
+        # TODO Add grid-related options
+        return f"OptionsFile(file='{self.file}')"
diff --git a/xbout/tests/data/options/BOUT.inp b/xbout/tests/data/options/BOUT.inp
index e69de29b..c5a0a84c 100644
--- a/xbout/tests/data/options/BOUT.inp
+++ b/xbout/tests/data/options/BOUT.inp
@@ -0,0 +1,143 @@
+timestep              = 80                 # Timestep length of outputted data
+nout                  = 2000                # Number of outputted timesteps
+
+#nxpe = 4
+
+MZ                    = 256                 # Number of Z points
+zmin                  = 0
+zmax                  = 23.873241           # 23.873241*2pi = 150, z is fracs of 2pi, so Lz = 150
+myg                   = 0                   # No need for Y communications
+
+[mesh:ddx]
+first                 = C2
+second                = C2
+upwind                = C2
+
+[mesh:ddz]
+first                 = C2
+second                = C2
+upwind                = C2
+
+[mesh]
+Ly                    = 2000.0             # ~4m
+Lx                    = 400.0
+
+nx                    = 484                # including 4 guard cells
+ny                    = 1                  # excluding guard cells
+dx                    = Lx/(nx-4)
+dy                    = Ly/ny 
+
+[laplace]
+type                  = cyclic
+global_flags          = 0 
+inner_boundary_flags  = 1
+outer_boundary_flags  = 16
+include_yguards       = false
+
+[solver]
+type=cvode
+#timestep = 0.00000001    # Suggested init timestep for numerical methods
+mxstep = 100000000       # max steps before result is deemed not to converge
+#ATOL =   1e-12          # Absolute tolerance
+#RTOL =   1e-5           # Relative tolerance
+
+[storm]
+B_0                   = 0.24                # Tesla
+T_e0                  = 15                  # eV
+T_i0                  = 30                  # eV
+m_i                   = 2                   # Atomic Units
+q                     = 6.2                 # Dimensionless
+R_c                   = 1.5                 # m
+n_0                   = 0.5e13              # m-3
+Z                     = 1                   # Dimensionless
+loglambda             = -1                  # Dimensionless
+
+# If these parameters are specified, they will be used instead of the values calculated from the primary parameters above.
+# Using ESEL-like dissipation parameters form Militello 2013
+#mu_n0                 = 0.0121  # Will use nominal value of ~0.01
+#mu_vort0              = 5.0
+mu_vort_core           = 5.0
+#mu                    = 1000
+#nu_parallel0          = 0.07
+#kappa0                = 0.0440
+#kappa0_perp           = 0.0440
+#g0                    = 0.002862
+
+ixsep                 = 200                  # Radial position of separatrix
+
+n_bg                  = 0.01
+T_bg                  = 0.01
+
+L                     = mesh:Ly
+bracket               = 2                   # 0 = std, 1 = simple, 2 = arakawa
+isothermal            = false                # switch for isothermal simulations
+boussinesq            = true                # switch for solution with boussinesq approx (without requires multigrid)
+blob_sim              = false               # switch for blob_sim
+SOL_closure           = heuristic           # choice of 2d closure
+sheath_linear         = false               # switch to linearise sheath terms
+unit_bg               = false                # loss terms cause decay of fluctuations back to a background of n=1
+initial_noise         = true                # switch to add random noise to trigger instabilities
+uniform_diss_paras    = true
+driftwaves            = false
+
+[All]
+#scale = 0.0 
+#xs_opt = 3
+#zs_opt = 3
+#bndry_all = neumann
+
+[n]
+scale = 1.0
+sep_loc = 0.33333
+#sep_loc = storm:ixsep/(mesh:nx-4)
+falloff = 0.20
+bg = storm:n_bg
+#         initial bg in core  initial bg in sep     initial exp decay in sep
+function = H(sep_loc-x)*1.0 + H(x-sep_loc)*(bg + (1-bg)*exp(-(x-sep_loc)/falloff))
+bndry_xin = neumann_o2(0.0)
+bndry_xout = neumann_o2(0.0)
+
+[T]
+scale = 1.0
+sep_loc = 0.33333
+#sep_loc = storm:ixsep/(mesh:nx-4)
+falloff = 0.20
+bg = storm:T_bg
+#   	     initial bg in core  initial bg in sep     initial exp decay in sep
+function = H(sep_loc-x)*1.0 + H(x-sep_loc)*(bg + (1-bg)*exp(-(x-sep_loc)/falloff))
+bndry_xin = neumann_o2(0.0)     # Fix logT = 0.0, so T = 1.0 on inner boundary
+bndry_xout = neumann_o2(0.0)
+
+[vort]
+scale = 1.0
+function = 0.0
+#xs_opt = 3
+#zs_opt = 3
+bndry_xin = neumann(0.0)
+bndry_xout = neumann(0.0) 
+
+[B]         # Magnetic field 
+scale = 1.0
+function = 1.0 ; 
+
+##################################################################
+[S]       # Density source
+Ly                    = mesh:Ly
+scale                 = 1.0             # No sources, instead n & T are fixed to non-zero values on inner boundary
+turb_smag             = 11.0
+turb_swidth           = 0.01666666
+s_loc                 = 0.1                                                                                    # Fraction of distance along x to place centre of turbulence source
+function              = turb_smag*exp(-((x-s_loc)/turb_swidth)^2.0) / Ly
+
+[S_E]       # Temperature source
+Ly                    = mesh:Ly
+scale                 = 1.0             # No sources, instead n	& T are	fixed to non-zero values on inner boundary
+turb_smag             = 20.0    
+turb_swidth           = 0.01666666
+s_loc                 = 0.1                                                                                    # Fraction of distance along x to place centre of turbulence source
+function              = turb_smag*exp(-((x-s_loc)/turb_swidth)^2.0) / Ly
+
+[uE]
+bndry_xin             = neumann_o2
+bndry_xout            = neumann_o2
+
diff --git a/xbout/tests/test_options.py b/xbout/tests/test_options.py
new file mode 100644
index 00000000..fdc2e7b0
--- /dev/null
+++ b/xbout/tests/test_options.py
@@ -0,0 +1,353 @@
+from pathlib import Path
+from textwrap import dedent
+
+import numpy as np
+
+import pytest
+
+from xbout.options import Section, OptionsTree, OptionsFile
+
+
+@pytest.fixture
+def example_section():
+    contents = {'type': 'cyclic',
+                'global_flags': '0',
+                'inner_boundary_flags': '1',
+                'outer_boundary_flags': '16  # dirichlet',
+                'include_yguards': 'false  ; another comment'}
+    return Section(name='laplace', data=contents)
+
+
+class TestSection:
+    def test_create_section(self, example_section):
+        sect = example_section
+        assert sect.name == 'laplace'
+        assert sect.data == {'type': 'cyclic',
+                             'global_flags': '0',
+                             'inner_boundary_flags': '1',
+                             'outer_boundary_flags': '16  # dirichlet',
+                             'include_yguards': 'false  ; another comment'}
+
+    def test_get(self, example_section):
+        sect = example_section
+        assert sect.get('type') == 'cyclic'
+
+    def test_get_comments(self, example_section):
+        sect = example_section
+        # Comments marked by '#'
+        assert sect.get('outer_boundary_flags') == '16'
+        with_comment = sect.get('outer_boundary_flags', evaluate=False,
+                                keep_comments=True)
+        assert with_comment == '16  # dirichlet'
+
+        # Comments marked by ';'
+        assert sect.get('include_yguards') == 'false'
+        with_other_comment = sect.get('include_yguards', evaluate=False,
+                                      keep_comments=True)
+        assert with_other_comment == 'false  ; another comment'
+
+    def test_getitem(self, example_section):
+        sect = example_section
+        assert sect['type'] == 'cyclic'
+        assert sect['outer_boundary_flags'] == '16'
+
+    def test_set(self, example_section):
+        sect = example_section
+        sect.set('max', '10000')
+        assert sect.get('max') == '10000'
+
+        sect['min'] = '100'
+        assert sect.get('min') == '100'
+
+    def test_set_nested(self, example_section):
+        sect = example_section
+        nested = Section('naulin', data={'iterations': '1000'})
+        sect['naulin'] = nested
+
+        assert isinstance(sect.get('naulin'), Section)
+        assert sect.get('naulin').get('iterations') == '1000'
+        assert sect['naulin']['iterations'] == '1000'
+
+    def test_set_nested_as_dict(self, example_section):
+        sect = example_section
+        sect['naulin'] = {'iterations': '1000'}
+
+        assert isinstance(sect.get('naulin'), Section)
+        assert sect.get('naulin').get('iterations', evaluate=False) == '1000'
+        assert sect['naulin']['iterations'] == '1000'
+        assert sect['naulin'].parent == 'laplace'
+
+    def test_find_parents(self, example_section):
+        sect = example_section
+        sect['naulin'] = {'iterations': '1000'}
+
+        assert sect['naulin'].lineage() == 'laplace:naulin'
+
+    def test_print(self, example_section):
+        sect = example_section
+        sect['naulin'] = {'iterations': '1000'}
+        expected = dedent("""\
+                          [laplace]
+                          |-- type = cyclic
+                          |-- global_flags = 0
+                          |-- inner_boundary_flags = 1
+                          |-- outer_boundary_flags = 16
+                          |-- include_yguards = false
+                          |-- [naulin]
+                          |-- |-- iterations = 1000
+                          """)
+        assert str(sect) == expected
+
+    def test_write(self, example_section, tmpdir_factory):
+        sect = example_section
+        sect['naulin'] = {'iterations': '1000'}
+        file = tmpdir_factory.mktemp("write_data").join('write_test.inp')
+
+        with open(file, 'w') as f:
+            sect._write(f)
+
+        with open(file, 'r') as f:
+            result = file.read()
+        expected = dedent("""\
+                          [laplace]
+                          type = cyclic
+                          global_flags = 0
+                          inner_boundary_flags = 1
+                          outer_boundary_flags = 16  # dirichlet
+                          include_yguards = false  ; another comment
+                          
+                          [laplace:naulin]
+                          iterations = 1000
+                          """)
+        assert result == expected
+
+    def test_write_root(self, tmpdir_factory):
+        sect = Section(name='root', parent=None,
+                       data={'timestep': '80  # Timestep length',
+                             'nout': '2000  # Number of timesteps'})
+        file = tmpdir_factory.mktemp("write_data").join('write_test.inp')
+
+        with open(file, 'w') as f:
+            sect._write(f)
+
+        with open(file, 'r') as f:
+            result = file.read()
+        expected = dedent("""\
+                          timestep = 80  # Timestep length
+                          nout = 2000  # Number of timesteps
+                          """)
+        assert result == expected
+
+
+@pytest.fixture
+def example_options_tree():
+    """
+    Mocks up a temporary example BOUT.inp (or BOUT.settings) tree
+    """
+
+    # Create options file
+    options = OptionsTree()
+
+    # Fill it with example data
+    options.data = {'timestep': '80',
+                    'nout': '2000    # Number of outputted timesteps'}
+    options['mesh'] = {'Ly': '2000.0   # ~4m',
+                       'Lx': '400.0',
+                       'nx': '484      # including 4 guard cells',
+                       'ny': '1        # excluding guard cells',
+                       'dx': 'Lx/(nx-4)',
+                       'dy': 'Ly/ny'}
+    options['mesh']['ddx'] = {'first': 'C2',
+                              'second': 'C2',
+                              'upwind': 'C2'}
+    options['laplace'] = {'type': 'cyclic',
+                          'global_flags': '0',
+                          'inner_boundary_flags': '1',
+                          'outer_boundary_flags': '16',
+                          'include_yguards': 'false'}
+    options['storm'] = {'B_0': '0.24          # Tesla',
+                        'T_e0': '15           # eV',
+                        'isothermal': 'true   # switch'}
+
+    return options
+
+"""
+def example_options_file():
+    # Create temporary directory
+    save_dir = tmpdir_factory.mktemp("inputdata")
+
+    # Save
+    optionsfilepath = save_dir.join('BOUT.inp')
+    options.write_to(optionsfilepath)
+
+    # Remove the first (and default) section headers to emulate BOUT.inp file format
+    # TODO do this without opening file 3 times
+    with open(optionsfilepath, 'r') as fin:
+        data = fin.read().splitlines(True)
+    with open(optionsfilepath, 'w') as fout:
+        fout.writelines(data[4:])
+
+    return Path(optionsfilepath)
+"""
+
+
+class TestAccess:
+    def test_get_sections(self, example_options_tree):
+        sections = OptionsTree(example_options_tree).sections()
+        lineages = [section.lineage() for section in sections]
+        assert lineages == ['root', 'root:mesh', 'root:mesh:ddx',
+                            'root:laplace', 'root:storm']
+
+    def test_get_str_values(self, example_options_tree):
+        opts = OptionsTree(example_options_tree)
+        assert opts['laplace']['type'] == 'cyclic'
+        assert opts['laplace'].get('type') == 'cyclic'
+
+    def test_get_nested_section_values(self, example_options_tree):
+        opts = OptionsTree(example_options_tree)
+        assert opts['mesh']['ddx']['upwind'] == 'C2'
+        assert opts['mesh']['ddx'].get('upwind') == 'C2'
+
+
+@pytest.mark.xfail
+class TestReadFile:
+    def test_no_file(self):
+        with pytest.raises(FileNotFoundError):
+            OptionsFile('./wrong.inp')
+
+    def test_open_real_example(self):
+        # TODO this absolute filepath is sensitive to the directory the tests are run from?
+        example_inp_path = Path.cwd() / './xbout/tests/data/options/BOUT.inp'
+        opts = OptionsFile(example_inp_path)
+        assert opts.filepath.name == 'BOUT.inp'
+
+    def test_open(self, example_options_file):
+        opts = OptionsFile(example_options_file)
+        assert opts.filepath.name == 'BOUT.inp'
+
+    def test_repr(self, example_options_file):
+        opts_repr = repr(OptionsFile(example_options_file))
+        assert opts_repr == f"OptionsFile('{example_options_file}')"
+
+
+class TestTypeConversion:
+    def test_get_int_values(self, example_options_tree):
+        opts = OptionsTree(example_options_tree)
+        assert opts['mesh'].get('nx', evaluate=True) == 484
+
+    def test_get_float_values(self, example_options_tree):
+        opts = OptionsTree(example_options_tree)
+        assert opts['mesh'].get('Lx', evaluate=True) == 400.0
+
+    def test_get_bool_values(self, example_options_tree):
+        opts = OptionsTree(example_options_tree)
+        assert opts['storm'].get('isothermal', evaluate=True) == True
+
+
+class TestEvaluation:
+    def test_eval_arithmetic(self):
+        params = Section(name='params', data={'T0': '5 * 10**18'})
+        assert params.get('T0', evaluate=True) == 5 * 10**18
+
+    def test_eval_powers(self):
+        params = Section(name='params', data={'n0': '10^18'})
+        assert params.get('n0', evaluate=True) == 10**18
+
+    # TODO generalise this to test other numpy imports
+    def test_eval_numpy(self):
+        params = Section(name='params', data={'angle': 'sin(pi/2)'})
+        assert params.get('angle', evaluate=True) == np.sin(np.pi/2)
+
+
+@pytest.fixture
+def example_options_tree_with_substitution():
+    """
+    Mocks up a temporary example BOUT.inp (or BOUT.settings) tree,
+    containing variables which need to be substituted with values from
+    other sections
+    """
+
+    # Create options file
+    options = OptionsTree()
+
+    # Fill it with example data
+    options.data = {'timestep': '80',
+                    'nout': '2000    # Number of outputted timesteps'}
+
+    # Substitutions only within section
+    options['mesh'] = {'Ly': '2000.0   # ~4m',
+                       'Lx': '400.0',
+                       'nx': '484      # including 4 guard cells',
+                       'dx': 'Lx/(nx-4)',
+                       'dy': 'Ly'}
+
+    # Substitute from section above
+    options['mesh']['ddx'] = {'first': 'C2',
+                              'length': 'mesh:Lx*2',
+                              'num': '3'}
+
+    # TODO Substitute from section below
+    options['model'] = {}
+    #options['model']['energy'] = ''
+    #options['model']['magnetic_energy'] = 'B_0^2'
+
+    # Substitute from arbitrary section
+    options['model']['storm'] = {'B_0': '0.24          # Tesla',
+                                 'T_e0': '15           # eV',
+                                 'L_parallel': 'mesh:ddx:num * 10',
+                                 'isothermal': 'true   # switch'}
+
+    # Self-referencing substution
+    options['apples'] = {'green': 'lemons:yellow'}
+    options['lemons'] = {'yellow': 'apples:green'}
+
+    return options
+
+
+@pytest.mark.xfail  
+class TestSubstitution:
+    def test_substitute_from_section(self, example_options_tree_with_substitution):
+        opts = example_options_tree_with_substitution
+        ly = opts['mesh']['Ly']
+        assert opts['mesh'].get('dy', substitute=True) == ly
+
+    def test_substitute_from_parent_section(self, example_options_tree_with_substitution):
+        opts = example_options_tree_with_substitution
+        lx = opts['mesh']['Lx']
+        assert opts['mesh']['ddx'].get('length', substitute=True) == f"{lx}*2"
+
+    @pytest.mark.skip
+    def test_substitute_from_child_section(self, example_options_tree_with_substitution):
+        opts = example_options_tree_with_substitution
+        lx = opts['mesh']['Lx']
+        assert opts['mesh']['ddx'].get('length', substitute=True) == f"{lx}*2"
+
+    def test_substitute_from_arbitrary_section(self, example_options_tree_with_substitution):
+        opts = example_options_tree_with_substitution
+        ddx = opts['mesh']['ddx']['num']
+        assert opts['model']['storm'].get('L_parallel', substitute=True) == f"{ddx} * 10"
+
+    def test_two_substitutions_required(self, example_options_tree_with_substitution):
+        opts = example_options_tree_with_substitution
+        assert opts['mesh'].get('dx', substitute=True) == '400.0/(484-4)'
+
+    def test_substitute_requires_substitution(self, example_options_tree_with_substitution):
+        ...
+
+    @pytest.mark.skip
+    def test_contains_colon_but_no_substitute_found(self, example_options_tree_with_substitution):
+        with pytest.raises(InterpolationError): #?
+            ...
+
+    def test_detect_infinite_substitution_cycle(self, example_options_tree_with_substitution):
+        ...
+
+
+@pytest.mark.skip
+class TestWriting:
+    ...
+
+
+@pytest.mark.skip
+class TestRoundTrip:
+    ...