Repository: myokit/myokit (language: python)

Issue:
ABF reader scaling wrong in some files
The scaling factor in the voltage signal (second channel, AD1) of `c260315_0000.abf` (model-cell-3) is off by a factor ten. It works properly in clampfit.



Discussion / hints:
Strange. The test case is a v1.65 file, so should be able to see exact code in support pack.

```
//==============================================================================================
// FUNCTION:   GetADCtoUUFactors
// PURPOSE:    Calculates the scaling factors used to convert ADC values to UserUnits.
// PARAMETERS:
//    nChannel        - The physical channel number to get the factors for.
//    pfADCToUUFactor - Pointers to return locations for scale and offset.
//    pfADCToUUShift    UserUnits = ADCValue * fADCToUUFactor + fADCToUUShift;
//
void WINAPI ABFH_GetADCtoUUFactors( const ABFFileHeader *pFH, int nChannel, 
                                    float *pfADCToUUFactor, float *pfADCToUUShift )
{
   ABFH_ASSERT(pFH);
   WPTRASSERT(pfADCToUUFactor);
   WPTRASSERT(pfADCToUUShift);
   ASSERT(nChannel < ABF_ADCCOUNT);

   float fTotalScaleFactor = pFH->fInstrumentScaleFactor[nChannel] *
                             pFH->fADCProgrammableGain[nChannel];
   if (pFH->nSignalType != 0)
      fTotalScaleFactor *= pFH->fSignalGain[nChannel];

   // Adjust for the telegraphed gain.
   if( pFH->nTelegraphEnable[nChannel] )
      fTotalScaleFactor *= pFH->fTelegraphAdditGain[nChannel];

   ASSERT(fTotalScaleFactor != 0.0F);
   if (fTotalScaleFactor==0.0F)
      fTotalScaleFactor = 1.0F;

   // InputRange and InputOffset is the range and offset of the signal in
   // user units when it hits the Analog-to-Digital converter

   float fInputRange = pFH->fADCRange / fTotalScaleFactor;
   float fInputOffset= -pFH->fInstrumentOffset[nChannel];
   if (pFH->nSignalType != 0)
      fInputOffset += pFH->fSignalOffset[nChannel];

   *pfADCToUUFactor = fInputRange / pFH->lADCResolution;
   *pfADCToUUShift  = -fInputOffset;
}
```

apart from data types etc that's the same as what the current code does

Need to use `h['nADCSamplingSeq'][i]` instead of `i` when reading meta data for channel `i` in v1

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index f5da2fd3d..04b3a7a49 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,23 +4,26 @@ This page lists the main changes made to Myokit in each release.
 
 ## Unreleased
 - Added
-  - [#985](https://github.com/myokit/myokit/pull/985) The `AbfFile` and `WcpFile` now implement a shared `SweepSource` interface.
+  - [#985](https://github.com/myokit/myokit/pull/985), [#998](https://github.com/myokit/myokit/pull/998) The `AbfFile` and `WcpFile` now implement a shared `SweepSource` interface.
   - [#985](https://github.com/myokit/myokit/pull/985) Added a method `AtfFile.path` that returns the full path.
   - [#985](https://github.com/myokit/myokit/pull/985) Added a method `WcpFile.info` that returns a string version of the file's meta data, and a method `WcpFile.sample_count` that returns the number of samples per channel.
 - Changed
-  - [#985](https://github.com/myokit/myokit/pull/985) The interface to the `AbfFile` class has changed significantly, and in a backwards incompatible ways. Please see the documentation for further details.
+  - [#985](https://github.com/myokit/myokit/pull/985), [#998](https://github.com/myokit/myokit/pull/998)  The interface to the `AbfFile` class has changed significantly, and in a backwards incompatible ways. Please see the documentation for further details.
   - [#991](https://github.com/myokit/myokit/pull/991) Myokit no longer uses the `logging` module but makes calls to `warnings.warn` instead.
   - [#992](https://github.com/myokit/myokit/pull/992) If a `Simulation` is created with `protocol={}`, then `pace` will not be set as a binding label. To stay backwards compatible it _will_ be set if `protocol=None`.
   - [#994](https://github.com/myokit/myokit/pull/994) Simulations can now start at arbitrary times, including negative ones. Note that protocol events have a default offset of 0, so no events are usually expected at negative times.
 - Deprecated
-  - [#985](https://github.com/myokit/myokit/pull/985) The method `AtfFile.myokit_log` is deprecated in favor of `AtfFile.log`.
   - [#985](https://github.com/myokit/myokit/pull/985) The method `WcpFile.myokit_log` is deprecated in favor of `WcpFile.log`.
   - [#985](https://github.com/myokit/myokit/pull/985) The methods `WcpFile.channels` and `WcpFile.records` are deprecated in favor of `WcpFile.channel_count` and `WcpFile.record_count`.
   - [#985](https://github.com/myokit/myokit/pull/985) The method `WcpFile.plot` is deprecated in favor of `WcpFile.matplotlib_figure`.
+  - [#985](https://github.com/myokit/myokit/pull/985) The method `AtfFile.myokit_log` is deprecated in favor of `AtfFile.log`.
+  - [#998](https://github.com/myokit/myokit/pull/998) The method `WcpFile.info` is deprecated in favor of `WcpFile.meta_str`.
+  - [#998](https://github.com/myokit/myokit/pull/998) The method `AtfFile.info` is deprecated in favor of `AtfFile.meta_str`.
 - Fixed
   - [#985](https://github.com/myokit/myokit/pull/985) `AbfFile.filename` now returns the filename, not the full path.
   - [#985](https://github.com/myokit/myokit/pull/985) `AtfFile.filename` now returns the filename, not the full path.
   - [#992](https://github.com/myokit/myokit/pull/992) `Simulation` now sets any variable bound to a pacing label to 0 if no protocol for that label is set.
+  - [#998](https://github.com/myokit/myokit/pull/998) Fixed a bug where channels in ABF v1 files were sometimes scaled using the wrong meta data.
 
 ## [1.35.0] - 2023-06-21
 - Added
diff --git a/docs/source/api_gui/datalog_viewer.rst b/docs/source/api_gui/datalog_viewer.rst
index 3902764d1..e6137b5c7 100644
--- a/docs/source/api_gui/datalog_viewer.rst
+++ b/docs/source/api_gui/datalog_viewer.rst
@@ -21,6 +21,13 @@ Widgets
 
 .. autoclass:: TabWidget
 
+.. autoclass:: GraphTabWidget
+
+.. autoclass:: SweepSourceTab
+
+Graphing widgets
+----------------
+
 .. autoclass:: AbfTab
 
 .. autoclass:: AtfTab
@@ -32,3 +39,4 @@ Widgets
 .. autoclass:: TxtTab
 
 .. autoclass:: WcpTab
+
diff --git a/docs/source/api_index/myokit.gui.rst b/docs/source/api_index/myokit.gui.rst
index 9d82d0323..402f65572 100644
--- a/docs/source/api_index/myokit.gui.rst
+++ b/docs/source/api_index/myokit.gui.rst
@@ -31,7 +31,9 @@ myokit.gui.datalog_viewer
 - :class:`myokit.gui.datalog_viewer.AtfTab`
 - :class:`myokit.gui.datalog_viewer.CsvTab`
 - :class:`myokit.gui.datalog_viewer.DataLogViewer`
+- :class:`myokit.gui.datalog_viewer.GraphTabWidget`
 - :class:`myokit.gui.datalog_viewer.MatTab`
+- :class:`myokit.gui.datalog_viewer.SweepSourceTab`
 - :class:`myokit.gui.datalog_viewer.TabWidget`
 - :class:`myokit.gui.datalog_viewer.TxtTab`
 - :class:`myokit.gui.datalog_viewer.WcpTab`
diff --git a/myokit/_unit.py b/myokit/_unit.py
index ded9309ef..0fcb0cdd5 100644
--- a/myokit/_unit.py
+++ b/myokit/_unit.py
@@ -258,7 +258,7 @@ def conversion_factor(unit1, unit2, helpers=None):
 
         Returns a :class:`myokit.Quantity`.
 
-        Raises a :class:`myokit.IncompatibleUnitError` if the units cannot be
+        Raises an :class:`myokit.IncompatibleUnitError` if the units cannot be
         converted.'
         """
         # Check unit1
diff --git a/myokit/formats/__init__.py b/myokit/formats/__init__.py
index b0187414c..9589dbb05 100644
--- a/myokit/formats/__init__.py
+++ b/myokit/formats/__init__.py
@@ -668,39 +668,50 @@ def register_external_ewriter(name, ewriter_class):
 class SweepSource:
     """
     Interface for classes that provide time-series data organised into *sweeps*
-    (or *records*) and *channels* (or *traces*).
+    (or *records* or *episodes*) and *channels* (or *traces*).
 
-    Each sweep contains the same number of channels, and each channel is a 1d
-    array of equal length.
+    The :class:`SweepSource` interface defines methods to get the number of
+    sweeps and channels, the names and units of the channels, and the data
+    stored in channels either as numpy arrays or in a :class:`myokit.DataLog`.
+
+    Each sweep contains the same number of channels, and each channel is
+    represented as a 1d array. In most cases these arrays have the same length
+    for every sweep, but whether this is the case for the current source can be
+    tested with :meth:`equal_length_sweeps`.
 
-    Two common ways to represent the data are (1) overlaid::
+    Data can be retrieved sweep by sweep::
 
-        time        n_t data points
-        sweep[0],   n_c channels, each with n_t data points
-        sweep[1],
+        (time[0], sweep[0]),        n_t data points, n_c channels of n_t points
+        (time[1], sweep[1]),
+        (time[2], sweep[2]),
         ...
-        sweep[n_s]
+        (time[n_s], sweep[n_s])
 
-    or (2) joined into a single time series::
+    This allows plotting in the common "overlaid" fashion, i.e. plotting every
+    ``sweep[i] against ``time[0]``.
 
-        time                            n_s*n_t data points
-        sweep[0] + sweep[1] + ...       n_c channels, each with n_s*n_t points
+    For other types of analysis (e.g. parameter estimation) the data can also
+    be returned as single time series::
+
+        time                            sum(n_t[i]) data points
+        sweep[0] + sweep[1] + ...       n_c channels, with sum(n_t[i]) points
+
+    Some formats can also contain information from which D/A output signals can
+    be reconstructed. These can be accessed using the :meth:`da`.
 
-    The :class:`SweepSource` interface defines methods to get the number of
-    sweeps and channels, the names of the channels, and the data stored in
-    channels either as numpy arrays or in a :class:`myokit.DataLog`.
     """
     def channel(self, channel_id, join_sweeps=False):
         """
         Returns the data for a single channel, identified by integer or string
         ``channel_id``.
 
-        With ``join_sweeps=False``, the data is returned as a tuple of numpy
-        arrays ``time, sweep_0, sweep_1, ...`` where ``time`` is the time
-        corresponding to the first sweep.
+        With ``join_sweeps=False``, the data is returned as a tuple
+        ``(times, sweeps)`` where ``times`` and ``sweeps`` are 2d numpy arrays
+        indexed so that ``times[i][j]`` is the ``j``-th time point for sweep
+        ``i``.
 
         If ``join_sweeps=True`` the sweeps are joined together, and a tuple
-        ``time, values`` is returned.
+        ``(times, values)`` is returned ``times`` and ``values`` are 1d arrays.
         """
         raise NotImplementedError
 
@@ -708,32 +719,123 @@ def channel_count(self):
         """ Returns the number of channels. """
         raise NotImplementedError
 
-    def channel_names(self):
-        """ Returns a list of names for each channel. """
+    def channel_names(self, index=None):
+        """
+        Returns the names of all channels or the name of a specific channel
+        ``index``.
+        """
+        raise NotImplementedError
+
+    def channel_units(self, index=None):
+        """
+        Returns the units (as :class:`myokit.Unit`) of all channels or the
+        units of a specific channel ``index``.
+        """
+        raise NotImplementedError
+
+    def da(self, output_id, join_sweeps=False):
+        """ Like :meth:`channel`, but returns reconstructed D/A signals. """
+        raise NotImplementedError
+
+    def da_count(self):
+        """
+        Returns the available number of reconstructed D/A output channels.
+
+        This should return 0 if D/A channels are not supported.
+        """
+        raise NotImplementedError
+
+    def da_names(self):
+        """
+        Returns the names of all reconstructed D/A output channels or the name
+        of a specific output channel ``index``.
+
+        This will raise a ``NotImplementedError`` if D/A channels are not
+        supported.
+        """
+        raise NotImplementedError
+
+    def da_protocol(self, output_id=None, tu='ms', vu='mV', cu='pA',
+                    n_digits=9):
+        """
+        Returns a :class:`myokit.Protocol` representation of the D/A output
+        channel specified by name or integer ``output_id``.
+
+        If no explicit output channel is set, a guess will be made.
+
+        Time units will be converted to ``tu``, voltage units to ``vu``, and
+        current units to ``cu``. Other unit types will be left unchanged. Units
+        may be specified as :class:`myokit.Unit` objects or strings.
+
+        By default, floating point numbers in the protocol will be rounded to
+        9 digits after the decimal point. This can be changed or disabled via
+        the ``n_digits`` argument (which is passed to the built-in method
+        ``round``).
+
+        If a D/A outputs cannot be converted to a :class:`myokit.Protocol`, a
+        ``ValueError`` is raised. A ``NotImplementedError`` is raised if D/A
+        channels are not supported at all.
+        """
         raise NotImplementedError
 
-    def log(self, join_sweeps=False, use_names=False, channels=None):
+    def da_units(self):
+        """
+        Returns the units (as :class:`myokit.Unit`) of all reconstructed D/A
+        output channels or the units of a specific output channel ``index``.
+
+        This will raise a ``NotImplementedError`` if D/A channels are not
+        supported.
+        """
+        raise NotImplementedError
+
+    def equal_length_sweeps(self):
+        """
+        Returns ``True`` only if each sweep in this source has the same length.
+        """
+        raise NotImplementedError
+
+    def log(self, join_sweeps=False, use_names=False, include_da=True):
         """
         Returns a :class:`myokit.DataLog` containing the data from all
-        channels or the subset specified by ``channels``.
+        channels.
+
+        The log will have a single entry ``time`` corresponding to the time of
+        the first sweep if ``join_sweeps`` is ``False``, or the time of all
+        points when ``join_sweeps`` is ``True``.
 
-        Sweeps can be joined together using ``join_sweeps=True``, and log
-        entries can use systematic names (``use_names=False``) or names
-        specified in the original source (``use_names=True``).
+        Names will have a systematic form ``sweep.channel.type``, for example
+        ``0.1.channel`` for the 1st sweep, second channel of A/D input, or
+        ``3.0.da`` for the 3d sweep, first channel of reconstructed D/A output.
+        To obtain a log with the user-specified names from the source instead,
+        set ``use_names`` to ``True``.
 
-        If ``join_sweeps=False``, the data for each sweep will be returned
-        separately, so that the log will have entries such as ``0.0.channel``,
-        ``1.0.channel``, and ``2.0.channel`` for the first three sweeps and the
-        first channel. Alternatively, if ``use_names=True`` and the first
-        channel is called "IN #0" its entries will be ``0.IN #0`` for the first
-        sweep, then ``1.IN #0`` etc.
+        To exclude D/A signal reconstructions, set ``include_da`` to ``False``.
 
-        To return only a subset of channels pass in a list of channel ids
-        (names or integers) as ``channels``.
+        A call with ``join_sweeps=False`` on a source where
+        :meth:`equal_length_sweeps()` returns ``False`` will raise a
+        ``ValueError``.
         """
         raise NotImplementedError
 
+    def meta_str(self, verbose=False):
+        """
+        Optional method that returns a multi-line string with unstructured meta
+        data about the source and its contents.
+
+        Will return ``None`` if no such string is available.
+        """
+        return None  # pragma: no cover
+
     def sweep_count(self):
-        """ Returns the number of sweeps. """
+        """
+        Returns the number of sweeps.
+
+        Note that a source with zero recorded channels may still report a
+        non-zero number of sweeps if it can provide D/A outputs.
+        """
+        raise NotImplementedError
+
+    def time_unit(self):
+        """ Returns the time unit used in this source. """
         raise NotImplementedError
 
diff --git a/myokit/formats/axon/_abf.py b/myokit/formats/axon/_abf.py
index 3b816d8d3..91960bfc1 100755
--- a/myokit/formats/axon/_abf.py
+++ b/myokit/formats/axon/_abf.py
@@ -1,8 +1,8 @@
 #
 # This module reads files in Axon Binary File format v1 or v2 used by Axon
 # Technologies and Molecular Devices.
-# The v1 format was used until Clampex version 9. Clampex 10 and onwards use
-# the v2 format.
+# The v1 format was used until Clampex version 9.
+# Clampex 10 and onwards use the v2 format.
 #
 # WARNING: This file hasn't been extensively tested.
 #
@@ -57,13 +57,6 @@
 #  an 'epoch'. The protocol section of a file defines one or more stimuli, each
 #  containing a list of epochs.
 #
-# Conversion to myokit formats
-# ----------------------------
-# There is no problem-free mapping of ABF data onto myokit structures, such as
-# the simulation log. A fundamental difference is that "sampling" during a
-# simulation happens at the same time for every signal. Channels in an ABF file
-# each have their own sampling rate.
-#
 #---------------------------------  license  ----------------------------------
 #
 # This file is part of Myokit.
@@ -71,12 +64,18 @@
 #
 #---------------------------------  credits  ----------------------------------
 #
-# This module is derived in part from code found in the Neo package for
-# representing electrophysiology data, specifically from a python module
-# authored by sgarcia and jnowacki.
+# A lot of this code uses the (sadly somewhat outdated) information made public
+# by Axon, e.g. at https://support.moleculardevices.com/s/article/
+#                               Axon-pCLAMP-ABF-File-Support-Pack-Download-Page
+# This information comes without a specific license, but states that
+# "Permission is granted to freely use, modify and copy the code in this file."
+#
+# In addition, this module was in part derived from an early version of the
+# Neo package for representing electrophysiology data, specifically from a
+# Python module authored by sgarcia and jnowacki.
 # Neo can be found at: http://neuralensemble.org/trac/neo
 #
-# The Neo package is licensed using the following BSD License:
+# The Neo package used was licensed using the following BSD License:
 #
 #----------------------------------  start  -----------------------------------
 # Copyright (c) 2010-2012, Neo authors and contributors
@@ -146,8 +145,6 @@
 
 import numpy as np
 
-from collections import OrderedDict
-
 import myokit
 import myokit.formats
 
@@ -168,7 +165,7 @@ class AbfFile(myokit.formats.SweepSource):
 
     Data in ABF files is recorded in *sweeps*, where each sweep contains one or
     more *channels* with recorded (A/D) data. In addition, zero or more output
-    channels may be defined (also called "protocol" or D/A channels). Where
+    waveforms may be defined (also called "protocol" or D/A channels). Where
     possible, the :class`AbfFile` class will convert these embedded protocols
     to time series and include them as additional channels.
 
@@ -186,9 +183,8 @@ class AbfFile(myokit.formats.SweepSource):
         10xVm
         Cmd 0
 
-    where the first two channels are A/D channels and the final one is a
-    converted D/A channel. *Some* A/D channels (see below) can be converted to
-    :class:`myokit.Protocol` objects, using the :meth:`protocol()` method.
+    where the first two channels are recorded A/D channels and the final one is
+    a reconstructed D/A output channel.
 
     Sweeps and channels are represented by :class:`Sweep` and :class:`Channel`
     objects respectively, and these can be used to obtain the data from a
@@ -200,15 +196,18 @@ class AbfFile(myokit.formats.SweepSource):
                 plot(channel.times(), channel.values())
 
     In addition the ``AbfFile`` class implements the
-    :class`myokit.formats.SweepSource` interface.
+    :class`myokit.formats.SweepSource` interface. Note that this interface
+    treats A/D and D/A as separate things, so :meth:`channel_count` returns the
+    number of A/D channels, not the total number of channels in a
+    :class:`Sweep` object (which can include D/A channels).
 
     Support notes:
 
     - Protocol (D/A) conversion is only supported for "episodic stimulation"
-      without "user lists".
+      with constant valued steps (so e.g. no ramps) and without "user lists".
     - Protocols with more than one sampling rate are not supported.
-    - The ABF format is not (openly) well documented, so there will be several
-      other issues and shortcomings.
+    - The publicly available information on the ABF format is not great, so
+      there will be several other issues and shortcomings.
 
     Arguments:
 
@@ -218,8 +217,8 @@ class AbfFile(myokit.formats.SweepSource):
     ``is_protocol_file``
         If set to ``True``, no attempt to read A/D data will be made and only
         D/A "protocol" information will be read. If left at its default value
-        of ``None`` files with the extension ``.pro`` will still be recognized
-        as protocol files.
+        of ``None`` files with the extension ``.pro`` will be recognized as
+        protocol files.
 
     """
     def __init__(self, filepath, is_protocol_file=None):
@@ -228,474 +227,87 @@ def __init__(self, filepath, is_protocol_file=None):
         self._filepath = os.path.abspath(filepath)
         self._filename = os.path.basename(filepath)
 
-        # Abf format version
-        self._version = None
-        self._version_str = None
-
-        # Protocol info
-        self._epoch_functions = None
-        self._number_of_trials = None
-        self._trial_start_to_start = None
-        self._runs_per_trial = None
-        self._run_start_to_start = None
-        self._sweeps_per_run = None
-        self._sweep_start_to_start = None
-
         # Read as protocol file yes?
         if is_protocol_file is None:
             self._is_protocol_file = os.path.splitext(filepath)[1] == '.pro'
         else:
             self._is_protocol_file = bool(is_protocol_file)
 
-        # The file header (an ordered dictionary)
-        self._header = self._read_header()
-
-        # Date/time of recording
-        self._datetime = self._read_datetime()
-
-        # Number of channels, sampling rate (Hz) and acquisition mode
-        # Note: Number of AD channels will be set to 0 if this is a
-        #       protocol-only file
-        # Note: Number of DA channels will be adjusted by _read_1_protocol
-        if self._version < 2:
-            self._nADC = int(self._header['nADCNumChannels'])
-            self._nDAC = len(self._header['sDACChannelName'])
-            self._rate = 1e6 / (
-                self._header['fADCSampleInterval'] * self._nADC)
-            self._mode = self._header['nOperationMode']
-        else:
-            self._nADC = int(self._header['sections']['ADC']['length'])
-            self._nDAC = int(self._header['sections']['DAC']['length'])
-            self._rate = 1e6 / self._header['protocol']['fADCSequenceInterval']
-            self._mode = self._header['protocol']['nOperationMode']
-        if self._mode not in acquisition_modes:
-            raise NotImplementedError(f'Unknown acquisition mode: {mode}')
-        if self._is_protocol_file:
-            self._nADC = 0
+        # Cached string-to-unit conversions
+        self._unit_cache = {}
 
-        # Conversion factors for integer data in the channels
-        self._adc_factors = None
-        self._adc_offsets = None
-        self._set_conversion_factors()
-
-        # Channel names: will be set on first call to channel_names()
-        self._channel_names = None
-
-        # Read the protocol info, and attempt to reconstruct the D/A signal as
-        # as list of sweeps. MUST be called before _read_ad_data.
-        self._sweeps = []
-        try:
-            self._read_1_protocol()
-        except Exception:  # pragma: no cover
-            # This is not something we _want_ to happen, so if we have test
-            # cases that trigger this error they should be resolved. At the
-            # same time, if it happens to a user we want it to "sort-of work"
-            # (an experimental rather than a production setting)
-            warnings.warn(f'Unable to read protocol from {self._filepath}')
-
-        # Create sweeps and add the AD data
-        if self._nADC:
-            self._read_2_data()
-
-        # Set final number of channels
-        self._nc = self._nADC + self._nDAC
-
-        # Store mapping from channel name to channel index
-        self._channel_names = {}
-        for sweep in self._sweeps:
-            for i, channel in enumerate(sweep):
-                self._channel_names[channel.name()] = i
-            break
-
-    def __getitem__(self, key):
-        return self._sweeps[key]
-
-    def __iter__(self):
-        return iter(self._sweeps)
-
-    def __len__(self):
-        return len(self._sweeps)
-
-    def ad_channel_count(self):
-        """ Returns the number of A/D channels in this file. """
-        return self._nADC
-
-    def channel(self, channel_id, join_sweeps=False):
-        # Docstring in SweepSource
-        channel_id = self._channel_id(channel_id)
-
-        if join_sweeps:
-            # Join sweeps
-            time, data = [], []
-            t = self._sweeps[0][channel_id].times()
-            for i, sweep in enumerate(self._sweeps):
-                time.append(t + i * self._sweep_start_to_start)
-                data.append(sweep[channel_id].values())
-            return (np.concatenate(time), np.concatenate(data))
-
-        else:
-            # Return individual sweeps
-            out = []
-            out.append(np.array(self._sweeps[0][channel_id].times()))
-            for sweep in self._sweeps:
-                out.append(np.array(sweep[channel_id].values()))
-            return tuple(out)
-
-    def channel_count(self):
-        # Docstring in SweepSource
-        return self._nADC + self._nDAC
-
-    def _channel_id(self, channel_id, must_be_da=False):
-        """ Checks an int or str channel id and returns a valid int. """
-        if len(self._sweeps) == 0:  # pragma: no cover
-            raise KeyError(f'Channel {channel_id} not found (empty file).')
-
-        # Handle string
-        if isinstance(channel_id, str):
-            int_id = self._channel_names[channel_id]  # Let KeyError go to user
-        else:
-            int_id = int(channel_id)  # TypeError for user
-            if int_id < 0 or int_id >= self._nc:
-                raise IndexError(f'channel_id out of range: {channel_id}')
-
-        # Check if this is a D/A protocol channel
-        if must_be_da and int_id < self._nADC:
-            raise ValueError(
-                'Given channel_id does not correspond to a D/A channel.')
-
-        return int_id
-
-    def channel_names(self):
-        # Docstring in SweepSource
-        return list(self._channel_names.keys())
-
-    def da_channel_count(self):
-        """ Returns the number of (non-empty) D/A channels in this file. """
-        return self._nDAC
-
-    def da_holding_level(self, channel_id):
-        """
-        Returns the holding level used by the requested output channel of the
-        embedded protocol.
-
-        Arguments:
-
-        ``channel_id``
-            The channel to return a protocol for, as an integer index or a
-            string name. If an integer is given it must be equal to or greater
-            than the value returned by :meth:`ad_channel_count()`.
-
-        """
-        ida = self._channel_id(channel_id, must_be_da=True) - self._nADC
-        dinfo, einfo_exists, einfo = self._epoch_functions
-        return dinfo(ida, 'fDACHoldingLevel')
-
-    def da_steps(self, channel_id):
-        """
-        Returns the steps seen during a D/A protocol, for episodic stimulation
-        (not including the holding potential).
-
-        For example, for a protocol that has holding level ``-120mV`` and
-        performs steps to ``-100mV``, then holding, then ``-80mV``, then
-        holding, and then ``-40mV`` the returned output will be::
-
-            ([-100, -80, -40])
-
-        For a more complicated protocol, where each step is followed by a step
-        down to ``-140mV``, the output would be::
-
-            ([-100, -80, -40], [-140, -140, -140])
-
-        Arguments:
-
-        ``channel_id``
-            The channel to return a protocol for, as an integer index or a
-            string name. If an integer is given it must be equal to or greater
-            than the value returned by :meth:`ad_channel_count()`.
-
-        """
-        ida = self._channel_id(channel_id, must_be_da=True) - self._nADC
-
-        # Get epoch functions
-        dinfo, einfo_exists, einfo = self._epoch_functions
-
-        # Create list of step lists
-        levels = []
-        for e in einfo(ida):
-            kind = e['type']
-            if kind not in epoch_types:
-                raise NotImplementedError(f'Unknown epoch type: {kind}')
-            if kind == EPOCH_DISABLED:
-                continue
-            elif kind == EPOCH_STEPPED:
-                levels.append([])
-            else:
-                raise NotImplementedError(
-                    'Unsupported epoch type: {epoch_types(kind)}')
-
-        # Gather steps
-        levels = tuple(levels)
-        for i in range(self._sweeps_per_run):
-            j = 0
-            for e in einfo(ida):
-                if e['type'] == EPOCH_STEPPED:
-                    levels[j].append(e['init_level'] + e['level_inc'] * i)
-                    j += 1
-        return levels
-
-    def filename(self):
-        """ Returns this AbfFile's filename. """
-        return self._filename
-
-    def info(self, show_header=False):
-        """
-        Returns a multi-line string with meta data about this file.
-
-        The optional argument ``show_header`` can be used to add the full
-        header contents to the output.
-        """
-        out = []
-
-        # File info
-        if self._is_protocol_file:
-            out.append(f'Axon Protocol File: {self._filename}')
-        else:
-            out.append(f'Axon Binary File: {self._filename}')
-        out.append(f'ABF Format version {self._version_str}')
-        out.append(f'Recorded on: {self._datetime}')
-
-        # AProtocol info
-        out.append(
-            f'Acquisition mode: {self._mode}: {acquisition_modes[self._mode]}')
-        if self._number_of_trials:
-            out.append(
-                f'Protocol set for {self._number_of_trials} trials,'
-                f' spaced {self._trial_start_to_start}s apart.')
-            out.append(
-                f'    with {self._runs_per_trial} runs per trial,'
-                f' spaced {self._run_start_to_start}s apart.')
-            out.append(
-                f'     and {self._sweeps_per_run} sweeps per run,'
-                f' spaced {self._sweep_start_to_start}s apart.')
-        else:   # pragma: no cover
-            out.append('Protocol data could not be determined.')
-        out.append(f'Sampling rate: {self._rate} Hz')
-
-        # Channel info
-        if self._sweeps:
-            # A/D channels
-            for i, c in enumerate(self._sweeps[0]._channels[:self._nADC]):
-                out.append(f'A/D Channel {i}: "{c._name}"')
-                if c._type:  # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    out.append(f'  Type: {type_mode_names[c._type]}')
-                out.append(f'  Unit: {c._unit.strip()}')
-                if c._lopass:
-                    out.append(f'  Low-pass filter: {c._lopass} Hz')
-                if c._cm:
-                    out.append(f'  Cm (telegraphed): {c._cm} pF')
-                if c._rs:   # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    out.append(f'  Rs (telegraphed): {c._rs}')
-            # Reconstructed D/A channels
-            for i, c in enumerate(self._sweeps[0]._channels[self._nADC:]):
-                out.append(f'D/A Channel {i}: "{c._name}"')
-                if c._type:  # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    out.append('  Type: {type_mode_names[c._type]}')
-                out.append(f'  Unit: {c._unit.strip()}')
-
-        # Show full header info
-        if show_header:
-            if self.strings:
-                dict_to_string(out, 'Strings', {'strings': self.strings})
-            dict_to_string(out, 'file header', self._header)
-
-        return '\n'.join(out)
-
-    def log(self, join_sweeps=False, use_names=False, channels=None):
-        # Docstring in SweepSource
+        # Read the file header as an ordered dictionary
+        self._header = None
+        self._strings = None
+        self._version = None
+        self._version_str = None
+        self._read_1_header()
 
-        # Select channels
-        if channels is None:
-            channels = list(range(self._nc))
-        else:
-            channels = [self._channel_id(i) for i in channels]
+        # Read the time of recording
+        self._datetime = None
+        self._read_2_time_of_recording()
 
-        # Create log
-        log = myokit.DataLog()
+        # Read the protocol information
+        self._n_adc = None
+        self._n_dac = None
 
-        # No data? Then return. Do this after channel check to tell user about
-        # invalid channel ids.
-        ns = len(self._sweeps)
-        if ns == 0:  # pragma: no cover
-            return log
+        self._rate = None
+        self._mode = None
 
-        # Get channel names
-        if use_names:
-            channel_names = [c.name() for c in self._sweeps[0]]
-        else:
-            channel_names = [f'{i}.channel' for i in range(self._nc)]
-        channel_names = [channel_names[i] for i in channels]
+        self._number_of_trials = None
+        self._trial_start_to_start = None
+        self._runs_per_trial = None
+        self._run_start_to_start = None
+        self._sweeps_per_run = None
+        self._sweep_start_to_start = None
+        self._samples_per_channel = None
 
-        # Populate log
-        if join_sweeps:
-            # Join sweeps
+        # To be able to treat v1 and v2 slightly more easily, we define 3
+        # functions read epoch info from episodic stimulation protocols
+        self._epoch_functions = None
 
-            # Set time
-            t = np.array(self._sweeps[0][0].times())
-            time = [t + i * self._sweep_start_to_start for i in range(ns)]
-            log['time'] = np.concatenate(time)
+        # Not all D/A channels can be converted, so we maintain an array with
+        # the original indices of the channels in the da_sweeps. (Note that
+        # this differs again from the "index" labels for user display, which
+        # are stored in the channels themselves).
+        self._dac_indices = []
 
-            # Set channels
-            for channel, name in zip(channels, channel_names):
-                log[name] = np.concatenate(
-                    [sweep[channel].values() for sweep in self._sweeps])
+        # Read protocol information and create empty sweep objects
+        # Sweeps contain both A/D channels and D/A reconstructions. Some files
+        # will have A/D but no (or no supported) D/A. Conversely protocol files
+        # will have D/A only. So all in one sweep is easiest.
+        self._sweeps = None
 
-        else:
-            # Return individual sweeps
-            log['time'] = np.array(self._sweeps[0][0].times())
-            for i, sweep in enumerate(self._sweeps):
-                for channel, name in zip(channels, channel_names):
-                    log[name, i] = sweep[channel].values()
+        self._read_3_protocol_information()
 
-        log.set_time_key('time')
-        return log
+        # Read and calculate conversion factors for integer data in ADC
+        self._adc_factors = None
+        self._adc_offsets = None
+        self._read_4_ad_conversion_factors()
 
-    def matplotlib_figure(self):
-        """ Creates and returns a matplotlib figure with this file's data. """
-        import matplotlib.pyplot as plt
-        f = plt.figure()
-        plt.suptitle(self.filename())
+        # Read the A/D channel data and add it to the sweeps
+        if self._n_adc:
+            self._read_5_ad_data()
 
-        # Plot AD channels
-        ax = plt.subplot(2, 1, 1)
-        ax.set_title('Measured data')
-        times = None
-        for sweep in self._sweeps:
-            for channel in sweep[:self._nADC]:
-                if times is None:
-                    times = channel.times()
-                plt.plot(times, channel.values())
+        # Reconstruct D/A signals as additional channels and add to the sweeps
+        self._read_6_da_reconstructions()
 
-        # Plot DA channels
-        n = self._nDAC
-        ax = [plt.subplot(2, n, n + 1 + i) for i in range(n)]
+        # Copy channel names and units, for easier SweepSource implementation
+        self._ad_names = {}
+        self._da_names = {}
+        self._ad_units = []
+        self._da_units = []
         for sweep in self._sweeps:
-            for i, channel in enumerate(sweep[self._nADC:]):
-                ax[i].set_title(channel.name())
-                ax[i].plot(times, channel.values())
-
-        return f
-
-    def path(self):
-        """ Returns the path to the underlying ABF file. """
-        return self._filepath
-
-    def protocol(self, channel_id=None, ms=True):
-        """
-        Returns the embedded protocol from the D/A channel identified by
-        ``channel_id`` as a :class:`myokit.Protocol`.
-
-        Only episodic stimulation is supported, without "user lists".
-
-        By default all times are converted to milliseconds, but this can be
-        disabled by setting ``ms=False``.
-
-        Arguments:
-
-        ``channel``
-            The channel index to return (as an integer).
-        ``ms``
-            Set to ``False`` to disable conversion to milliseconds.
-
-        """
-        # Only episodic stimulation is supported.
-        if self._mode != ACMODE_EPISODIC_STIMULATION:  # pragma: no cover
-            raise NotImplementedError(
-                'Only episodic stimulation is supported.')
-
-        # Get D/A channel index
-        if channel_id is None:
-            ida = self._nADC
-        else:
-            ida = self._channel_id(channel_id, must_be_da=True) - self._nADC
-
-        # User lists are not supported
-        if self._version < 2:  # pragma: no cover
-            if self._header['nULEnable'][ida]:
-                raise NotImplementedError('User lists are not supported.')
-        else:   # pragma: no cover
-            for userlist in self._header['listUserListInfo']:
-                en1 = 'nULEnable' in userlist and userlist['nULEnable']
-                en2 = 'nConditEnable' in userlist and userlist['nConditEnable']
-                if en1 or en2:
-                    raise NotImplementedError('User lists are not supported.')
-
-        # Create protocol
-        p = myokit.Protocol()
-
-        dinfo, einfo_exists, einfo = self._epoch_functions
-        start = 0
-        next_start = 0
-        f = 1e3 if ms else 1
-        for i_sweep in range(self._sweeps_per_run):
-            for e in einfo(ida):
-                kind = e['type']
-
-                if kind not in epoch_types:
-                    raise NotImplementedError(f'Unknown epoch type: {kind}')
-
-                if kind == EPOCH_DISABLED:
-                    continue
-
-                elif kind == EPOCH_STEPPED:
-                    # Event at step
-                    dur = f * e['init_duration'] / self._rate
-                    inc = f * e['duration_inc'] / self._rate
-                    e_level = e['init_level'] + e['level_inc'] * i_sweep
-                    e_start = start
-                    e_length = dur + i_sweep * inc
-                    p.schedule(e_level, e_start, e_length)
-                    start += e_length
-
-                else:   # pragma: no cover
-                    raise NotImplementedError(
-                        'Usupported epoch type: ' + epoch_types(kind))
-
-            # Event at holding potential
-            next_start += f * self._sweep_start_to_start
-            e_level = dinfo(ida, 'fDACHoldingLevel')
-            e_start = start
-            e_length = next_start - start
-            p.schedule(e_level, e_start, e_length)
-            start = next_start
-
-        return p
-
-    def _read_datetime(self):
-        """ Reads the date/time this file was recorded """
-        # Get date and time
-        if self._version < 2:
-            t1 = str(self._header['lFileStartDate'])
-            t2 = float(self._header['lFileStartTime'])
-        else:
-            t1 = str(self._header['uFileStartDate'])
-            t2 = float(self._header['uFileStartTimeMS']) / 1000
-
-        YY = int(t1[0:4])
-        MM = int(t1[4:6])
-        DD = int(t1[6:8])
-        hh = int(t2 / 3600)
-        mm = int((t2 - hh * 3600) / 60)
-        ss = t2 - hh * 3600 - mm * 60
-        ms = int((ss % 1) * 1e6)
-        ss = int(ss)
-
-        return datetime.datetime(YY, MM, DD, hh, mm, ss, ms)
+            for i, channel in enumerate(sweep[:self._n_adc]):
+                self._ad_names[channel.name()] = i
+                self._ad_units.append(channel.unit())
+            for i, channel in enumerate(sweep[self._n_adc:]):
+                self._da_names[channel.name()] = i
+                self._da_units.append(channel.unit())
+            break
 
-    def _read_header(self):
-        """ Reads the file's header. """
+    def _read_1_header(self):
+        """ Read the file header. """
 
         def read_f(f, form, offset=None):
             """ Read and unpack a file section in the format ``form``. """
@@ -718,7 +330,7 @@ def ups(val):
             values = [0] * len(val)
             for i, v in enumerate(val):
                 if isinstance(v, bytes):
-                    v = v.decode(_ENC)
+                    v = v.decode(_ENC).strip()
                     # Handle long \x00 lists
                     if v and ord(v[0]) == 0:
                         return None
@@ -743,8 +355,8 @@ def ups(val):
                 raise NotImplementedError('Unknown ABF Format "{sig}".')
 
             # Gather header fields
-            header = OrderedDict()
-            for key, offset, form in headerFields[version]:
+            header = {}
+            for key, offset, form in HEADER_FIELDS[version]:
                 header[key] = ups(read_f(f, form, offset))
 
             # Get uniform file version number
@@ -770,21 +382,21 @@ def ups(val):
                 for i in range(header['lNumTagEntries']):  # pragma: no cover
                     # Cover pragma: Don't have appropriate test file
                     f.seek(header['lTagSectionPtr'] + i * 64)
-                    tag = OrderedDict()
-                    for key, form in TagInfoDescription:
+                    tag = {}
+                    for key, form in ABF2_TAG_INFO_DESCRIPTION:
                         tag[key] = ups(read_f(f, form))
                     tags.append(tag)
                 header['tags'] = tags
-                self.strings = []
+                self._strings = []
 
             else:
 
                 # Version 2
                 # Find location of file sections
-                sections = OrderedDict()
-                for i, s in enumerate(abf2FileSections):
+                sections = {}
+                for i, s in enumerate(ABF2_FILE_SECTIONS):
                     index, data, length = read_f(f, 'IIl', 76 + i * 16)
-                    sections[s] = OrderedDict()
+                    sections[s] = {}
                     sections[s]['index'] = index
                     sections[s]['data'] = data
                     sections[s]['length'] = length
@@ -806,8 +418,8 @@ def ups(val):
 
                 # C-style string termination
                 strings = strings.split(b'\x00')
-                strings = [s.decode(_ENC) for s in strings]
-                self.strings = strings
+                strings = [s.decode(_ENC).strip() for s in strings]
+                self._strings = strings
 
                 # Read tag section
                 tags = []
@@ -816,17 +428,17 @@ def ups(val):
                 for i in range(sections['Tag']['length']):  # pragma: no cover
                     # Cover pragma: Don't have appropriate test file
                     f.seek(offs + i * size)
-                    tag = OrderedDict()
-                    for key, form in TagInfoDescription:
+                    tag = {}
+                    for key, form in ABF2_TAG_INFO_DESCRIPTION:
                         tag[key] = ups(read_f(f, form))
                     tags.append(tag)
                 header['tags'] = tags
 
                 # Read protocol section
-                protocol = OrderedDict()
+                protocol = {}
                 offs = sections['Protocol']['index'] * BLOCKSIZE
                 f.seek(offs)
-                for key, form in protocolFields:
+                for key, form in ABF2_PROTOCOL_FIELDS:
                     protocol[key] = ups(read_f(f, form))
                 header['protocol'] = protocol
 
@@ -835,9 +447,9 @@ def ups(val):
                 offs = sections['ADC']['index'] * BLOCKSIZE
                 size = sections['ADC']['data']
                 for i in range(sections['ADC']['length']):
-                    ADC = OrderedDict()
+                    ADC = {}
                     f.seek(offs + i * size)
-                    for key, form in ADCFields:
+                    for key, form in ABF2_ADC_FIELDS:
                         ADC[key] = ups(read_f(f, form))
                     # Get channel name and unit
                     ADC['ADCChNames'] = (
@@ -852,8 +464,8 @@ def ups(val):
                 size = sections['DAC']['data']
                 for i in range(sections['DAC']['length']):
                     f.seek(offs + size * i)
-                    DAC = OrderedDict()
-                    for key, form in DACFields:
+                    DAC = {}
+                    for key, form in ABF2_DAC_FIELDS:
                         DAC[key] = ups(read_f(f, form))
                     DAC['sDACChannelName'] = \
                         strings[DAC['lDACChannelNameIndex'] - 1]
@@ -863,62 +475,111 @@ def ups(val):
                 header['listDACInfo'] = dac
 
                 # Read UserList section
-                userlists = []
+                user_lists = []
+                offs = sections['UserList']['index'] * BLOCKSIZE
+                size = sections['UserList']['data']
                 r = range(sections['UserList']['length'])
                 for i in r:  # pragma: no cover
                     # Cover pragma: User lists are not supported
                     f.seek(offs + size * i)
-                    UserList = OrderedDict()
-                    for key, form in UserListFields:
-                        UserList[key] = ups(read_f(f, form))
-                    userlists.append(DAC)
-                header['listUserListInfo'] = userlists
+                    user_list = {}
+                    for key, form in ABF2_USER_LIST_FIELDS:
+                        user_list[key] = ups(read_f(f, form))
+                    user_lists.append(user_list)
+                header['listUserListInfo'] = user_lists
 
                 # Read epoch-per-DAC section
                 # The resulting OrderedDict has the following structure:
                 #  - the first index is the DAC number
                 #  - the second index is the epoch number
-                header['epochInfoPerDAC'] = OrderedDict()
+                header['epochInfoPerDAC'] = {}
                 offs = sections['EpochPerDAC']['index'] * BLOCKSIZE
                 size = sections['EpochPerDAC']['data']
-                info = OrderedDict()
+                info = {}
                 for i in range(sections['EpochPerDAC']['length']):
                     f.seek(offs + size * i)
-                    einf = OrderedDict()
-                    for key, form in EpochInfoPerDACFields:
+                    einf = {}
+                    for key, form in ABF2_EPOCH_INFO_PER_DAC_FIELD:
                         einf[key] = ups(read_f(f, form))
                     DACNum = einf['nDACNum']
                     EpochNum = einf['nEpochNum']
                     if DACNum not in info:
-                        info[DACNum] = OrderedDict()
+                        info[DACNum] = {}
                     info[DACNum][EpochNum] = einf
                 header['epochInfoPerDAC'] = info
 
-        return header
+        self._header = header
 
-    def _read_1_protocol(self):
-        """
-        Reads the protocols stored in the ABF file and converts it to an analog
-        signal if possible.
+    def _read_2_time_of_recording(self):
+        """ Read and process the time when this file was recorded. """
+
+        if self._version < 2:
+            t1 = str(self._header['lFileStartDate'])
+            t2 = float(self._header['lFileStartTime'])
+        else:
+            t1 = str(self._header['uFileStartDate'])
+            t2 = float(self._header['uFileStartTimeMS']) / 1000
+
+        YY = int(t1[0:4])
+        MM = int(t1[4:6])
+        DD = int(t1[6:8])
+        hh = int(t2 / 3600)
+        mm = int((t2 - hh * 3600) / 60)
+        ss = t2 - hh * 3600 - mm * 60
+        ms = int((ss % 1) * 1e6)
+        ss = int(ss)
 
-        Only works for episodic stimulation, without any user lists.
+        self._datetime = datetime.datetime(YY, MM, DD, hh, mm, ss, ms)
 
-        The resulting analog signal has the same size as the recorded signals,
-        so not the full length of the protocol!
+    def _read_3_protocol_information(self):
         """
+        Reads the header fields detailing the number of runs, sweeps, and the
+        type of protocol used. Create empty sweeps.
+        """
+        h = self._header
+
+        # Number of channels, sampling rate (Hz) and acquisition mode
+        # Note: Number of A/D channels will be set to 0 if this is a
+        #       protocol-only file
+        # Note: Number of D/A channels will be adjusted after checking support
+        if self._version < 2:
+            # In (newer versions of) version 1.x, only 2 D/A channels have
+            # full "waveform" support. There are still 4 D/A channels but I
+            # don't understand what the other 2 do.
+            # 1.x versions only seem to have 1 DAC channel, but this is not
+            # supported here.
+            self._n_adc = int(h['nADCNumChannels'])
+            self._n_dac = min(len(h['sDACChannelName']), 2)
+            self._rate = 1e6 / (h['fADCSampleInterval'] * self._n_adc)
+            self._mode = h['nOperationMode']
+        else:
+            # In version 2, there are up to 8 "waveform" D/A channels
+
+            self._n_adc = int(h['sections']['ADC']['length'])
+            self._n_dac = int(h['sections']['DAC']['length'])
+            self._rate = 1e6 / h['protocol']['fADCSequenceInterval']
+            self._mode = h['protocol']['nOperationMode']
+
+        if self._mode not in acquisition_modes:  # pragma: no cover
+            raise NotImplementedError(f'Unknown acquisition mode: {mode}')
+
+        # Protocol files don't have A/D channels by definition
+        if self._is_protocol_file:
+            self._n_adc = 0
+
         # Only episodic stimulation is supported.
         if self._mode != ACMODE_EPISODIC_STIMULATION:  # pragma: no cover
             warnings.warn(
                 'Unsupported acquisition method '
                 + acquisition_modes[self._mode] + '; unable to read D/A'
                 ' channels.')
-            self._nDAC = 0
-            return
 
-        # Start reading
-        h = self._header
+            # Remaining code is all about reading D/A info for episodic
+            # stimulation, so return
+            self._n_dac = 0
+            return
 
-        # Step 1: Gather information about the protocol
+        # Gather protocol information
         if self._version < 2:
 
             # Before version 2: Sections are fixed length, locations absolute
@@ -931,17 +592,24 @@ def _read_1_protocol(self):
 
             # Number of samples in a channel for each sweep
             # (Only works for fixed-length, high-speed-osc or episodic)
-            nSam = h['lNumSamplesPerEpisode'] // h['nADCNumChannels']
+            self._samples_per_channel = \
+                h['lNumSamplesPerEpisode'] // h['nADCNumChannels']
 
             def dinfo(index, name):
+                """ Return DAC channel info, ABF1 version. """
                 return h[name][index]
 
             def einfo_exists(index):
-                return True
+                """ Check that epoch info exists for a DAC, ABF1 version. """
+                # Fields always exist for 2 channels, not always set.
+                # But not useful to look at unset ones, so using n_dac instead
+                # of hardcoded 2!
+                return 0 <= index < self._n_dac
 
             def einfo(index):
+                """ Return epoch info for a DAC, ABF1 version. """
                 lo = index * 8
-                hi = index + 8
+                hi = lo + 8
                 for i in range(lo, hi):
                     yield {
                         'type': h['nEpochType'][i],
@@ -967,16 +635,20 @@ def einfo(index):
             self._sweep_start_to_start = p['fSweepStartToStart']
 
             # Number of samples in a channel in a single sweep
-            nSam = p['lNumSamplesPerEpisode'] // h['sections']['ADC']['length']
+            self._samples_per_channel = \
+                p['lNumSamplesPerEpisode'] // h['sections']['ADC']['length']
 
             # Compatibility functions
             def dinfo(index, name):
+                """ Return DAC info, ABF2 version. """
                 return h['listDACInfo'][index][name]
 
             def einfo_exists(index):
+                """ Check that epoch info exists for a DAC, ABF2 version. """
                 return index in h['epochInfoPerDAC']
 
             def einfo(index):
+                """ Return epoch info for a DAC, ABF2 version. """
                 for e in h['epochInfoPerDAC'][index].values():
                     yield {
                         'type': e['nEpochType'],
@@ -989,165 +661,192 @@ def einfo(index):
 
         # If sweepStartToStart == 0, we set it to the duration of a sweep
         if self._sweep_start_to_start == 0:    # pragma: no cover
-            self._sweep_start_to_start = nSam / self._rate
+            self._sweep_start_to_start = self._samples_per_channel / self._rate
 
-        # Step 2: Generate analog signals corresponding to the waveforms
-        # suggested by the 'epochs' in the protocol
-
-        # User lists are not supported
+        # Create empty sweeps
+        n = h['lActualSweeps']
+        if self._is_protocol_file:
+            n = self._sweeps_per_run
+        self._sweeps = [Sweep() for i in range(n)]
+
+        # User lists are not supported for D/A reconstruction
+        # I haven't been able to figure out how you see if a user list is
+        # being used, or which channel is using it. There is an 'enable' field
+        # but that's been 0 in files that definitely used a UserList...
+        # So for now not reading ANY DAC if a userlist even exists.
         user_lists = False
         if self._version < 2:   # pragma: no cover
             user_lists = any(self._header['nULEnable'])
-        else:  # pragma: no cover
-            for userlist in self._header['listUserListInfo']:
-                if 'nULEnable' in userlist and userlist['nULEnable']:
-                    user_lists = True
-                    break
-                if 'nConditEnable' in userlist and userlist['nConditEnable']:
-                    user_lists = True
-                    break
-        if user_lists:  # pragma: no cover
+        else:                   # pragma: no cover
+            user_lists = len(self._header['listUserListInfo']) > 0
+        if user_lists:          # pragma: no cover
             warnings.warn(
                 'Unsupported acquisition method: episodic with user lists;'
                 ' unable to read D/A channels.')
-            self._nDAC = 0
+            self._n_dac = 0
             return
 
-        # Number of ADC channels
-        nADC = self._nADC
+        # Get indices of enabled and supported DAC reconstructions
+        supported = {EPOCH_DISABLED, EPOCH_STEPPED}
+        for i_dac in range(self._n_dac):
+            if einfo_exists(i_dac):
+                i = einfo(i_dac)
+                use = False
+
+                # Check for unsupported features (or disabled waveforms/epochs)
+                # Version 1 files can only have two waveform channels
+                if self._version < 2 and i_dac > 1:  # pragma: no cover
+                    source = DAC_DISABLED
+                else:
+                    source = dinfo(i_dac, 'nWaveformSource')
+                if source == DAC_EPOCHTABLEWAVEFORM:
+                    # Any epoch types besides disabled/stepped? Then don't use
+                    # Also don't use if exclusively disabled
+                    for e in i:
+                        t = e['type']
+                        if t == EPOCH_STEPPED:
+                            use = True
+                        elif t != EPOCH_DISABLED:  # pragma: no cover
+                            use = False
+                            warnings.warn(
+                                f'Unsupported epoch type: {epoch_types(t)}')
+                            break
+                elif source == DAC_DACFILEWAVEFORM:  # pragma: no cover
+                    # Stimulus file? Then don't use
+                    warnings.warn('Stimulus file D/A channel not supported.')
 
-        # Number of used DAC channels
-        channels = []
-        for iDAC in range(self._nDAC):
-            if einfo_exists(iDAC):
-                if any(e['type'] != EPOCH_DISABLED for e in einfo(iDAC)):
-                    channels.append(iDAC)
-        self._nDAC = nDAC = len(channels)
+                if use:
+                    self._dac_indices.append(i_dac)
 
-        # Number of sweeps
-        if self._is_protocol_file:
-            n_sweeps = self._sweeps_per_run
-        else:
-            n_sweeps = h['lActualSweeps']
+        # Set true number of D/A outputs
+        self._n_dac = len(self._dac_indices)
 
-        # Read
-        start = 0
-        for i_sweep in range(n_sweeps):
-            sweep = Sweep(nADC + nDAC)
+    def _read_4_ad_conversion_factors(self):
+        """ Calculate the factors to convert any integer data to float. """
+        self._adc_factors = []
+        self._adc_offsets = []
+        h = self._header
+        if self._version < 2:
+            for i in range(self._n_adc):
+                j = h['nADCSamplingSeq'][i]
 
-            # Create channels for this sweep
-            i_channel = nADC
-            for iDAC in channels:
+                # Multiplier
+                f = (
+                    h['fInstrumentScaleFactor'][j]
+                    * h['fADCProgrammableGain'][j]
+                    * h['lADCResolution']
+                    / h['fADCRange'])
 
-                # Convert, where possible
-                c = Channel(self)
-                c._name = dinfo(iDAC, 'sDACChannelName').strip()
-                c._unit = dinfo(iDAC, 'sDACChannelUnits').strip()
-                if self._version < 2:
-                    c._numb = iDAC
-                else:
-                    c._numb = int(dinfo(iDAC, 'lDACChannelNameIndex'))
-                c._data = np.ones(nSam) * dinfo(iDAC, 'fDACHoldingLevel')
-                c._rate = self._rate
-                c._start = start
-                sweep[i_channel] = c
-                i_channel += 1
+                # Signal conditioner used?
+                if h['nSignalType'] != 0:   # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    f *= h['fSignalGain'][j]
 
-                # Save last sample index
-                i_last = int(nSam * 15625 / 1e6)  # TODO: What's this?
+                # Additional gain?
+                if h['nTelegraphEnable'][j]:
+                    f *= h['fTelegraphAdditGain'][j]
 
-                # For each 'epoch' in the stimulation signal
-                for e in einfo(iDAC):
-                    kind = e['type']
-                    if kind == EPOCH_DISABLED:
+                # Set final gain factor
+                self._adc_factors.append(1 / f)
 
-                        continue
+                # Shift
+                s = h['fInstrumentOffset'][j]
 
-                    elif kind == EPOCH_STEPPED:
+                # Signal conditioner used?
+                if h['nSignalType'] != 0:   # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    s -= h['fSignalOffset'][j]
 
-                        dur = e['init_duration']
-                        inc = e['duration_inc']
-                        i1 = i_last
-                        i2 = i_last + dur + i_sweep * inc
-                        if i2 > nSam:
-                            # The protocol may extend beyond the number of
-                            # samples in the recording
-                            i2 = nSam
-                        level = e['init_level'] + e['level_inc'] * i_sweep
-                        c._data[i1:i2] = level * np.ones(len(range(i2 - i1)))
-                        i_last += dur
-                        if i_last > nSam:
-                            # The protocol may extend beyond the number of
-                            # samples in the recording
-                            break
+                # Set final offset
+                self._adc_offsets.append(s)
+
+        else:
+
+            a = h['listADCInfo']
+            p = h['protocol']
+            for i in range(self._n_adc):
+                # Multiplier
+                f = (
+                    a[i]['fInstrumentScaleFactor']
+                    * a[i]['fADCProgrammableGain']
+                    * p['lADCResolution']
+                    / p['fADCRange'])
 
-                    else:  # pragma: no cover
+                # Signal conditioner used?
+                if h.get('nSignalType', 0) != 0:  # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    f *= a[i]['fSignalGain']
 
-                        warnings.warn(
-                            'Unsupported epoch type: ' + epoch_types(kind))
-                        continue
+                # Additional gain?
+                if a[i]['nTelegraphEnable']:
+                    f *= a[i]['fTelegraphAdditGain']
 
-            self._sweeps.append(sweep)
-            start += self._sweep_start_to_start
+                # Set final gain factor
+                self._adc_factors.append(1 / f)
+
+                # Shift
+                s = a[i]['fInstrumentOffset']
+
+                # Signal conditioner used?
+                if h.get('nSignalType', 0) != 0:  # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    s -= a[i]['fSignalOffset']
+
+                # Set final offset
+                self._adc_offsets.append(s)
 
-    def _read_2_data(self):
+    def _read_5_ad_data(self):
         """ Reads the A/D data and appends it to the list of sweeps. """
 
-        header = self._header
-        version = self._version
-        nc = self._nADC
+        h = self._header
 
         # Sampling rate is constant for all sweeps and channels
         # TODO: This won't work for 2-rate protocols
         rate = self._rate
 
         # Get binary integer format
-        dt = np.dtype('i2') if header['nDataFormat'] == 0 else np.dtype('f4')
+        dt = np.dtype('i2') if h['nDataFormat'] == 0 else np.dtype('f4')
 
         # Get number of channels, create a numpy memory map
-        if version < 2:
+        if self._version < 2:
             # Old files, get info from fields stored directly in header
-            o = header['lDataSectionPtr'] * BLOCKSIZE \
-                + header['nNumPointsIgnored'] * dt.itemsize
-            n = header['lActualAcqLength']
+            o = h['lDataSectionPtr'] * BLOCKSIZE \
+                + h['nNumPointsIgnored'] * dt.itemsize
+            n = h['lActualAcqLength']
         else:
             # New files, get info from appropriate header section
-            o = header['sections']['Data']['index'] * BLOCKSIZE
-            n = header['sections']['Data']['length']
+            o = h['sections']['Data']['index'] * BLOCKSIZE
+            n = h['sections']['Data']['length']
         data = np.memmap(self._filepath, dt, 'r', shape=(n,), offset=o)
 
         # Load list of sweeps (Sweeps are called 'episodes' in ABF < 2)
-        if version < 2:
-            n = header['lSynchArraySize']
-            o = header['lSynchArrayPtr'] * BLOCKSIZE
+        if self._version < 2:
+            n = h['lSynchArraySize']
+            o = h['lSynchArrayPtr'] * BLOCKSIZE
         else:
-            n = header['sections']['SynchArray']['length']
-            o = header['sections']['SynchArray']['index'] * BLOCKSIZE
+            n = h['sections']['SynchArray']['length']
+            o = h['sections']['SynchArray']['index'] * BLOCKSIZE
         if n > 0:
             dt = [(str('offset'), str('i4')), (str('len'), str('i4'))]
-            sdata = np.memmap(self._filepath, dt, 'r', shape=(n,), offset=o)
+            sweep_data = np.memmap(
+                self._filepath, dt, 'r', shape=(n,), offset=o)
         else:   # pragma: no cover
             # Cover pragma: Don't have appropriate test file
-            sdata = np.empty((1), dt)
-            sdata[0]['len'] = data.size
-            sdata[0]['offset'] = 0
+            sweep_data = np.empty((1), dt)
+            sweep_data[0]['len'] = data.size
+            sweep_data[0]['offset'] = 0
 
         # Number of sweeps must equal n
-        if n != header['lActualSweeps']:
+        if n != h['lActualSweeps']:
             raise NotImplementedError(
                 'Unable to read file with different sizes per sweep.')
 
-        # Sweeps should have been created by read_protocol, but if the protocol
-        # type was unsupported they might not have been
-        if len(self._sweeps) == 0:  # pragma: no cover
-            self._sweeps = [Sweep(nc) for i in range(n)]
-
-        # Time-offset at start of sweep
-        start = sdata[0]['offset'] / rate
+        # Time-offset at start of first sweep
+        start = sweep_data[0]['offset'] / rate
 
         # Get data
         pos = 0
-        for sweep, sdat in zip(self._sweeps, sdata):
+        for i_sweep, sdat in enumerate(sweep_data):
 
             # Get the number of data points
             size = sdat['len']
@@ -1155,42 +854,51 @@ def _read_2_data(self):
             # Calculate the correct size for variable-length event mode
             if self._mode == ACMODE_VARIABLE_LENGTH_EVENTS:  # pragma: no cover
                 # Cover pragma: Only episodic stimulus is supported.
-                if version < 2:
-                    f = float(header['fSynchTimeUnit'])
+                if self._version < 2:
+                    f = float(h['fSynchTimeUnit'])
                 else:
-                    f = float(header['protocol']['fSynchTimeUnit'])
+                    f = float(h['protocol']['fSynchTimeUnit'])
                 if f != 0:
                     size /= f
 
             # Get a memory map to the relevant part of the data
             part = data[pos: pos + size]
             pos += size
-            part = part.reshape((part.size // nc, nc)).astype('f')
+            part = part.reshape(
+                (part.size // self._n_adc, self._n_adc)).astype('f')
 
             # If needed, reformat the integers
-            if header['nDataFormat'] == 0:
+            if h['nDataFormat'] == 0:
                 # Data given as integers? Convert to floating point
-                for i in range(nc):
+
+                for i in range(self._n_adc):
                     part[:, i] *= self._adc_factors[i]
                     part[:, i] += self._adc_offsets[i]
 
-            # Create channel
+            # Get start in other modes
             if self._mode != ACMODE_EPISODIC_STIMULATION:  # pragma: no cover
                 # All modes except episodic stimulation
-                start = sdata[j]['offset'] / rate
+                start = data['offset'] / rate
 
-            for i in range(nc):
+            # Create and populate sweep
+            sweep = self._sweeps[i_sweep]
+            for i in range(self._n_adc):
                 c = Channel(self)
-                c._data = part[:, i]
-                if version < 2:
-                    c._name = str(header['sADCChannelName'][i]).strip()
-                    c._unit = str(header['sADCUnits'][i]).strip()
-                    c._numb = int(header['nADCPtoLChannelMap'][i])
+                c._data = part[:, i]    # Actually store the data
+                c._rate = rate
+                c._start = start
+
+                if self._version < 2:
+                    j = h['nADCSamplingSeq'][i]
+
+                    c._name = h['sADCChannelName'][j]
+                    c._index = int(h['nADCPtoLChannelMap'][j])
+                    c._unit = self._unit(h['sADCUnits'][j])
 
                     # Get telegraphed info
                     def get(field):
                         try:
-                            return float(header[field][i])
+                            return float(h[field][j])
                         except KeyError:
                             return None
 
@@ -1201,7 +909,7 @@ def get(field):
                         c._lopass = get('fTelegraphFilter')
 
                     # Updated low-pass cutoff
-                    if header['nSignalType'] != 0:  # pragma: no cover
+                    if h['nSignalType'] != 0:  # pragma: no cover
                         # Cover pragma: Don't have appropriate test file
                         # If a signal conditioner is used, the cutoff frequency
                         # is an undescribed "complex function" of both low-pass
@@ -1209,122 +917,433 @@ def get(field):
                         c._lopass = None
 
                 else:
-
-                    c._name = str(header['listADCInfo'][i]['ADCChNames']
-                                  ).strip()
-                    c._unit = str(header['listADCInfo'][i]['ADCChUnits']
-                                  ).strip()
-                    c._numb = int(header['listADCInfo'][i]['nADCNum'])
+                    c._name = h['listADCInfo'][i]['ADCChNames']
+                    c._index = int(h['listADCInfo'][i]['nADCNum'])
+                    c._unit = self._unit(h['listADCInfo'][i]['ADCChUnits'])
 
                     # Get telegraphed info
-                    if header['listADCInfo'][i]['nTelegraphEnable']:
-                        c._type = int(
-                            header['listADCInfo'][i]['nTelegraphMode'])
+                    if h['listADCInfo'][i]['nTelegraphEnable']:
+                        c._type = int(h['listADCInfo'][i]['nTelegraphMode'])
                         c._cm = float(
-                            header['listADCInfo'][i]['fTelegraphMembraneCap'])
+                            h['listADCInfo'][i]['fTelegraphMembraneCap'])
                         c._rs = float(
-                            header['listADCInfo'][i][
-                                'fTelegraphAccessResistance'])
+                            h['listADCInfo'][i]['fTelegraphAccessResistance'])
                         c._lopass = float(
-                            header['listADCInfo'][i]['fTelegraphFilter'])
+                            h['listADCInfo'][i]['fTelegraphFilter'])
 
                     # Updated low-pass cutoff
-                    if 'nSignalType' in header['protocol']:  # pragma: no cover
+                    if 'nSignalType' in h['protocol']:  # pragma: no cover
                         # Cover pragma: Don't have appropriate test file
-                        if header['protocol']['nSignalType'] != 0:
+                        if h['protocol']['nSignalType'] != 0:
                             # If a signal conditioner is used, the cutoff
                             # frequency is an undescribed "complex function" of
                             # both low-pass settings...
                             c._lopass = None
-                c._rate = rate
-                c._start = start
-                sweep[i] = c
+
+                sweep._channels.append(c)
 
             if self._mode == ACMODE_EPISODIC_STIMULATION:
                 # Increase time according to sweeps in episodic stim. mode
                 start += self._sweep_start_to_start
 
-    def _set_conversion_factors(self):
+    def _read_6_da_reconstructions(self):
         """
-        Calculates the conversion factors to convert integer data from the ABF
-        file to floats.
+        Convert supported D/A waveforms to channels.
+
+        Only works for episodic stimulation, with step protocols and no
+        user lists.
+
+        The resulting analog signal has the same size as the recorded
+        signals, so not always the full length of the protocol!
+
         """
-        self._adc_factors = []
-        self._adc_offsets = []
-        h = self._header
-        if self._version < 2:
-            for i in range(self._nADC):
-                # Multiplier
-                f = (
-                    h['fInstrumentScaleFactor'][i]
-                    * h['fADCProgrammableGain'][i]
-                    * h['lADCResolution']
-                    / h['fADCRange'])
+        dinfo, einfo_exists, einfo = self._epoch_functions
 
-                # Signal conditioner used?
-                if h['nSignalType'] != 0:   # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    f *= h['fSignalGain'][i]
+        ns = self._samples_per_channel
+        start = 0
+        for i_sweep, sweep in enumerate(self._sweeps):
+            for i_dac in self._dac_indices:
 
-                # Additional gain?
-                if h['nTelegraphEnable'][i]:
-                    f *= h['fTelegraphAdditGain'][i]
+                # Create a channel
+                c = Channel(self)
+                c._name = dinfo(i_dac, 'sDACChannelName')
+                if self._version < 2:
+                    c._index = i_dac
+                else:
+                    c._index = int(dinfo(i_dac, 'lDACChannelNameIndex'))
+                c._data = np.ones(ns) * dinfo(i_dac, 'fDACHoldingLevel')
+                c._rate = self._rate
+                c._start = start
+                c._unit = self._unit(dinfo(i_dac, 'sDACChannelUnits'))
+
+                # Find start of first epoch. This is defined as being at t=0
+                # but axon likes to add some samples before the first and after
+                # the last epoch. We can find out the number of samples using
+                # a procedure found in ABF v1's _GetHoldingLength()
+                if self._is_protocol_file:
+                    i2 = 0
+                else:
+                    i2 = ns // 64  # ABFH_HOLDINGFRACTION = 64
+                    i2 -= i2 % self._n_adc
+                    if (i2 < self._n_adc):  # pragma: no cover
+                        i2 = self._n_adc
 
-                # Set fina gain factor
-                self._adc_factors.append(1.0 / f)
+                # For each 'epoch' in the stimulation signal
+                for e in einfo(i_dac):
+                    if e['type'] == EPOCH_STEPPED:
+                        dur = e['init_duration']
+                        inc = e['duration_inc']
+                        i1 = i2
+                        i2 += dur + i_sweep * inc
+                        level = e['init_level'] + e['level_inc'] * i_sweep
+                        c._data[i1:i2] = level * np.ones(len(range(i2 - i1)))
 
-                # Shift
-                s = h['fInstrumentOffset'][i]
+                # Store channel
+                sweep._channels.append(c)
 
-                # Signal conditioner used?
-                if h['nSignalType'] != 0:   # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    s -= h['fSignalOffset'][i]
+            # Update start for next sweep
+            start += self._sweep_start_to_start
 
-                # Set final offset
-                self._adc_offsets.append(s)
+    def __getitem__(self, key):
+        return self._sweeps[key]
+
+    def __iter__(self):
+        return iter(self._sweeps)
+
+    def __len__(self):
+        return len(self._sweeps)
+
+    def _channel_id(self, channel_id):
+        """ Checks an int or str channel id and returns a valid int. """
+        if len(self._sweeps) == 0:  # pragma: no cover
+            raise KeyError(f'Channel {channel_id} not found (empty file).')
+
+        # Handle string
+        if isinstance(channel_id, str):
+            int_id = self._ad_names[channel_id]  # Bubble KeyError to user
+        else:
+            int_id = int(channel_id)    # Propagate TypeError
+            if int_id < 0 or int_id >= self._n_adc:
+                raise IndexError(f'channel_id out of range: {channel_id}')
+
+        return int_id
+
+    def channel(self, channel_id, join_sweeps=False):
+        # Docstring in SweepSource
+        channel_id = self._channel_id(channel_id)
+        time, data = [], []
+        for i, sweep in enumerate(self._sweeps):
+            time.append(sweep[channel_id].times())
+            data.append(sweep[channel_id].values())
+        if join_sweeps:
+            return (np.concatenate(time), np.concatenate(data))
+        return time, data
+
+    def channel_count(self):
+        # Docstring in SweepSource
+        return self._n_adc
+
+    def channel_names(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._ad_names.keys())
+        return list(self._ad_names.keys())[index]
+
+    def channel_units(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._ad_units)
+        return self._ad_units[index]
+
+    def _da_id(self, output_id):
+        """
+        Checks an int or str D/A channel id and returns a valid int.
+
+        Note: The integer here is from 0 to da_count(), so not equal to the
+        channel :meth:`index()` shown in pclamp.
+        """
+        if len(self._sweeps) == 0:  # pragma: no cover
+            raise KeyError(f'D/A output {output_id} not found (empty file).')
 
+        # Handle string
+        if isinstance(output_id, str):
+            int_id = self._da_names[output_id]  # Propagate KeyError to user
         else:
+            int_id = int(output_id)  # Propagate TypeError
+            if int_id < 0 or int_id >= self._n_dac:
+                raise IndexError(f'output_id out of range: {output_id}')
 
-            a = h['listADCInfo']
-            p = h['protocol']
-            for i in range(self._nADC):
-                # Multiplier
-                f = (
-                    a[i]['fInstrumentScaleFactor']
-                    * a[i]['fADCProgrammableGain']
-                    * p['lADCResolution']
-                    / p['fADCRange'])
+        return int_id
 
-                # Signal conditioner used?
-                if 'nSignalType' in h:  # pragma: no cover
-                    # Cover pragma: Don't have appropriate test file
-                    if h['nSignalType'] != 0:
-                        f *= a[i]['fSignalGain']
+    def da(self, output_id, join_sweeps=False):
+        # Docstring in SweepSource
+        channel_id = self._n_adc + self._da_id(output_id)
+        time, data = [], []
+        for i, sweep in enumerate(self._sweeps):
+            time.append(sweep[channel_id].times())
+            data.append(sweep[channel_id].values())
+        if join_sweeps:
+            return (np.concatenate(time), np.concatenate(data))
+        return time, data
 
-                # Additional gain?
-                if a[i]['nTelegraphEnable']:
-                    f *= a[i]['fTelegraphAdditGain']
+    def da_count(self):
+        # Docstring in SweepSource
+        return self._n_dac
 
-                # Set final gain factor
-                self._adc_factors.append(1.0 / f)
+    def da_names(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._da_names.keys())
+        return list(self._da_names.keys())[index]
 
-                # Shift
-                s = a[i]['fInstrumentOffset']
+    def da_protocol(self, output_id=0, tu='ms', vu='mV', cu='pA', n_digits=9,
+                    include_initial_holding=False):
+        """
+        See :meth:`myokit.formats.SweepSource.da_protocol()`.
 
-                # Signal conditioner used?
-                if 'nSignalType' in h:  # pragma: no cover
+        This implementation adds a keyword argument ``include_initial_holding``
+        that lets you switch between the declared protocol (``False``) and the
+        protocol as actually implemented (``True``). In the latter case, a
+        short holding time is added before the first epoch in every sweep.
+        """
+
+        # Check the output id. This also raises an error if no supported D/A
+        # channels are present.
+        output_id = self._da_id(output_id)
+
+        # Get the index in dinfo
+        i_dac = self._dac_indices[output_id]
+        dinfo, einfo_exists, einfo = self._epoch_functions
+
+        # Get the time and data conversion factors
+        units = myokit.units
+        tf = myokit.Unit.conversion_factor(units.s, tu)
+        if myokit.Unit.can_convert(self._da_units[output_id], units.V):
+            df = myokit.Unit.conversion_factor(self._da_units[output_id], vu)
+        elif myokit.Unit.can_convert(
+                self._da_units[output_id], units.A):  # pragma: no cover
+            df = myokit.Unit.conversion_factor(self._da_units[output_id], cu)
+        else:  # pragma: no cover
+            # Not a voltage or current? Then don't convert
+            df = 1
+        tf, df = float(tf), float(df)
+
+        # Axon has the annoying habit of adding some extra holding at the start
+        # We can include this if we want. See _read_6 for details.
+        if self._is_protocol_file:
+            offset = 0
+        else:
+            offset = self._samples_per_channel // 64
+            offset -= offset % self._n_adc
+            if (offset < self._n_adc):  # pragma: no cover
+                # Don't have a test for this, but this is part of the
+                # established procedure.
+                offset = self._n_adc
+            offset /= self._rate
+
+        # Holding level (converted and rounded)
+        holding = round(df * dinfo(i_dac, 'fDACHoldingLevel'), n_digits)
+
+        # Create protocol
+        p = myokit.Protocol()
+        start = 0
+        next_start = self._sweep_start_to_start
+        for i_sweep in range(self._sweeps_per_run):
+            # Start of sweep: secret event at holding potential
+            if include_initial_holding:
+                e_start = round(tf * start, n_digits)
+                e_length = round(tf * offset, n_digits)
+                p.schedule(holding, e_start, e_length)
+                start += offset
+
+            for e in einfo(i_dac):
+                if e['type'] == EPOCH_STEPPED:
+                    dur = e['init_duration'] / self._rate
+                    inc = e['duration_inc'] / self._rate
+                    duration = dur + i_sweep * inc
+                    level = e['init_level'] + e['level_inc'] * i_sweep
+
+                    e_level = round(df * level, n_digits)
+                    e_start = round(tf * start, n_digits)
+                    e_length = round(tf * duration, n_digits)
+                    p.schedule(e_level, e_start, e_length)
+
+                    start += duration
+                # Note: Only other type can be EPOCH_DISABLED at this point
+
+            # End of sweep: event at holding potential
+            e_start = round(tf * start, n_digits)
+            e_length = round(tf * (next_start - start), n_digits)
+            p.schedule(holding, e_start, e_length)
+            start = next_start
+            next_start += self._sweep_start_to_start
+
+        return p
+
+    def da_units(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._da_units)
+        return self._da_units[index]
+
+    def equal_length_sweeps(self):
+        # Always true for ABF
+        return True
+
+    def filename(self):
+        """ Returns this ABF file's filename. """
+        return self._filename
+
+    def log(self, join_sweeps=False, use_names=False, include_da=True):
+        # Docstring in SweepSource
+
+        # Create log, return if no sweeps or channels
+        log = myokit.DataLog()
+        ns = len(self._sweeps)
+        if ns == 0 or (self._n_adc + self._n_dac) == 0:  # pragma: no cover
+            return log
+
+        # Get channel names
+        if use_names:
+            nc = self._n_adc + (self._n_dac if include_da else 0)
+            names = [c.name() for c in self._sweeps[0][:nc]]
+        else:
+            names = [f'{i}.channel' for i in range(self._n_adc)]
+            if include_da:
+                names += [f'{i}.da' for i in range(self._n_dac)]
+
+        # Gather data and return
+        t = self._sweeps[0][0].times()
+        if not join_sweeps:
+            log['time'] = t
+            for i_sweep, sweep in enumerate(self._sweeps):
+                for channel, name in zip(sweep, names):
+                    log[name, i_sweep] = channel.values()
+        else:
+            log['time'] = np.concatenate(
+                [t + i * self._sweep_start_to_start for i in range(ns)])
+            for i_channel, name in enumerate(names):
+                log[name] = np.concatenate(
+                    [sweep[i_channel].values() for sweep in self._sweeps])
+        log.set_time_key('time')
+        return log
+
+    def matplotlib_figure(self):
+        """ Creates and returns a matplotlib figure with this file's data. """
+        import matplotlib.pyplot as plt
+        f = plt.figure()
+        plt.suptitle(self.filename())
+
+        # Plot AD channels
+        ax = plt.subplot(2, 1, 1)
+        ax.set_title('Measured data')
+        times = None
+        for sweep in self._sweeps:
+            for channel in sweep[:self._n_adc]:
+                if times is None:
+                    times = channel.times()
+                plt.plot(times, channel.values())
+
+        # Plot DA channels
+        n = self._n_dac
+        ax = [plt.subplot(2, n, n + 1 + i) for i in range(n)]
+        for sweep in self._sweeps:
+            for i, channel in enumerate(sweep[self._n_adc:]):
+                ax[i].set_title(channel.name())
+                ax[i].plot(times, channel.values())
+
+        return f
+
+    def meta_str(self, show_header=False):
+        """
+        Returns a multi-line string with meta data about this file.
+
+        The optional argument ``show_header`` can be used to add the full
+        header contents to the output.
+        """
+        out = []
+
+        # File info
+        if self._is_protocol_file:
+            out.append(f'Axon Protocol File: {self._filename}')
+        else:
+            out.append(f'Axon Binary File: {self._filename}')
+        out.append(f'ABF Format version {self._version_str}')
+        out.append(f'Recorded on: {self._datetime}')
+
+        # AProtocol info
+        out.append(
+            f'Acquisition mode: {self._mode}: {acquisition_modes[self._mode]}')
+        if self._number_of_trials:
+            out.append(
+                f'Protocol set for {self._number_of_trials} trials,'
+                f' spaced {self._trial_start_to_start}s apart.')
+            out.append(
+                f'    with {self._runs_per_trial} runs per trial,'
+                f' spaced {self._run_start_to_start}s apart.')
+            out.append(
+                f'     and {self._sweeps_per_run} sweeps per run,'
+                f' spaced {self._sweep_start_to_start}s apart.')
+        else:   # pragma: no cover
+            out.append('Protocol data could not be determined.')
+        out.append(f'Sampling rate: {self._rate} Hz')
+
+        # Channel info
+        if len(self._sweeps) > 0:
+
+            # A/D recordings
+            for i, c in enumerate(self._sweeps[0][:self._n_adc]):
+                out.append(f'A/D Channel {i}: "{c._name}"')
+                if c._type:  # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    out.append(f'  Type: {type_mode_names[c._type]}')
+                out.append(f'  Unit: {c._unit}')
+                if c._lopass:
+                    out.append(f'  Low-pass filter: {c._lopass} Hz')
+                if c._cm:
+                    out.append(f'  Cm (telegraphed): {c._cm} pF')
+                if c._rs:   # pragma: no cover
+                    # Cover pragma: Don't have appropriate test file
+                    out.append(f'  Rs (telegraphed): {c._rs}')
+
+            # Reconstructed D/A outputs
+            for i, c in enumerate(self._sweeps[0][self._n_adc:]):
+                out.append(f'D/A Channel {i}: "{c._name}"')
+                if c._type:  # pragma: no cover
                     # Cover pragma: Don't have appropriate test file
-                    if h['nSignalType'] != 0:
-                        s -= a[i]['fSignalOffset']
+                    out.append('  Type: {type_mode_names[c._type]}')
+                out.append(f'  Unit: {c._unit}')
 
-                # Set final offset
-                self._adc_offsets.append(s)
+        # Add full header info
+        if show_header:
+            if self._strings:
+                dict_to_string(out, 'Strings', {'strings': self._strings})
+            dict_to_string(out, 'file header', self._header)
+
+        return '\n'.join(out)
+
+    def path(self):
+        """ Returns the path to the underlying ABF file. """
+        return self._filepath
 
     def sweep_count(self):
-        # See parent for docstring
+        # Docstring in SweepSource
         return len(self._sweeps)
 
+    def time_unit(self):
+        # Docstring in SweepSource
+        # For ABF, this is always seconds
+        return myokit.units.s
+
+    def _unit(self, unit_string):
+        """ Parses a unit string and returns a :class:`myokit.Unit`. """
+        try:
+            return self._unit_cache[unit_string]
+        except KeyError:
+            unit = myokit.parse_unit(unit_string.replace(MU, 'u'))
+            self._unit_cache[unit_string] = unit
+            return unit
+
     def version(self):
         """ Returns a string representation of this file's version number. """
         return self._version_str
@@ -1334,15 +1353,10 @@ class Sweep:
     """
     Represents a single sweep (also called an *episode*).
 
-    Each sweep contains a number of :class:`channels<Channel>`.
+    Each sweep contains a fixed number of :class:`channels<Channel>`.
     """
-    def __init__(self, n):
-        super().__init__()
-        n = int(n)
-        if n < 0:   # pragma: no cover
-            raise ValueError('Number channels cannot be negative.')
-        self._nc = n    # Number of channels
-        self._channels = [None] * n
+    def __init__(self):
+        self._channels = []
 
     def __getitem__(self, key):
         return self._channels[key]  # Handles slices etc.
@@ -1351,31 +1365,26 @@ def __iter__(self):
         return iter(self._channels)
 
     def __len__(self):
-        return self._nc
-
-    def __setitem__(self, key, value):
-        if type(key) == slice:
-            raise NotImplementedError(
-                'Assignment with slices is not supported.')
-        self._channels[key] = value
+        return len(self._channels)
 
 
 class Channel:
     """
-    Represents an analog signal for a single channel.
+    Represents a signal for a single channel.
 
-    To obtain this channel's data, use :meth:`times` and :meth:`trace`.
+    To obtain its data, use :meth:`times` and :meth:`values`.
     """
     def __init__(self, parent_file):
-        super().__init__()
         self._parent_file = parent_file  # The abf file this channel is from
         self._type = TYPE_UNKNOWN   # Type of recording
 
         # This channel's name
         self._name = None
 
-        # This channel's index (see note below)
-        self._numb = None
+        # This channel's index in the file. This is basically a name, and does
+        # not correspond to e.g. its index in the ADC/DAC info or its index in
+        # the sweep's list of channels.
+        self._index = None
 
         # The units this channel's data is in
         self._unit = None
@@ -1398,36 +1407,32 @@ def __init__(self, parent_file):
         # The reported low-pass filter cut-off frequency
         self._lopass = None
 
-        # Note that the channel indices are not necessarily sequential! So a
-        # file with 2 channels can have indices 0 and 3.
+    def index(self):
+        """ Returns the index set for this channel. """
+        return self._index
 
     def name(self):
         """ Returns the name set for this channel. """
         return self._name
 
-    def number(self):
-        """
-        Returns the channel index used by pClamp.
-
-        Note that this does not necessarily equal its index in the Python sweep
-        data.
-        """
-        return self._numb
-
     def __str__(self):
         return (
-            f'Channel({self._numb} "{self._name}"); {len(self._data)}'
-            f' points sampled at {self._rate}Hz, starts at t={self._start}.')
+            f'Channel({self._index} "{self._name}"); {len(self._data)} points'
+            f' sampled at {self._rate}Hz, starts at t={self._start}.')
 
     def times(self):
         """ Returns a copy of the values on the time axis. """
         n = len(self._data)
-        f = 1.0 / self._rate
+        f = 1 / self._rate
         return np.arange(self._start, self._start + n * f, f)[0:n]
 
+    def unit(self):
+        """ Returns the units this channel is in. """
+        return self._unit
+
     def values(self):
         """ Returns a copy of the values on the data axis. """
-        return np.array(self._data, copy=True)
+        return np.copy(self._data)
 
 
 def dict_to_string(out, name, d, tab=''):
@@ -1436,7 +1441,7 @@ def dict_to_string(out, name, d, tab=''):
     out.append(f'{tab}{"-" * m}  {name}  {"-" * m}')
     for n, v in d.items():
         n = str(n)
-        if type(v) == OrderedDict:
+        if type(v) == dict:
             dict_to_string(out, n, v, f'{tab}  ')
         elif type(v) == list:
             list_to_string(out, n, v, tab)
@@ -1450,7 +1455,7 @@ def list_to_string(out, name, d, tab=''):
     """ Used by AbfFile.info(). """
     for index, item in enumerate(d):
         n = f'{name}[{index}]'
-        if type(item) == OrderedDict:
+        if type(item) == dict:
             dict_to_string(out, n, item, tab)
         elif type(item) == list:    # pragma: no cover
             # Cover pragma: Don't have appropriate test file
@@ -1468,11 +1473,17 @@ def list_to_string(out, name, d, tab=''):
 BLOCKSIZE = 512
 
 
+# A mu, sometimes found in unit strings
+MU = '\u00b5'
+
+
 # Header fields for versions 1 and 2
 # Stored as (key, offset, format) where format corresponds to a struct
 #  unpacking format as documented in:
 #  http://docs.python.org/library/struct.html#format-characters
-headerFields = {
+HEADER_FIELDS = {
+    # Note that a lot of the groups in the version 1 header start with obsolete
+    # fields, followed later by their newer equivalents.
     1: [
         ('fFileSignature', 0, '4s'),       # Coarse file version indication
         # Group 1, File info and sizes
@@ -1569,7 +1580,7 @@ def list_to_string(out, name, d, tab=''):
         ('nDigitalHolding', 1584, 'h'),
         ('nDigitalInterEpisode', 1586, 'h'),
         ('nDigitalValue', 2588, '10h'),
-        ('lDACFilePtr', 2048, '2i'),            # Pointer to protocol?
+        ('lDACFilePtr', 2048, '2i'),
         ('lDACFileNumEpisodes', 2056, '2i'),
         ('fDACCalibrationFactor', 2074, '4f'),
         ('fDACCalibrationOffset', 2090, '4f'),
@@ -1581,13 +1592,12 @@ def list_to_string(out, name, d, tab=''):
         ('fEpochLevelInc', 2428, '20f'),
         ('lEpochInitDuration', 2508, '20i'),
         ('lEpochDurationInc', 2588, '20i'),
-        # Group 10, DAC Output file
+        # Group 10, DAC Output file (Stimulus file)
         ('fDACFileScale', 2708, 'd'),
         ('fDACFileOffset', 2716, 'd'),
         ('lDACFileEpisodeNum', 2724, 'i'),
         ('nDACFileADCNum', 2732, '2h'),
-        # 256 * 2char = utf8? or 2 strings?
-        ('sDACFilePath', 2736, '256s' * 2),
+        ('sDACFilePath', 2736, '256s' * 2),     # Two strings
         # Group 11,
         # Group 12, User list parameters
         ('nULEnable', 3360, '4h'),
@@ -1606,7 +1616,7 @@ def list_to_string(out, name, d, tab=''):
         # Group 22
         # Group 23 Post-processing
         # Group 24 Legacy stuff
-        # Group 6 again?
+        # Group 6 extended
         ('nTelegraphEnable', 4512, '16h'),
         ('fTelegraphAdditGain', 4576, '16f'),
     ],
@@ -1634,7 +1644,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 File sections
-abf2FileSections = [
+ABF2_FILE_SECTIONS = [
     'Protocol',
     'ADC',
     'DAC',
@@ -1657,7 +1667,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 Fields in the tag section
-TagInfoDescription = [
+ABF2_TAG_INFO_DESCRIPTION = [
     ('lTagTime', 'i'),
     ('sComment', '56s'),
     ('nTagType', 'h'),
@@ -1666,7 +1676,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 Fields in the protocol section
-protocolFields = [
+ABF2_PROTOCOL_FIELDS = [
     ('nOperationMode', 'h'),                # 0
     ('fADCSequenceInterval', 'f'),          # 2
     ('bEnableFileCompression', 'b'),        # 6
@@ -1742,7 +1752,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 Fields in the ADC section
-ADCFields = [
+ABF2_ADC_FIELDS = [
     ('nADCNum', 'h'),
     ('nTelegraphEnable', 'h'),
     ('nTelegraphInstrument', 'h'),
@@ -1774,7 +1784,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 Fields in the DAC section
-DACFields = [
+ABF2_DAC_FIELDS = [
     ('nDACNum', 'h'),
     ('nTelegraphDACScaleFactorEnable', 'h'),
     ('fInstrumentHoldingLevel', 'f'),
@@ -1821,7 +1831,7 @@ def list_to_string(out, name, d, tab=''):
 
 
 # ABF2 Fields in the DAC-Epoch section
-EpochInfoPerDACFields = [
+ABF2_EPOCH_INFO_PER_DAC_FIELD = [
     ('nEpochNum', 'h'),
     ('nDACNum', 'h'),
     ('nEpochType', 'h'),
@@ -1833,7 +1843,9 @@ def list_to_string(out, name, d, tab=''):
     ('lEpochPulseWidth', 'i'),
     ('sUnused', '18s'),
 ]
-UserListFields = [
+
+# ABF2 User list fields
+ABF2_USER_LIST_FIELDS = [
     ('nListNum', 'h'),
     ('nULEnable', 'h'),
     ('nULParamToVary', 'h'),
@@ -1912,6 +1924,12 @@ def list_to_string(out, name, d, tab=''):
 }
 
 
+# DAC waveform types
+DAC_DISABLED = 0
+DAC_EPOCHTABLEWAVEFORM = 1  # Epochs
+DAC_DACFILEWAVEFORM = 2     # Stimulus file
+
+
 # User list parameter to vary
 '''
 CONDITNUMPULSES 0
diff --git a/myokit/formats/axon/_atf.py b/myokit/formats/axon/_atf.py
index 5f1a13acd..8acb38773 100644
--- a/myokit/formats/axon/_atf.py
+++ b/myokit/formats/axon/_atf.py
@@ -61,7 +61,11 @@ def filename(self):
 
     def info(self):
         """ Returns this ATF's header/meta data. """
-        return self._meta
+        # Deprecated since 2023-07-12
+        import warnings
+        warnings.warn(
+            'The method `info` is deprecated. Please use `meta_str` instead.')
+        return self.meta_str()
 
     def items(self):
         """ Returns an iterator over all ``(key, value)`` pairs. """
@@ -82,6 +86,15 @@ def log(self):
             log[k] = v
         return log
 
+    def meta_str(self, verbose=True):
+        """
+        Returns this ATF's header data as an unstructured multi-line string.
+
+        Note that the ``verbose`` argument doesn't do anything, but provides
+        compatibility with similar methods in other files.
+        """
+        return self._meta
+
     def myokit_log(self):
         """ Deprecated alias of :meth:`log`. """
         # Deprecated since 2023-06-22
diff --git a/myokit/formats/axon/_importer.py b/myokit/formats/axon/_importer.py
index 630a9c47c..8b2d3977e 100644
--- a/myokit/formats/axon/_importer.py
+++ b/myokit/formats/axon/_importer.py
@@ -16,7 +16,7 @@ class AbfImporter(myokit.formats.Importer):
     def supports_protocol(self):
         return True
 
-    def protocol(self, filename, channel=None):
+    def protocol(self, filename, channel=0):
         """
         Attempts to load the protocol from the file at ``filename``.
 
@@ -25,4 +25,4 @@ def protocol(self, filename, channel=None):
         """
         from myokit.formats.axon import AbfFile
         abf = AbfFile(filename)
-        return abf.protocol(channel)
+        return abf.da_protocol(channel)
diff --git a/myokit/formats/wcp/_wcp.py b/myokit/formats/wcp/_wcp.py
index f0983faa8..b1bc03830 100755
--- a/myokit/formats/wcp/_wcp.py
+++ b/myokit/formats/wcp/_wcp.py
@@ -14,7 +14,7 @@
 _ENC = 'ascii'
 
 
-class WcpFile:
+class WcpFile(myokit.formats.SweepSource):
     """
     Represents a read-only WinWCP file (``.wcp``), stored at ``path``.
 
@@ -37,23 +37,30 @@ def __init__(self, path):
         self._filename = os.path.basename(path)
 
         # Header info, per file, per record, and per channel
+        self._version_str = None
         self._header = {}
         self._header_raw = {}
         self._record_headers = []
         self._channel_headers = []
 
         # Records
-        self._records = None
+        self._records = []
         self._nr = None     # Records in file
         self._nc = None     # Channels per record
         self._np = None     # Samples per channel
         self._dt = None     # Sampling interval (s)
 
+        # Channels
+        self._channel_names = []
+        self._channel_name_map = {}
+
         # Time signal
         self._time = None
 
-        # Version string
-        self._version_str = None
+        # Units
+        self._time_unit = None
+        self._channel_units = []
+        self._unit_cache = {}
 
         # Read the file
         with open(path, 'rb') as f:
@@ -119,9 +126,14 @@ def _read(self, f):
         except KeyError:
             self._np = (self._header['nbd'] * 512) // (2 * self._nc)
 
+        # Get time units
+        self._time_unit = myokit.units.s
+        # Time as set by dt (which is what we need) is _always_ in seconds.
+        # John Dempster says: "The TU= value referred to the time units which
+        # were displayed in early versions of WinWCP and no longer exists in
+        # recent WCP data files." (Email to michael, 2023-07-12).
+
         # Get channel-specific fields
-        self._channel_names = []
-        self._channel_name_map = {}
         for i in range(self._nc):
             c = {}
             for k, t in HEADER_CHANNEL_FIELDS.items():
@@ -129,6 +141,7 @@ def _read(self, f):
             self._channel_headers.append(c)
             self._channel_names.append(c['yn'])
             self._channel_name_map[c['yn']] = i
+            self._channel_units.append(self._unit(c['yu']))
 
         # Analysis block size and data block size
         # Data is stored as 16 bit integers (little-endian)
@@ -145,7 +158,6 @@ def _read(self, f):
         adcmax = self._header['adcmax']
 
         # Read data records
-        records = []
         offset = header_size
         for i in range(self._nr):
             # Read analysis block
@@ -163,9 +175,15 @@ def _read(self, f):
             # Time of recording, as float, not sure how to interpret
             rtime = struct.unpack(str('<f'), f.read(4))[0]
 
-            # Sampling interval: pretty sure this should be the same as the
-            # file wide one in header['dt']
+            # Sampling interval
+            # It is technically possible to have different dts for different
+            # records (see email J.D. 2023-07-12), but rare.
+            # Not supported here!
             rint = round(struct.unpack(str('<f'), f.read(4))[0], 6)
+            if rint != self._dt:  # pragma: no cover
+                raise ValueError(
+                    'Unsupported feature: WCP file contains more than one'
+                    ' sampling rate.')
 
             # Maximum positive limit of A/D converter voltage range
             vmax = struct.unpack(
@@ -179,7 +197,6 @@ def _read(self, f):
                 'status': rstatus,
                 'type': rtype,
                 'rtime': rtime,
-                'dt': rint,
                 'marker': marker,
             })
 
@@ -201,13 +218,11 @@ def _read(self, f):
                 vmx / (adcmax * h['yg']) * data[:, h['yo']].astype('f4')
                 for h, vmx in zip(self._channel_headers, vmax)]
 
-            records.append(record)
+            self._records.append(record)
 
             # Increase offset beyong data block
             offset += rdb_size
 
-        self._records = records
-
         # Create time signal
         self._time = np.arange(self._np) * self._dt
 
@@ -220,25 +235,30 @@ def __iter__(self):
     def __len__(self):
         return self._nr
 
+    def _channel_id(self, channel_id):
+        """ Checks an int or str channel id and returns a valid int. """
+        if self._nr == 0:  # pragma: no cover
+            raise KeyError(f'Channel {channel_id} not found (empty file).')
+
+        # Handle string
+        if isinstance(channel_id, str):
+            return self._channel_name_map[channel_id]  # Pass KeyError to user
+
+        int_id = int(channel_id)  # TypeError for user
+        if int_id < 0 or int_id >= self._nc:
+            raise IndexError(f'channel_id out of range: {channel_id}')
+        return int_id
+
     def channel(self, channel_id, join_sweeps=False):
         # Docstring in SweepSource
         channel_id = self._channel_id(channel_id)
-
+        time, data = [], []
+        for r, h in zip(self._records, self._record_headers):
+            time.append(self._time + h['rtime'])
+            data.append(r[channel_id])
         if join_sweeps:
-            # Join sweeps
-            time, data = [], []
-            for r, h in zip(self._records, self._record_headers):
-                time.append(self._time + h['rtime'])
-                data.append(r[channel_id])
             return (np.concatenate(time), np.concatenate(data))
-
-        else:
-            # Return individual sweeps
-            out = []
-            out.append(np.array(self._time))
-            for r in self._records:
-                out.append(np.array(r[channel_id]))
-            return tuple(out)
+        return time, data
 
     def channels(self):
         """ Deprecated alias of :meth:`channel_count`. """
@@ -247,80 +267,49 @@ def channels(self):
         warnings.warn(
             'The method `channels` is deprecated. Please use'
             ' WcpFile.channel_count() instead.')
-
         return self._nc
 
     def channel_count(self):
         # Docstring in SweepSource
         return self._nc
 
-    def _channel_id(self, channel_id):
-        """ Checks an int or str channel id and returns a valid int. """
-        if self._nr == 0:  # pragma: no cover
-            raise KeyError(f'Channel {channel_id} not found (empty file).')
+    def channel_names(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._channel_names)
+        return self._channel_names[index]
 
-        # Handle string
-        if isinstance(channel_id, str):
-            return self._channel_name_map[channel_id]  # Pass KeyError to user
+    def channel_units(self, index=None):
+        # Docstring in SweepSource
+        if index is None:
+            return list(self._channel_units)
+        return self._channel_units[index]
 
-        int_id = int(channel_id)  # TypeError for user
-        if int_id < 0 or int_id >= self._nc:
-            raise IndexError(f'channel_id out of range: {channel_id}')
-        return int_id
+    def da_count(self):
+        # Docstring in SweepSource. Rest is allowed raise NotImplementedError
+        return 0
 
-    def channel_names(self):
+    def equal_length_sweeps(self):
         # Docstring in SweepSource
-        return list(self._channel_names)
+        return True
 
     def filename(self):
         """ Returns this file's name. """
         return self._filename
 
     def info(self):
-        """ Returns a multi-line string with meta data about this file. """
-        h = self._header
-        out = []
-
-        # Add file info
-        out.append(f'WinWCP file: {self._filename}')
-        out.append(f'WinWCP Format version {self._version_str}')
-        out.append(f'Recorded on: {h["rtime"]}')
-
-        # Basic records info
-        out.append(f'  Number of records: {self._nr}')
-        out.append(f'  Channels per record: {self._nc}')
-        out.append(f'  Samples per channel: {self._np}')
-        out.append(f'  Sampling interval: {self._dt} s')
-        if 'tu' in h:  # pragma: no cover
-            out.append(f'Time units: {h["tu"]}')
-
-        # Channel info
-        for c in self._channel_headers:
-            out.append(f'A/D channel: {c["yn"]}')
-            out.append(f'  Unit: {c["yu"]}')
-
-        # Record info
-        out.append('Records: Type, Status, Sampling Interval, Start, Marker')
-        for i, r in enumerate(self._record_headers):
-            out.append(f'Record {i}: {r["type"]}, {r["status"]}, {r["dt"]},'
-                       f' {r["rtime"]}, "{r["marker"]}"')
-
-        return '\n'.join(out)
+        """ Deprecated alias of :meth:`meta_str` """
+        # Deprecated since 2023-07-12
+        import warnings
+        warnings.warn(
+            'The method `info` is deprecated. Please use `meta_str` instead.')
+        return self.meta_str(False)
 
-    def log(self, join_sweeps=False, use_names=False, channels=None):
+    def log(self, join_sweeps=False, use_names=False, include_da=True):
         # Docstring in SweepSource
 
-        # Select channels
-        if channels is None:
-            channels = list(range(self._nc))
-        else:
-            channels = [self._channel_id(i) for i in channels]
-
         # Create log
         log = myokit.DataLog()
-
-        # No data? Then return. Do this after channel check to tell user about
-        # invalid channel ids.
         if self._nr == 0:  # pragma: no cover
             return log
 
@@ -328,24 +317,19 @@ def log(self, join_sweeps=False, use_names=False, channels=None):
         channel_names = self._channel_names
         if not use_names:
             channel_names = [f'{i}.channel' for i in range(self._nc)]
-        channel_names = [channel_names[i] for i in channels]
 
-        # Populate log
+        # Gather data and return
         if join_sweeps:
-            # Join sweeps
-            time = []
-            for h in self._record_headers:
-                time.append(self._time + h['rtime'])
-            log['time'] = np.concatenate(time)
-            for c, name in zip(channels, channel_names):
+            log['time'] = np.concatenate(
+                [self._time + h['rtime'] for h in self._record_headers])
+            for c, name in enumerate(channel_names):
                 log[name] = np.concatenate([r[c] for r in self._records])
-
         else:
             # Return individual sweeps
-            log['time'] = np.array(self._time)
-            for i, record in enumerate(self._records):
-                for channel, name in zip(channels, channel_names):
-                    log[name, i] = record[channel]
+            log['time'] = self._time
+            for r, record in enumerate(self._records):
+                for c, name in enumerate(channel_names):
+                    log[name, r] = record[c]
 
         log.set_time_key('time')
         return log
@@ -360,6 +344,46 @@ def matplotlib_figure(self):
                 ax.plot(self._time, channel)
         return f
 
+    def meta_str(self, verbose=False):
+        # Docstring in SweepSource
+        h = self._header
+        out = []
+
+        # Add file info
+        out.append(f'WinWCP file: {self._filename}')
+        out.append(f'WinWCP Format version {self._version_str}')
+        out.append(f'Recorded on: {h["rtime"]}')
+
+        # Basic records info
+        out.append(f'  Number of records: {self._nr}')
+        out.append(f'  Channels per record: {self._nc}')
+        out.append(f'  Samples per channel: {self._np}')
+        out.append(f'  Sampling interval: {self._dt} s')
+
+        # Channel info
+        for c in self._channel_headers:
+            out.append(f'A/D channel: {c["yn"]}')
+            out.append(f'  Unit: {c["yu"]}')
+
+        # Record info
+        out.append(
+            'Records: Type, Status, Sampling Interval, Start, Marker')
+        for i, r in enumerate(self._record_headers):
+            out.append(f'Record {i}: {r["type"]}, {r["status"]},'
+                       f' {r["rtime"]}, "{r["marker"]}"')
+
+        # Parsed and unparsed header
+        if verbose:
+            out.append(f'{"-" * 35} header {"-" * 35}')
+            for k, v in self._header.items():
+                out.append(f'{k}: {v}')
+
+            out.append(f'{"-" * 33} raw header {"-" * 33}')
+            for k, v in self._header_raw.items():
+                out.append(f'{k}: {v}')
+
+        return '\n'.join(out)
+
     def myokit_log(self):
         """
         Deprecated method. Please use ``WcpFile.log(use_names=True)`` instead.
@@ -369,6 +393,7 @@ def myokit_log(self):
         warnings.warn(
             'The method `myokit_log` is deprecated. Please use'
             ' WcpFile.log(use_names=True) instead.')
+        return self.log(self)
 
     def path(self):
         """ Returns the path to this file. """
@@ -416,10 +441,23 @@ def sweep_count(self):
         # Docstring in SweepSource
         return self._nr
 
+    def time_unit(self):
+        # Docstring in SweepSource
+        return self._time_unit
+
     def times(self):
         """ Returns the time points sampled at. """
         return np.array(self._time)
 
+    def _unit(self, unit_string):
+        """ Parses a unit string and returns a :class:`myokit.Unit`. """
+        try:
+            return self._unit_cache[unit_string]
+        except KeyError:
+            unit = myokit.parse_unit(unit_string)
+            self._unit_cache[unit_string] = unit
+            return unit
+
     def values(self, record, channel):
         """
         Returns the values of channel ``channel``, recorded in record
@@ -446,7 +484,6 @@ def version(self):
     'np': int,         # No. of A/D samples per channel
     'dt': float,       # A/D sampling interval (s)
     'nz': int,         # No. of samples averaged to calculate a zero level.
-    'tu': str,         # Time units
     'id': str,         # Experiment identification line
 }
 
diff --git a/myokit/gui/datalog_viewer.py b/myokit/gui/datalog_viewer.py
index 30db8492d..0c8edc782 100755
--- a/myokit/gui/datalog_viewer.py
+++ b/myokit/gui/datalog_viewer.py
@@ -546,9 +546,8 @@ def fileTabChangeEvent(self, index):
 
 
 class TabWidget(QtWidgets.QTabWidget):
-    """
-    Tab widget that can move up and down when asked.
-    """
+    """ Generic tab widget with first/last next/previous methods. """
+
     def first(self):
         """ Select the first widget. """
         self.setCurrentIndex(0)
@@ -574,86 +573,144 @@ def previous(self):
         self.setCurrentIndex(n - 1 if i < 0 else i)
 
 
-class AbfTab(TabWidget):
-    """
-    A widget displaying an ABF file.
-    """
-    def __init__(self, parent, abf):
+class GraphTabWidget(TabWidget):
+    """ Tab widget to graph a data source. """
+
+    def __init__(self, parent):
         super().__init__(parent)
+
         self.setTabsClosable(False)
         self.setTabPosition(self.TabPosition.East)
-        self._abf = abf
+
         self._figures = []
         self._axes = []
-        for i in range(self._abf.channel_count()):
-            tab, name = self.create_graph_tab(i)
-            self.addTab(tab, name)
-        self.addTab(self.create_info_tab(), 'Info')
-        del self._abf
 
-    def create_graph_tab(self, channel):
-        """
-        Creates a widget displaying the main data.
-        """
+    def deleteLater(self):
+        """ Deletes this tab (later). """
+        for figure in self._figures:
+            figure.clear()
+        for axes in self._axes:
+            axes.cla()
+        del self._figures, self._axes
+        gc.collect()
+        super().deleteLater()
+
+
+class SweepSourceTab(GraphTabWidget):
+    """ A tab widget for sources implementing the SweepSource interface. """
+
+    def __init__(self, parent, source):
+        super().__init__(parent)
+
+        # Add A/D
+        for i in range(source.channel_count()):
+            self._add_graph_tab(source, i)
+
+        # Add D/A
+        for i in range(source.da_count()):
+            self._add_graph_tab(source, i, True)
+
+        # Add meta data
+        self._add_meta_tab(source)
+
+    def _add_graph_tab(self, source, index, da=False):
+        """ Adds a tab for a graph. """
+
+        # Create widget
         widget = QtWidgets.QWidget(self)
 
         # Create figure
         figure = matplotlib.figure.Figure()
-        figure.suptitle(self._abf.filename())
         canvas = backend.FigureCanvasQTAgg(figure)
         canvas.setParent(widget)
+        toolbar = backend.NavigationToolbar2QT(canvas, widget)
+
+        # Draw signal
+        join_sweeps = not source.equal_length_sweeps()
+        if da:
+            name = source.da_names(index)
+            units = source.da_units(index)
+            times, values = source.da(index, join_sweeps)
+        else:
+            name = source.channel_names(index)
+            units = source.channel_units(index)
+            times, values = source.channel(index, join_sweeps)
+
         axes = figure.add_subplot(1, 1, 1)
+        axes.set_xlabel(f'Time {source.time_unit()}')
+        axes.set_ylabel(f'{name} {units}')
+        if join_sweeps:
+            axes.plot(times, values)
+        else:
+            for v in values:
+                axes.plot(times[0], v)
+
+        # Store for later deletion
+        self._figures.append(figure)
+        self._axes.append(axes)
+
+        # Create a layout
+        vbox = QtWidgets.QVBoxLayout()
+        vbox.addWidget(canvas)
+        vbox.addWidget(toolbar)
+        widget.setLayout(vbox)
+
+        # Add tab
+        self.addTab(widget, name)
+
+    '''
+    def debug_tab(self, channel_index, da=True):
+        """ Add a tab graphing data using the DataLog method. """
+        widget = QtWidgets.QWidget(self)
+
+        # Create widget
+        widget = QtWidgets.QWidget(self)
+
+        # Create figure
+        figure = matplotlib.figure.Figure()
+        figure.suptitle(self._abf.filename())
+        canvas = backend.FigureCanvasQTAgg(figure)
+        canvas.setParent(widget)
         toolbar = backend.NavigationToolbar2QT(canvas, widget)
 
         # Draw lines
-        name = f'Channel {channel}'  # Default if no data is present
-        times = None
-        for i, sweep in enumerate(self._abf):
-            if times is None:
-                name = sweep[channel].name()
-                times = sweep[channel].times()
-            axes.plot(times, sweep[channel].values())
+        p = self._abf.da_protocol(channel_index, tu='s', vu='V')
+        times, _ = self._abf.da(channel_index)
+        for t in times:
+            d = p.log_for_interval(t[0], t[-1], for_drawing=True).npview()
+            axes.plot(d['time'] - t[0], d['pace'])
 
         # Create a layout
         vbox = QtWidgets.QVBoxLayout()
         vbox.addWidget(canvas)
         vbox.addWidget(toolbar)
         widget.setLayout(vbox)
-        self._figures.append(figure)
-        self._axes.append(axes)
-        return widget, name
 
-    def create_info_tab(self):
-        """ Creates a tab displaying information about the file. """
-        widget = QtWidgets.QTextEdit(self)
-        widget.setText(self._abf.info(show_header=True))
-        widget.setReadOnly(True)
-        return widget
+        # Add tab
+        self.addTab(widget, name)
+'''
 
-    def deleteLater(self):
-        """ Deletes this tab (later). """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
+    def _add_meta_tab(self, source):
 
+        meta = source.meta_str(True)
+        if meta:
+            widget = QtWidgets.QTextEdit(self)
+            widget.setText(meta)
+            widget.setReadOnly(True)
+            self.addTab(widget, 'info')
 
-class AtfTab(TabWidget):
-    """
-    A widget displaying an AGF file.
-    """
+
+class AbfTab(SweepSourceTab):
+    """ A widget displaying an ABF file. """
+    pass
+
+
+class AtfTab(GraphTabWidget):
+    """ A widget displaying an ATF file. """
     def __init__(self, parent, atf):
         super().__init__(parent)
         self._atf = atf
 
-        self.setTabsClosable(False)
-        self.setTabPosition(self.TabPosition.East)
-
-        self._figures = []
-        self._axes = []
         keys = list(self._atf.keys())
         if len(keys) > 1:
             time = keys[0]  # Time is always first (and regularly sampled)
@@ -697,18 +754,8 @@ def create_info_tab(self):
         widget.setReadOnly(True)
         return widget
 
-    def deleteLater(self):
-        """ Deletes this tab (later). """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
-
 
-class CsvTab(TabWidget):
+class CsvTab(GraphTabWidget):
     """
     A widget displaying a CSV file.
 
@@ -716,12 +763,9 @@ class CsvTab(TabWidget):
     """
     def __init__(self, parent, log, filename):
         super().__init__(parent)
-        self.setTabsClosable(False)
-        self.setTabPosition(self.TabPosition.East)
+
         self._log = log.npview()
         self._filename = filename
-        self._figures = []
-        self._axes = []
 
         # Check time key was found
         time = log.time_key()
@@ -792,30 +836,13 @@ def create_graph_tab(self, key, indices=None):
         self._axes.append(axes)
         return widget
 
-    def deleteLater(self):
-        """
-        Deletes this tab (later).
-        """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
-
 
-class MatTab(TabWidget):
-    """
-    A widget displaying a MAT file.
-    """
+class MatTab(GraphTabWidget):
+    """ A widget displaying a .mat file. """
     def __init__(self, parent, mat, filename):
         super().__init__(parent)
-        self.setTabsClosable(False)
-        self.setTabPosition(self.TabPosition.East)
-        self._figures = []
+
         self._filename = filename
-        self._axes = []
 
         # Find usable data
         for key in mat.keys():
@@ -880,30 +907,14 @@ def create_graph_tab(self, time, data):
         self._axes.append(axes)
         return widget
 
-    def deleteLater(self):
-        """
-        Deletes this tab (later).
-        """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
 
+class TxtTab(GraphTabWidget):
+    """ A widget displaying a .txt file (with lots of heuristics!). """
 
-class TxtTab(TabWidget):
-    """
-    A widget displaying a TXT file (with lots of heuristics!).
-    """
     def __init__(self, parent, data, filename):
         super().__init__(parent)
-        self.setTabsClosable(False)
-        self.setTabPosition(self.TabPosition.East)
-        self._figures = []
+
         self._filename = filename
-        self._axes = []
 
         # Find usable data
         if np.prod(data.shape) == np.max(data.shape):
@@ -962,76 +973,7 @@ def create_graph_tab(self, time, data):
         self._axes.append(axes)
         return widget
 
-    def deleteLater(self):
-        """
-        Deletes this tab (later).
-        """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
-
-
-class WcpTab(TabWidget):
-    """
-    A widget displaying a WCP file.
-    """
-    def __init__(self, parent, wcp):
-        super().__init__(parent)
-        self.setTabsClosable(False)
-        self.setTabPosition(self.TabPosition.East)
-        self._wcp = wcp
-        self._figures = []
-        self._axes = []
-        for i in range(self._wcp.record_count()):
-            self.addTab(self.create_graph_tab(i), 'Record ' + str(i))
-        del self._wcp
-
-    def create_graph_tab(self, record):
-        """ Creates a widget displaying the data in the given ``record``. """
-        widget = QtWidgets.QWidget(self)
-
-        # Create figure
-        figure = matplotlib.figure.Figure()
-        figure.suptitle(self._wcp.filename())
-        canvas = backend.FigureCanvasQTAgg(figure)
-        canvas.setParent(widget)
-        axes = figure.add_subplot(1, 1, 1)
-        toolbar = backend.NavigationToolbar2QT(canvas, widget)
-
-        # Draw lines
-        for i in range(self._wcp.channel_count()):
-            axes.plot(
-                np.array(self._wcp.times(), copy=True),
-                np.array(self._wcp.values(record, i), copy=True),
-            )
-
-        # Create a layout
-        vbox = QtWidgets.QVBoxLayout()
-        vbox.addWidget(canvas)
-        vbox.addWidget(toolbar)
-        widget.setLayout(vbox)
-        self._figures.append(figure)
-        self._axes.append(axes)
-        return widget
-
-    def create_info_tab(self):
-        """ Creates a tab displaying information about the file. """
-        widget = QtWidgets.QTextEdit(self)
-        widget.setText(self._wcp.info())
-        widget.setReadOnly(True)
-        return widget
 
-    def deleteLater(self):
-        """ Deletes this tab (later). """
-        for figure in self._figures:
-            figure.clear()
-        for axes in self._axes:
-            axes.cla()
-        del self._figures, self._axes
-        gc.collect()
-        super().deleteLater()
+class WcpTab(SweepSourceTab):
+    pass
 


Test patch:
diff --git a/myokit/tests/data/formats/abf-v2.abf b/myokit/tests/data/formats/abf-v2.abf
index a98880e54..5d4502343 100644
Binary files a/myokit/tests/data/formats/abf-v2.abf and b/myokit/tests/data/formats/abf-v2.abf differ
diff --git a/myokit/tests/test_formats_axon.py b/myokit/tests/test_formats_axon.py
index 386171553..ad35785a9 100755
--- a/myokit/tests/test_formats_axon.py
+++ b/myokit/tests/test_formats_axon.py
@@ -26,9 +26,9 @@
      and 9 sweeps per run, spaced 0.5s apart.
 Sampling rate: 10000.0 Hz
 A/D Channel 0: "IN 0"
-  Unit: pA
+  Unit: [pA]
 D/A Channel 0: "OUT 0"
-  Unit: mV
+  Unit: [mV]
 '''.strip()
 
 V1_PROTO_INFO = '''
@@ -41,7 +41,7 @@
      and 30 sweeps per run, spaced 5.0s apart.
 Sampling rate: 20000.0 Hz
 D/A Channel 0: "Cmd 0"
-  Unit: mV
+  Unit: [mV]
 '''.strip()
 
 V1_PROTOCOL = '''
@@ -67,28 +67,124 @@
 0.0      4100.0   400.0    0.0      0
 '''.strip()
 
+V1_PROTOCOL_VS = '''
+[[protocol]]
+# Level  Start    Length   Period   Multiplier
+-0.1     0.0      0.1      0.0      0
+0.0      0.1      0.4      0.0      0
+-0.08    0.5      0.1      0.0      0
+0.0      0.6      0.4      0.0      0
+-0.06    1.0      0.1      0.0      0
+0.0      1.1      0.4      0.0      0
+-0.04    1.5      0.1      0.0      0
+0.0      1.6      0.4      0.0      0
+-0.02    2.0      0.1      0.0      0
+0.0      2.1      0.4      0.0      0
+0.0      2.5      0.1      0.0      0
+0.0      2.6      0.4      0.0      0
+0.02     3.0      0.1      0.0      0
+0.0      3.1      0.4      0.0      0
+0.04     3.5      0.1      0.0      0
+0.0      3.6      0.4      0.0      0
+0.06     4.0      0.1      0.0      0
+0.0      4.1      0.4      0.0      0
+'''.strip()
+
 V2_INFO = '''
 Axon Binary File: abf-v2.abf
 ABF Format version 2.0.0.0
-Recorded on: 2014-10-01 14:03:55.980999
+Recorded on: 2016-01-07 10:51:55.345000
 Acquisition mode: 5: Episodic stimulation mode
 Protocol set for 1 trials, spaced 0.0s apart.
     with 1 runs per trial, spaced 0.0s apart.
-     and 1 sweeps per run, spaced 5.0s apart.
-Sampling rate: 10000.0 Hz
+     and 37 sweeps per run, spaced 5.0s apart.
+Sampling rate: 20000.0 Hz
 A/D Channel 0: "IN 0"
-  Unit: pA
+  Unit: [pA]
   Low-pass filter: 10000.0 Hz
-  Cm (telegraphed): 6.34765625 pF
+  Cm (telegraphed): 11.083984375 pF
 D/A Channel 0: "Cmd 0"
-  Unit: mV
+  Unit: [mV]
+
 '''.strip()
 
 V2_PROTOCOL = '''
 [[protocol]]
 # Level  Start    Length   Period   Multiplier
-12.0     0.0      1000.0   0.0      0
--120.0   1000.0   4000.0   0.0      0
+-100.0   0.0      25.0     0.0      0
+-120.0   25.0     4975.0   0.0      0
+-95.0    5000.0   25.0     0.0      0
+-120.0   5025.0   4975.0   0.0      0
+-90.0    10000.0  25.0     0.0      0
+-120.0   10025.0  4975.0   0.0      0
+-85.0    15000.0  25.0     0.0      0
+-120.0   15025.0  4975.0   0.0      0
+-80.0    20000.0  25.0     0.0      0
+-120.0   20025.0  4975.0   0.0      0
+-75.0    25000.0  25.0     0.0      0
+-120.0   25025.0  4975.0   0.0      0
+-70.0    30000.0  25.0     0.0      0
+-120.0   30025.0  4975.0   0.0      0
+-65.0    35000.0  25.0     0.0      0
+-120.0   35025.0  4975.0   0.0      0
+-60.0    40000.0  25.0     0.0      0
+-120.0   40025.0  4975.0   0.0      0
+-55.0    45000.0  25.0     0.0      0
+-120.0   45025.0  4975.0   0.0      0
+-50.0    50000.0  25.0     0.0      0
+-120.0   50025.0  4975.0   0.0      0
+-45.0    55000.0  25.0     0.0      0
+-120.0   55025.0  4975.0   0.0      0
+-40.0    60000.0  25.0     0.0      0
+-120.0   60025.0  4975.0   0.0      0
+-35.0    65000.0  25.0     0.0      0
+-120.0   65025.0  4975.0   0.0      0
+-30.0    70000.0  25.0     0.0      0
+-120.0   70025.0  4975.0   0.0      0
+-25.0    75000.0  25.0     0.0      0
+-120.0   75025.0  4975.0   0.0      0
+-20.0    80000.0  25.0     0.0      0
+-120.0   80025.0  4975.0   0.0      0
+-15.0    85000.0  25.0     0.0      0
+-120.0   85025.0  4975.0   0.0      0
+-10.0    90000.0  25.0     0.0      0
+-120.0   90025.0  4975.0   0.0      0
+-5.0     95000.0  25.0     0.0      0
+-120.0   95025.0  4975.0   0.0      0
+0.0      100000.0 25.0     0.0      0
+-120.0   100025.0 4975.0   0.0      0
+5.0      105000.0 25.0     0.0      0
+-120.0   105025.0 4975.0   0.0      0
+10.0     110000.0 25.0     0.0      0
+-120.0   110025.0 4975.0   0.0      0
+15.0     115000.0 25.0     0.0      0
+-120.0   115025.0 4975.0   0.0      0
+20.0     120000.0 25.0     0.0      0
+-120.0   120025.0 4975.0   0.0      0
+25.0     125000.0 25.0     0.0      0
+-120.0   125025.0 4975.0   0.0      0
+30.0     130000.0 25.0     0.0      0
+-120.0   130025.0 4975.0   0.0      0
+35.0     135000.0 25.0     0.0      0
+-120.0   135025.0 4975.0   0.0      0
+40.0     140000.0 25.0     0.0      0
+-120.0   140025.0 4975.0   0.0      0
+45.0     145000.0 25.0     0.0      0
+-120.0   145025.0 4975.0   0.0      0
+50.0     150000.0 25.0     0.0      0
+-120.0   150025.0 4975.0   0.0      0
+55.0     155000.0 25.0     0.0      0
+-120.0   155025.0 4975.0   0.0      0
+60.0     160000.0 25.0     0.0      0
+-120.0   160025.0 4975.0   0.0      0
+65.0     165000.0 25.0     0.0      0
+-120.0   165025.0 4975.0   0.0      0
+70.0     170000.0 25.0     0.0      0
+-120.0   170025.0 4975.0   0.0      0
+75.0     175000.0 25.0     0.0      0
+-120.0   175025.0 4975.0   0.0      0
+80.0     180000.0 25.0     0.0      0
+-120.0   180025.0 4975.0   0.0      0
 '''.strip()
 
 
@@ -108,78 +204,97 @@ def test_read_v1(self):
         self.assertIsInstance(abf, myokit.formats.SweepSource)
 
         # Check getting info
-        self.assertIn('version 1.65', abf.info())
+        self.assertIn('version 1.65', abf.meta_str())
         self.maxDiff = None
-        self.assertEqual(abf.info(), V1_INFO)
+        self.assertEqual(abf.meta_str(), V1_INFO)
 
         # Test getting full header runs without crashing
-        abf.info(True)
+        abf.meta_str(True)
 
         # Get version
         self.assertEqual(abf.version(), '1.65')
 
-        # Test len returns number of sweeps
+        # Sweep count
         self.assertEqual(len(abf), 9)
+        self.assertEqual(abf.sweep_count(), 9)
+        self.assertEqual(len([s for s in abf]), 9)
 
-        # Test access to A/D channels via native API
-        self.assertEqual(abf.ad_channel_count(), 1)    # 1 data channel
+        # Test access to A/D and D/A channels via sequence interface
         self.assertIsInstance(abf[0], axon.Sweep)
+        self.assertIsInstance(abf[8], axon.Sweep)
+        self.assertFalse(abf[0] is abf[1])
+        self.assertEqual(len(abf[0]), 2)
+        self.assertEqual(len(abf[8]), 2)
         self.assertIsInstance(abf[0][0], axon.Channel)
-        self.assertEqual(abf.da_channel_count(), 1)    # 1 D/A channel
         self.assertIsInstance(abf[0][1], axon.Channel)
-        self.assertEqual(len([s for s in abf]), 9)
-        self.assertEqual(len(abf[0]), abf.channel_count())
-
-        # Test conversion to Myokit protocol
-        p = abf.protocol(1)
-        self.assertEqual(len(p), 18)
-        self.assertEqual(p.code(), V1_PROTOCOL)
 
-        # Test other D/A channel methods
-        self.assertEqual(abf.da_holding_level(1), 0)
-        self.assertEqual(len(abf.da_steps(1)[0]), 9)
-        self.assertEqual(
-            list(abf.da_steps(1)[0]),
-            [-100, -80, -60, -40, -20, 0, 20, 40, 60])
-
-        # Test Channel methods
+        # Test abf.Channel methods
         channel = abf[0][0]
-        self.assertIsInstance(channel.number(), int)
+        self.assertEqual(channel.index(), 0)
         self.assertEqual(channel.name(), 'IN 0')
+        self.assertIsInstance(channel.unit(), myokit.Unit)
+        self.assertEqual(channel.unit(), myokit.units.pA)
         self.assertEqual(str(channel),
                          'Channel(0 "IN 0"); 5000 points sampled at 10000.0Hz,'
                          ' starts at t=0.0.')
         self.assertEqual(len(channel.times()), len(channel.values()))
         self.assertFalse(np.all(channel.times() == channel.values()))
 
-        # Test SweepSource interface
-        self.assertEqual(abf.sweep_count(), 9)
-        self.assertEqual(abf.channel_count(), 2)
-        self.assertEqual(abf.channel_names(), ['IN 0', 'OUT 0'])
-
-        # Channel without joining
-        x = abf.channel(0)
-        self.assertEqual(len(x), 1 + len(abf))
-        self.assertEqual(len(x[0]), len(x[1]))
-        self.assertEqual(len(x[0]), len(x[2]))
-        self.assertEqual(len(x[0]), len(x[3]))
-        self.assertEqual(len(x[0]), len(x[4]))
-        self.assertEqual(len(x[0]), len(x[5]))
-        self.assertEqual(len(x[0]), len(x[6]))
-        self.assertEqual(len(x[0]), len(x[7]))
-        self.assertEqual(len(x[0]), len(x[8]))
-        self.assertEqual(len(x[0]), len(x[9]))
-
-        # Channel with joining
-        y = abf.channel(0, join_sweeps=True)
-        self.assertEqual(len(y), 2)
-        self.assertEqual(len(y[0]), 9 * len(x[0]))
-        self.assertEqual(len(y[0]), len(y[1]))
+        # Test SweepSource info
+        self.assertTrue(abf.equal_length_sweeps())
+        self.assertEqual(abf.time_unit(), myokit.units.s)
+
+        # Test A/D channel access via SweepSource interface
+        self.assertEqual(abf.channel_count(), 1)
+        self.assertEqual(abf.channel_names(), ['IN 0'])
+        self.assertEqual(abf.channel_names(0), 'IN 0')
+        self.assertEqual(abf.channel_units(), [myokit.units.pA])
+        self.assertEqual(abf.channel_units(0), myokit.units.pA)
+
+        times, values = abf.channel(0)
+        self.assertEqual(9, len(times), len(values))
+        self.assertEqual(len(times[0]), len(values[0]))
+        self.assertEqual(len(times[0]), len(times[1]), len(values[1]))
+        self.assertEqual(len(times[0]), len(times[2]), len(values[2]))
+        self.assertEqual(len(times[0]), len(times[3]), len(values[3]))
+        self.assertEqual(len(times[0]), len(times[4]), len(values[4]))
+        self.assertEqual(len(times[0]), len(times[5]), len(values[5]))
+        self.assertEqual(len(times[0]), len(times[6]), len(values[6]))
+        self.assertEqual(len(times[0]), len(times[7]), len(values[7]))
+        self.assertEqual(len(times[0]), len(times[8]), len(values[8]))
+
+        tj, vj = abf.channel(0, join_sweeps=True)
+        self.assertTrue(np.all(tj == np.concatenate(times)))
+        self.assertTrue(np.all(vj == np.concatenate(values)))
 
         # Channel doesn't exist
         self.assertRaises(IndexError, abf.channel, -1)
-        self.assertRaises(IndexError, abf.channel, 2)
-        self.assertRaises(KeyError, abf.channel, 'Tom')
+        self.assertRaises(IndexError, abf.channel, 1)
+        self.assertRaises(KeyError, abf.channel, 'OUT 0')
+
+        # Test D/A output access via SweepSource interface
+        self.assertEqual(abf.da_count(), 1)
+        self.assertEqual(abf.da_names(), ['OUT 0'])
+        self.assertEqual(abf.da_names(0), 'OUT 0')
+        self.assertEqual(abf.da_units(), [myokit.units.mV])
+        self.assertEqual(abf.da_units(0), myokit.units.mV)
+
+        da_times, da_values = abf.da(0)
+        self.assertEqual(9, len(da_times), len(da_values))
+        self.assertEqual(len(da_times[0]), len(da_values[0]))
+        self.assertEqual(len(da_times[0]), len(da_times[1]), len(da_values[1]))
+        self.assertEqual(len(da_times[0]), len(da_times[8]), len(da_values[8]))
+        self.assertEqual(len(da_times[0]), len(times[0]))
+        self.assertTrue(np.all(times[0] == da_times[0]))
+        self.assertFalse(np.all(values[0] == da_values[0]))
+        self.assertTrue(np.all(da_values[3] == abf.da('OUT 0')[1][3]))
+        self.assertRaises(IndexError, abf.da, -1)
+        self.assertRaises(IndexError, abf.da, 1)
+        self.assertRaises(KeyError, abf.da, 'hiya')
+
+        tj, vj = abf.da(0, join_sweeps=True)
+        self.assertTrue(np.all(tj == np.concatenate(da_times)))
+        self.assertTrue(np.all(vj == np.concatenate(da_values)))
 
         # Conversion to data log without joining
         x = abf.log()
@@ -195,17 +310,32 @@ def test_read_v1(self):
         self.assertIn('6.0.channel', k)
         self.assertIn('7.0.channel', k)
         self.assertIn('8.0.channel', k)
-        self.assertIn('0.1.channel', k)
-        self.assertIn('1.1.channel', k)
-        self.assertIn('2.1.channel', k)
-        self.assertIn('3.1.channel', k)
-        self.assertIn('4.1.channel', k)
-        self.assertIn('5.1.channel', k)
-        self.assertIn('6.1.channel', k)
-        self.assertIn('7.1.channel', k)
-        self.assertIn('8.1.channel', k)
+        self.assertIn('0.0.da', k)
+        self.assertIn('1.0.da', k)
+        self.assertIn('2.0.da', k)
+        self.assertIn('3.0.da', k)
+        self.assertIn('4.0.da', k)
+        self.assertIn('5.0.da', k)
+        self.assertIn('6.0.da', k)
+        self.assertIn('7.0.da', k)
+        self.assertIn('8.0.da', k)
+        self.assertEqual(len(x['time']), len(x['5.0.channel']))
+        self.assertEqual(len(x['7.0.da']), len(x['5.0.channel']))
+
+        x = abf.log(include_da=False)
+        k = list(x.keys())
+        self.assertEqual(len(x), 1 + 1 * 9)
+        self.assertIn('time', x)
+        self.assertIn('0.0.channel', k)
+        self.assertIn('1.0.channel', k)
+        self.assertIn('2.0.channel', k)
+        self.assertIn('3.0.channel', k)
+        self.assertIn('4.0.channel', k)
+        self.assertIn('5.0.channel', k)
+        self.assertIn('6.0.channel', k)
+        self.assertIn('7.0.channel', k)
+        self.assertIn('8.0.channel', k)
         self.assertEqual(len(x['time']), len(x['5.0.channel']))
-        self.assertEqual(len(x['7.1.channel']), len(x['5.0.channel']))
 
         x = abf.log(use_names=True)
         k = list(x.keys())
@@ -230,33 +360,15 @@ def test_read_v1(self):
         self.assertIn('7.OUT 0', k)
         self.assertIn('8.OUT 0', k)
 
-        y = abf.log(use_names=True, channels=[1])
-        k = list(y.keys())
-        self.assertEqual(len(y), 1 + 9)
-        self.assertIn('time', y)
-        self.assertIn('0.OUT 0', y)
-        x = abf.log(use_names=True, channels=['OUT 0'])
-        y = list(x.keys())
-        self.assertEqual(len(y), 1 + 9)
-        self.assertIn('time', y)
-        self.assertIn('8.OUT 0', y)
-        x = abf.log(use_names=True, channels=['IN 0', 'OUT 0'])
-        y = list(x.keys())
-        self.assertEqual(len(y), 1 + 2 * 9)
-        self.assertIn('time', y)
-        self.assertIn('0.OUT 0', y)
-        self.assertIn('2.IN 0', y)
-        self.assertEqual(len(x['time']), len(x['2.IN 0']))
-
         # Conversion to data log with joining
         y = abf.log(join_sweeps=True)
         self.assertEqual(len(y), 3)
         self.assertIn('time', y.keys())
         self.assertIn('0.channel', y.keys())
-        self.assertIn('1.channel', y.keys())
+        self.assertIn('0.da', y.keys())
         self.assertEqual(len(y['time']), 9 * len(x['time']))
         self.assertEqual(len(y['time']), len(y['0.channel']))
-        self.assertEqual(len(y['time']), len(y['1.channel']))
+        self.assertEqual(len(y['time']), len(y['0.da']))
 
         y = abf.log(join_sweeps=True, use_names=True)
         self.assertEqual(len(y), 3)
@@ -267,12 +379,23 @@ def test_read_v1(self):
         self.assertEqual(len(y['time']), len(y['IN 0']))
         self.assertEqual(len(y['time']), len(y['OUT 0']))
 
-        y = abf.log(join_sweeps=True, use_names=True, channels=['IN 0'])
+        y = abf.log(join_sweeps=True, use_names=True, include_da=False)
         self.assertEqual(len(y), 2)
         self.assertIn('time', y.keys())
         self.assertIn('IN 0', y.keys())
+        self.assertEqual(len(y['time']), 9 * len(x['time']))
+        self.assertEqual(len(y['time']), len(y['IN 0']))
+
+        # Test conversion to Myokit protocol
+        p = abf.da_protocol(0)
+        self.assertEqual(len(p), 18)
+        self.assertEqual(p.code(), V1_PROTOCOL)
 
-    def test_read_protocol_v1(self):
+        p = abf.da_protocol(0, vu=myokit.units.V, tu='s')
+        self.assertEqual(len(p), 18)
+        self.assertEqual(p.code(), V1_PROTOCOL_VS)
+
+    def test_read_protocol_file_v1(self):
         # Test reading a v1 protocol file.
 
         # Load file
@@ -280,28 +403,32 @@ def test_read_protocol_v1(self):
         abf = axon.AbfFile(path)
 
         # Check version info
-        self.assertIn('version 1.65', abf.info())
-        self.assertEqual(abf.info(), V1_PROTO_INFO)
-        abf.info(True)  # Test getting full header runs without crashing
-        self.assertEqual(abf.channel_count(), 1)
-        self.assertEqual(abf.ad_channel_count(), 0)
-        self.assertEqual(abf.da_channel_count(), 1)
-        self.assertEqual(abf.channel_names(), ['Cmd 0'])
+        self.assertIn('version 1.65', abf.meta_str())
+        self.assertEqual(abf.meta_str(), V1_PROTO_INFO)
+        abf.meta_str(True)  # Test getting full header runs without crashing
+        self.assertEqual(abf.channel_count(), 0)
+        self.assertEqual(abf.channel_names(), [])
+        self.assertEqual(abf.channel_units(), [])
+        self.assertEqual(abf.da_count(), 1)
+        self.assertEqual(abf.da_names(), ['Cmd 0'])
+        self.assertEqual(abf.da_names(0), 'Cmd 0')
+        self.assertEqual(abf.da_units(), [myokit.units.mV])
+        self.assertEqual(abf.da_units(0), myokit.units.mV)
         #self.assertEqual(len(abf.log(join_sweeps=True)), 2)
 
         # Load, force as protocol
         path = os.path.join(DIR_FORMATS, 'abf-protocol.pro')
         abf = axon.AbfFile(path, is_protocol_file=True)
-        self.assertEqual(abf.channel_count(), 1)
-        self.assertEqual(abf.ad_channel_count(), 0)
-        self.assertEqual(abf.da_channel_count(), 1)
+        self.assertEqual(abf.channel_count(), 0)
+        self.assertEqual(abf.channel_count(), 0)
+        self.assertEqual(abf.da_count(), 1)
 
         # Check version info
-        self.assertIn('version 1.65', abf.info())
-        self.assertIn('Axon Protocol File', abf.info())
+        self.assertIn('version 1.65', abf.meta_str())
+        self.assertIn('Axon Protocol File', abf.meta_str())
 
         # Test protocol extraction
-        p = abf.protocol()
+        p = abf.da_protocol()
         self.assertEqual(len(p), 60)
 
     def test_read_v2(self):
@@ -314,127 +441,184 @@ def test_read_v2(self):
         self.assertEqual(abf.filename(), 'abf-v2.abf')
 
         # Check getting info
-        self.assertIn('version 2.0', abf.info())
+        self.assertIn('version 2.0', abf.meta_str())
         self.maxDiff = None
-        self.assertEqual(abf.info(), V2_INFO)
+        self.assertEqual(abf.meta_str(), V2_INFO)
 
         # Test getting full header runs without crashing
-        abf.info(True)
+        abf.meta_str(True)
 
         # Get version
         self.assertEqual(abf.version(), '2.0.0.0')
 
-        # Test len returns number of sweeps
-        self.assertEqual(len(abf), 1)
+        # Sweep count
+        self.assertEqual(len(abf), 37)
+        self.assertEqual(abf.sweep_count(), 37)
 
         # Test access to A/D channels via native API
-        self.assertEqual(abf.ad_channel_count(), 1)    # 1 data channel
         self.assertIsInstance(abf[0], axon.Sweep)
+        self.assertIsInstance(abf[1], axon.Sweep)
+        self.assertIsInstance(abf[36], axon.Sweep)
+        self.assertFalse(abf[0] is abf[1])
+        self.assertEqual(len(abf[0]), 2)
+        self.assertEqual(len(abf[36]), 2)
         self.assertIsInstance(abf[0][0], axon.Channel)
-        self.assertEqual(abf.da_channel_count(), 1)    # 1 D/A channel
         self.assertIsInstance(abf[0][1], axon.Channel)
-        self.assertEqual(len([s for s in abf]), 1)
-        self.assertEqual(len(abf[0]), abf.channel_count())
-
-        # Test conversion to Myokit protocol
-        p = abf.protocol(1)
-        self.assertEqual(len(p), 2)
-        self.assertEqual(p.code(), V2_PROTOCOL)
-        self.assertRaises(ValueError, abf.protocol, 0)  # Not a D/A channel
-
-        # Test other D/A channel methods
-        self.assertEqual(abf.da_holding_level(1), -120)
-        self.assertEqual(len(abf.da_steps(1)[0]), 1)
-        self.assertEqual(abf.da_steps(1)[0], [12])
+        self.assertFalse(abf[0][0] is abf[0][1])
 
-        # Test Channel methods
+        # Test abf.Channel methods
         channel = abf[0][0]
-        self.assertIsInstance(channel.number(), int)
+        self.assertEqual(channel.index(), 0)
         self.assertEqual(channel.name(), 'IN 0')
+        self.assertIsInstance(channel.unit(), myokit.Unit)
+        self.assertEqual(channel.unit(), myokit.units.pA)
         self.assertEqual(str(channel),
-                         'Channel(0 "IN 0"); 3100 points sampled at 10000.0Hz,'
-                         ' starts at t=76.5501.')
+                         'Channel(0 "IN 0"); 516 points sampled at 20000.0Hz,'
+                         ' starts at t=0.0.')
         self.assertEqual(len(channel.times()), len(channel.values()))
         self.assertFalse(np.all(channel.times() == channel.values()))
 
-        # Test SweepSource interface
-        self.assertEqual(abf.sweep_count(), 1)
-        self.assertEqual(abf.channel_count(), 2)
-        self.assertEqual(abf.channel_names(), ['IN 0', 'Cmd 0'])
+        # Test SweepSource info
+        self.assertTrue(abf.equal_length_sweeps())
+        self.assertEqual(abf.time_unit(), myokit.units.s)
 
-        # Channel without joining
-        x = abf.channel(0)
-        self.assertEqual(len(x), 1 + len(abf))
-        self.assertEqual(len(x[0]), len(x[1]))
+        # Test A/D channel access via SweepSource interface
+        self.assertEqual(abf.channel_count(), 1)
+        self.assertEqual(abf.channel_names(), ['IN 0'])
+        self.assertEqual(abf.channel_names(0), 'IN 0')
+        self.assertEqual(abf.channel_units(), [myokit.units.pA])
+        self.assertEqual(abf.channel_units(0), myokit.units.pA)
+
+        times, values = abf.channel(0)
+        self.assertEqual(37, len(times), len(values))
+        self.assertEqual(len(times[0]), len(values[0]))
+        self.assertEqual(len(times[0]), len(times[1]), len(values[1]))
+        self.assertEqual(len(times[0]), len(times[2]), len(values[2]))
+        self.assertEqual(len(times[0]), len(times[17]), len(values[13]))
+        self.assertEqual(len(times[0]), len(times[36]), len(values[36]))
+
+        tj, vj = abf.channel(0, join_sweeps=True)
+        self.assertTrue(np.all(tj == np.concatenate(times)))
+        self.assertTrue(np.all(vj == np.concatenate(values)))
 
-        # Channel with joining
-        y = abf.channel(0, join_sweeps=True)
-        self.assertEqual(len(y), 2)
-        self.assertEqual(len(y[0]), len(x[0]))
-        self.assertEqual(len(y[0]), len(y[1]))
+        # Channel doesn't exist
+        self.assertRaises(IndexError, abf.channel, -1)
+        self.assertRaises(IndexError, abf.channel, 1)
+        self.assertRaises(KeyError, abf.channel, 'OUT 0')
+
+        # Test D/A output access via SweepSource interface
+        self.assertEqual(abf.da_count(), 1)
+        self.assertEqual(abf.da_names(), ['Cmd 0'])
+        self.assertEqual(abf.da_names(0), 'Cmd 0')
+        self.assertEqual(abf.da_units(), [myokit.units.mV])
+        self.assertEqual(abf.da_units(0), myokit.units.mV)
+
+        da_times, da_values = abf.da(0)
+        self.assertEqual(37, len(da_times), len(da_values))
+        self.assertEqual(len(da_times[0]), len(da_values[0]))
+        self.assertEqual(len(da_times[0]), len(da_times[1]), len(da_values[1]))
+        self.assertEqual(len(da_times[0]), len(da_times[4]), len(da_values[9]))
+        self.assertEqual(len(da_times[0]), len(times[0]))
+        self.assertTrue(np.all(times[0] == da_times[0]))
+        self.assertFalse(np.all(values[0] == da_values[0]))
+
+        tj, vj = abf.da(0, join_sweeps=True)
+        self.assertTrue(np.all(tj == np.concatenate(da_times)))
+        self.assertTrue(np.all(vj == np.concatenate(da_values)))
 
         # Conversion to data log without joining
         x = abf.log()
         k = list(x.keys())
-        self.assertEqual(len(x), 3)
+        self.assertEqual(len(x), 1 + 2 * 37)
         self.assertIn('time', x)
         self.assertIn('0.0.channel', k)
-        self.assertIn('0.1.channel', k)
-        self.assertEqual(len(x['time']), len(x['0.0.channel']))
-        self.assertEqual(len(x['time']), len(x['0.1.channel']))
-        self.assertTrue(np.all(x['time'] == abf.channel(0)[0]))
-        self.assertTrue(np.all(x['0.0.channel'] == abf.channel(0)[1]))
+        self.assertIn('1.0.channel', k)
+        self.assertIn('2.0.channel', k)
+        self.assertIn('3.0.channel', k)
+        self.assertIn('36.0.channel', k)
+        self.assertIn('0.0.da', k)
+        self.assertIn('5.0.da', k)
+        self.assertIn('16.0.da', k)
+        self.assertIn('36.0.da', k)
+        self.assertEqual(len(x['time']), len(x['5.0.channel']))
+        self.assertEqual(len(x['7.0.da']), len(x['5.0.channel']))
+
+        x = abf.log(include_da=False)
+        k = list(x.keys())
+        self.assertEqual(len(x), 1 + 1 * 37)
+        self.assertIn('time', x)
+        self.assertIn('0.0.channel', k)
+        self.assertIn('1.0.channel', k)
+        self.assertIn('2.0.channel', k)
+        self.assertIn('3.0.channel', k)
+        self.assertIn('36.0.channel', k)
+        self.assertNotIn('0.0.da', k)
+        self.assertNotIn('7.0.da', k)
+        self.assertNotIn('36.0.da', k)
+        self.assertEqual(len(x['time']), len(x['6.0.channel']))
 
         x = abf.log(use_names=True)
         k = list(x.keys())
-        self.assertEqual(len(x), 3)
-        self.assertIn('time', k)
+        self.assertEqual(len(x), 1 + 2 * 37)
+        self.assertIn('time', x)
         self.assertIn('0.IN 0', k)
+        self.assertIn('1.IN 0', k)
+        self.assertIn('2.IN 0', k)
+        self.assertIn('36.IN 0', k)
         self.assertIn('0.Cmd 0', k)
-
-        y = abf.log(use_names=True, channels=[1])
-        k = list(y.keys())
-        self.assertEqual(len(y), 2)
-        self.assertIn('time', y)
-        self.assertIn('0.Cmd 0', y)
-        x = abf.log(use_names=True, channels=['Cmd 0'])
-        y = list(x.keys())
-        self.assertEqual(len(y), 2)
-        self.assertIn('time', y)
-        self.assertIn('0.Cmd 0', y)
-        x = abf.log(use_names=True, channels=[0, 'Cmd 0'])
-        y = list(x.keys())
-        self.assertEqual(len(y), 3)
-        self.assertIn('time', y)
-        self.assertIn('0.Cmd 0', y)
-        self.assertIn('0.IN 0', y)
-        self.assertEqual(len(x['time']), len(x['0.IN 0']))
+        self.assertIn('5.Cmd 0', k)
+        self.assertIn('16.Cmd 0', k)
+        self.assertIn('36.Cmd 0', k)
+        self.assertEqual(len(x['time']), len(x['5.IN 0']))
+        self.assertEqual(len(x['7.Cmd 0']), len(x['5.IN 0']))
 
         # Conversion to data log with joining
         y = abf.log(join_sweeps=True)
         self.assertEqual(len(y), 3)
         self.assertIn('time', y.keys())
         self.assertIn('0.channel', y.keys())
-        self.assertIn('1.channel', y.keys())
-        self.assertEqual(len(y['time']), 1 * len(x['time']))
+        self.assertIn('0.da', y.keys())
+        self.assertEqual(len(y['time']), 37 * len(x['time']))
         self.assertEqual(len(y['time']), len(y['0.channel']))
-        self.assertEqual(len(y['time']), len(y['1.channel']))
-        self.assertTrue(np.all(y['time'] == abf.channel(0, True)[0]))
-        self.assertTrue(np.all(y['0.channel'] == abf.channel(0, True)[1]))
+        self.assertEqual(len(y['time']), len(y['0.da']))
 
         y = abf.log(join_sweeps=True, use_names=True)
         self.assertEqual(len(y), 3)
         self.assertIn('time', y.keys())
         self.assertIn('IN 0', y.keys())
         self.assertIn('Cmd 0', y.keys())
-        self.assertEqual(len(y['time']), 1 * len(x['time']))
+        self.assertEqual(len(y['time']), 37 * len(x['time']))
         self.assertEqual(len(y['time']), len(y['IN 0']))
         self.assertEqual(len(y['time']), len(y['Cmd 0']))
 
-        y = abf.log(join_sweeps=True, use_names=True, channels=['IN 0'])
+        y = abf.log(join_sweeps=True, use_names=True, include_da=False)
         self.assertEqual(len(y), 2)
         self.assertIn('time', y.keys())
         self.assertIn('IN 0', y.keys())
+        self.assertEqual(len(y['time']), 37 * len(x['time']))
+        self.assertEqual(len(y['time']), len(y['IN 0']))
+
+        # Test conversion to Myokit protocol
+        p = abf.da_protocol(0)
+        self.assertEqual(len(p), 37 * 2)
+        self.assertEqual(p.code(), V2_PROTOCOL)
+
+        # Test conversion with initial holding as in the real experiment
+        p = abf.da_protocol(0, include_initial_holding=True)
+        e = p.events()
+        self.assertEqual(e[0].level(), -120)    # Pre step
+        self.assertEqual(e[0].start(), 0)
+        self.assertEqual(e[0].duration(), 0.4)
+        self.assertEqual(e[1].level(), -100)    # Real step
+        self.assertEqual(e[1].start(), 0.4)
+        self.assertEqual(e[1].duration(), 25)
+        self.assertEqual(e[2].level(), -120)    # Shortened post-step
+        self.assertEqual(e[2].start(), 25.4)
+        self.assertEqual(e[2].duration(), 4974.6)
+        self.assertEqual(e[3].level(), -120)    # Next pre step
+        self.assertEqual(e[3].start(), 5000)
+        self.assertEqual(e[3].duration(), 0.4)
+        self.assertEqual(len(e), 37 * 3)
 
     def test_matplotlib_figure(self):
         # Test figure drawing method (doesn't inspect output).
@@ -659,7 +843,11 @@ def test_accessors(self):
             self.assertEqual(len(atf), 3)
 
             # Test info
-            self.assertIn('myokit', atf.info())
+            meta = atf.meta_str()
+            self.assertIn('myokit', meta)
+            with WarningCollector() as w:
+                self.assertEqual(meta, atf.info())
+            self.assertIn('deprecated', w.text())
 
             # Test version
             self.assertEqual(atf.version(), '1.0')
diff --git a/myokit/tests/test_formats_wcp.py b/myokit/tests/test_formats_wcp.py
index 95a01afbb..dc0d1bc07 100755
--- a/myokit/tests/test_formats_wcp.py
+++ b/myokit/tests/test_formats_wcp.py
@@ -10,6 +10,7 @@
 
 import numpy as np
 
+import myokit
 import myokit.formats.wcp as wcp
 
 from myokit.tests import DIR_FORMATS, WarningCollector
@@ -28,20 +29,60 @@
 A/D channel: Vm
   Unit: mV
 Records: Type, Status, Sampling Interval, Start, Marker
-Record 0: TEST, ACCEPTED, 0.001, 0.0, ""
-Record 1: TEST, ACCEPTED, 0.001, 0.5, ""
-Record 2: TEST, ACCEPTED, 0.001, 1.0, ""
-Record 3: TEST, ACCEPTED, 0.001, 1.5, ""
-Record 4: TEST, ACCEPTED, 0.001, 2.0615234375, ""
-Record 5: TEST, ACCEPTED, 0.001, 3.0615234375, ""
-Record 6: TEST, ACCEPTED, 0.001, 3.5615234375, ""
-Record 7: TEST, ACCEPTED, 0.001, 4.0615234375, ""
-Record 8: TEST, ACCEPTED, 0.001, 4.5615234375, ""
-Record 9: TEST, ACCEPTED, 0.001, 5.125, ""
-Record 10: TEST, ACCEPTED, 0.001, 5.625, ""
+Record 0: TEST, ACCEPTED, 0.0, ""
+Record 1: TEST, ACCEPTED, 0.5, ""
+Record 2: TEST, ACCEPTED, 1.0, ""
+Record 3: TEST, ACCEPTED, 1.5, ""
+Record 4: TEST, ACCEPTED, 2.0615234375, ""
+Record 5: TEST, ACCEPTED, 3.0615234375, ""
+Record 6: TEST, ACCEPTED, 3.5615234375, ""
+Record 7: TEST, ACCEPTED, 4.0615234375, ""
+Record 8: TEST, ACCEPTED, 4.5615234375, ""
+Record 9: TEST, ACCEPTED, 5.125, ""
+Record 10: TEST, ACCEPTED, 5.625, ""
 '''.strip()
 
 
+INFO_LONG = INFO + '''
+----------------------------------- header -----------------------------------
+ver: 9
+ctime: 21/11/2014 12:43:08
+rtime: 21/11/2014 14:18:28
+nbh: 1024
+adcmax: 32677
+nc: 2
+nba: 2
+nbd: 2
+ad: 10.0
+nr: 11
+dt: 0.001
+nz: 20
+id: ''' + '''
+--------------------------------- raw header ---------------------------------
+rtimesecs: 9062.11
+yo0: 0
+yu0: pA
+yn0: Im
+yg0: 0.0005
+yz0: 0
+yr0: 0
+yo1: 1
+yu1: mV
+yn1: Vm
+yg1: 0.01
+yz1: 0
+yr1: 0
+txperc: 0
+pkpavg: 1
+nsvchan: 0
+nsvalign: 0
+nsvtypr: 0
+nsvs2p: F
+nsvcur0: 0
+nsvcur1: 0
+'''.rstrip()
+
+
 class WcpTest(unittest.TestCase):
     """
     Tests the WCP format module.
@@ -50,7 +91,7 @@ class WcpTest(unittest.TestCase):
     def test_data_file(self):
         # Test basic wcp file reading.
 
-        # Load old file from Maastricht
+        # Load test file
         fname = 'wcp-file.wcp'
         path = os.path.join(DIR_FORMATS, fname)
         w = wcp.WcpFile(path)
@@ -65,7 +106,11 @@ def test_data_file(self):
 
         # Test meta data
         self.maxDiff = None
-        self.assertEqual(w.info(), INFO)
+        self.assertEqual(w.meta_str(), INFO)
+        self.assertEqual(w.meta_str(True), INFO_LONG)
+        with WarningCollector() as c:
+            self.assertEqual(w.info(), INFO)
+        self.assertIn('deprecated', c.text())
 
         # Test Sequence interface
         self.assertEqual(len(w), w.record_count())
@@ -75,19 +120,28 @@ def test_data_file(self):
         # Test SweepSource interface
         self.assertEqual(w.channel_count(), 2)
         self.assertEqual(w.channel_names(), ['Im', 'Vm'])
+        self.assertEqual(w.channel_units(), [myokit.units.pA, myokit.units.mV])
+        self.assertEqual(w.channel_names(0), 'Im')
+        self.assertEqual(w.channel_names(1), 'Vm')
+        self.assertEqual(w.channel_units(0), myokit.units.pA)
+        self.assertEqual(w.channel_units(1), myokit.units.mV)
         self.assertEqual(w.sweep_count(), w.record_count())
+        self.assertTrue(w.equal_length_sweeps())
+        self.assertEqual(w.time_unit(), myokit.units.s)
 
         # Test SweepSource.channel()
         # Without joining
-        x = w.channel(0)
-        y = w.channel(1)
-        self.assertEqual(len(x), len(y), w.sweep_count())
-        self.assertEqual(len(x[0]), w.sample_count())
-        self.assertEqual(len(x[1]), w.sample_count())
-        self.assertEqual(len(x[2]), w.sample_count())
-        self.assertTrue(np.all(x[1] == w.channel('Im')[1]))
-        self.assertTrue(np.all(x[1] != w.channel('Vm')[1]))
-        self.assertTrue(np.all(y[1] == w.channel('Vm')[1]))
+        t0, v0 = w.channel(0)
+        t1, v1 = w.channel(1)
+        self.assertEqual(len(t0), len(t1), w.sweep_count())
+        self.assertEqual(len(v0), len(v1), w.sweep_count())
+        self.assertEqual(len(t0[0]), w.sample_count())
+        self.assertEqual(len(t1[1]), w.sample_count())
+        self.assertEqual(len(v0[2]), w.sample_count())
+        self.assertEqual(len(v1[3]), w.sample_count())
+        self.assertTrue(np.all(v0[0] == w.channel('Im')[1][0]))
+        self.assertTrue(np.all(v0[0] != w.channel('Vm')[1][0]))
+        self.assertTrue(np.all(v1[1] == w.channel('Vm')[1][1]))
 
         # With joining
         x = w.channel(0, join_sweeps=True)
@@ -113,8 +167,8 @@ def test_data_file(self):
             for c in range(w.channel_count()):
                 k.append(f'{r}.{c}.channel')
         self.assertEqual(set(k), set(d.keys()))
-        self.assertTrue(np.all(d['time'] == w.channel(0)[0]))
-        self.assertTrue(np.all(d['1.0.channel'] == w.channel(0)[2]))
+        self.assertTrue(np.all(d['time'] == w.channel(0)[0][0]))
+        self.assertTrue(np.all(d['1.0.channel'] == w.channel(0)[1][1]))
 
         # Without joining, use names
         d = w.log(use_names=True)
@@ -123,8 +177,8 @@ def test_data_file(self):
             for r in range(w.record_count()):
                 k.append(f'{r}.{c}')
         self.assertEqual(set(k), set(d.keys()))
-        self.assertTrue(np.all(d['time'] == w.channel(0)[0]))
-        self.assertTrue(np.all(d['1.Im'] == w.channel(0)[2]))
+        self.assertTrue(np.all(d['time'] == w.channel(0)[0][0]))
+        self.assertTrue(np.all(d['1.Im'] == w.channel(0)[1][1]))
 
         # With joining
         d = w.log(join_sweeps=True)
@@ -138,12 +192,6 @@ def test_data_file(self):
         self.assertTrue(np.all(d['time'] == w.channel(0, True)[0]))
         self.assertTrue(np.all(d['Vm'] == w.channel(1, True)[1]))
 
-        # Selected channels
-        d = w.log(join_sweeps=True, use_names=True, channels=['Im'])
-        self.assertEqual(list(d.keys()), ['time', 'Im'])
-        d = w.log(join_sweeps=True, use_names=True, channels=[1])
-        self.assertEqual(list(d.keys()), ['time', 'Vm'])
-
         # Deprecated methods
         with WarningCollector() as c:
             self.assertEqual(w.records(), w.record_count())
@@ -155,6 +203,13 @@ def test_data_file(self):
             e = w.myokit_log()
         self.assertIn('deprecated', c.text())
 
+        # Unsupported da methods
+        self.assertEqual(w.da_count(), 0)
+        self.assertRaises(NotImplementedError, w.da, 0)
+        self.assertRaises(NotImplementedError, w.da_names)
+        self.assertRaises(NotImplementedError, w.da_units)
+        self.assertRaises(NotImplementedError, w.da_protocol)
+
     def test_figure_method(self):
         # Tests matplotlib_figure
         # Select matplotlib backend that doesn't require a screen