import numbers
import struct
import warnings

import h5py
import numpy as np

import discretisedfield as df
import discretisedfield.plotting as dfp
import discretisedfield.util as dfu
import ubermagutil.typesystem as ts

# TODO: tutorials, line operations


@ts.typesystem(mesh=ts.Typed(expected_type=df.Mesh, const=True),
               dim=ts.Scalar(expected_type=int, positive=True, const=True))
class Field:
    """Finite-difference field.

    This class specifies a finite-difference field and defines operations for
    its analysis and visualisation. The field is defined on a finite-difference
    mesh (`discretisedfield.Mesh`) passed using ``mesh``. Another value that
    must be passed is the dimension of the field's value using ``dim``. For
    instance, for a scalar field, ``dim=1`` and for a three-dimensional vector
    field ``dim=3`` must be passed. The value of the field can be set by
    passing ``value``. For details on how the value can be defined, refer to
    ``discretisedfield.Field.value``. Similarly, if the field has ``dim>1``,
    the field can be normalised by passing ``norm``. For details on setting the
    norm, please refer to ``discretisedfield.Field.norm``.

    Parameters
    ----------
    mesh : discretisedfield.Mesh

        Finite-difference rectangular mesh.

    dim : int

        Dimension of the field's value. For instance, if `dim=3` the field is a
        three-dimensional vector field and for `dim=1` the field is a scalar
        field.

    value : array_like, callable, dict, optional

        Please refer to ``discretisedfield.Field.value`` property. Defaults to
        0, meaning that if the value is not provided in the initialisation,
        "zero-field" will be defined.

    norm : numbers.Real, callable, optional

        Please refer to ``discretisedfield.Field.norm`` property. Defaults to
        ``None`` (``norm=None`` defines no norm).

    Examples
    --------
    1. Defining a uniform three-dimensional vector field on a nano-sized thin
    film.

    >>> import discretisedfield as df
    ...
    >>> p1 = (-50e-9, -25e-9, 0)
    >>> p2 = (50e-9, 25e-9, 5e-9)
    >>> cell = (1e-9, 1e-9, 0.1e-9)
    >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)
    >>> dim = 3
    >>> value = (0, 0, 1)
    ...
    >>> field = df.Field(mesh=mesh, dim=dim, value=value)
    >>> field
    Field(mesh=...)
    >>> field.average
    (0.0, 0.0, 1.0)

    2. Defining a scalar field.

    >>> p1 = (-10, -10, -10)
    >>> p2 = (10, 10, 10)
    >>> n = (1, 1, 1)
    >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
    >>> dim = 1
    >>> value = 3.14
    ...
    >>> field = df.Field(mesh=mesh, dim=dim, value=value)
    >>> field
    Field(mesh=...)
    >>> field.average
    3.14

    3. Defining a uniform three-dimensional normalised vector field.

    >>> import discretisedfield as df
    ...
    >>> p1 = (-50e9, -25e9, 0)
    >>> p2 = (50e9, 25e9, 5e9)
    >>> cell = (1e9, 1e9, 0.1e9)
    >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)
    >>> dim = 3
    >>> value = (0, 0, 8)
    >>> norm = 1
    ...
    >>> field = df.Field(mesh=mesh, dim=dim, value=value, norm=norm)
    >>> field
    Field(mesh=...)
    >>> field.average
    (0.0, 0.0, 1.0)

    .. seealso:: :py:func:`~discretisedfield.Mesh`

    """

    def __init__(self, mesh, dim, value=0, norm=None, components=None):
        self.mesh = mesh
        self.dim = dim
        self.value = value
        self.norm = norm

        self._components = None  # required in here for correct initialisation
        self.components = components

    @property
    def value(self):
        """Field value representation.

        This property returns a representation of the field value if it exists.
        Otherwise, ``discretisedfield.Field.array`` containing all field values
        is returned.

        The value of the field can be set using a scalar value for ``dim=1``
        fields (e.g. ``value=3``) or ``array_like`` value for ``dim>1`` fields
        (e.g. ``value=(1, 2, 3)``). Alternatively, the value can be defined
        using a callable object, which takes a point tuple as an input argument
        and returns a value of appropriate dimension. Internally, callable
        object is called for every point in the mesh on which the field is
        defined. For instance, callable object can be a Python function or
        another ``discretisedfield.Field``. Finally, ``numpy.ndarray`` with
        shape ``(*self.mesh.n, dim)`` can be passed.

        Parameters
        ----------
        value : numbers.Real, array_like, callable, dict

            For scalar fields (``dim=1``) ``numbers.Real`` values are allowed.
            In the case of vector fields, ``array_like`` (list, tuple,
            numpy.ndarray) value with length equal to `dim` should be used.
            Finally, the value can also be a callable (e.g. Python function or
            another field), which for every coordinate in the mesh returns a
            valid value. If ``value=0``, all values in the field will be set to
            zero independent of the field dimension.

            If subregions are defined value can be initialised with a dict.
            Allowed keys are names of all subregions and ``default``. Items
            must be either ``numbers.Real`` for ``dim=1`` or ``array_like``
            for ``dim=3``. If subregion names are missing, the value of
            ``default`` is used if given. If parts of the region are not
            contained within one subregion ``default`` is used if specified,
            else these values are set to 0.

        Returns
        -------
        array_like, callable, numbers.Real, numpy.ndarray

            The value used (representation) for setting the field is returned.
            However, if the actual value of the field does not correspond to
            the initially used value anymore, a ``numpy.ndarray`` is returned
            containing all field values.

        Raises
        ------
        ValueError

            If unsupported type is passed.

        Examples
        --------
        1. Different ways of setting and getting the field value.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> value = (0, 0, 1)
        ...
        >>> # if value is not specified, zero-field is defined
        >>> field = df.Field(mesh=mesh, dim=3)
        >>> field.value
        0
        >>> field.value = (0, 0, 1)
        >>> field.value
        (0, 0, 1)
        >>> # Setting the field value using a Python function (callable).
        >>> def value_function(point):
        ...     x, y, z = point
        ...     if x <= 1:
        ...         return (0, 0, 1)
        ...     else:
        ...         return (0, 0, -1)
        >>> field.value = value_function
        >>> field.value
        <function value_function at ...>
        >>> # We now change the value of a single cell so that the
        >>> # representation used for initialising field is not valid
        >>> # anymore.
        >>> field.array[0, 0, 0, :] = (0, 0, 0)
        >>> field.value
        array(...)
        >>> field.value.shape
        (2, 2, 1, 3)

        2. Field with subregions in mesh
        >>> import discretisedfield as df
        ...
        >>> p1 = (0,0,0)
        >>> p2 = (2,2,2)
        >>> cell = (1,1,1)
        >>> sub1 = df.Region(p1=(0,0,0), p2=(2,2,1))
        >>> sub2 = df.Region(p1=(0,0,1), p2=(2,2,2))
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell,\
                           subregions={'s1': sub1, 's2': sub2})
        >>> field = df.Field(mesh, dim=1, value={'s1': 1, 's2': 1})
        >>> (field.array == 1).all()
        True
        >>> field = df.Field(mesh, dim=1, value={'s1': 1})
        Traceback (most recent call last):
        ...
        KeyError: ...
        >>> field = df.Field(mesh, dim=1, value={'s1': 2, 'default': 1})
        >>> (field.array == 1).all()
        False
        >>> (field.array == 0).any()
        False
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell, subregions={'s': sub1})
        >>> field = df.Field(mesh, dim=1, value={'s': 1})
        Traceback (most recent call last):
        ...
        KeyError: ...
        >>> field = df.Field(mesh, dim=1, value={'default': 1})
        >>> (field.array == 1).all()
        True

        .. seealso:: :py:func:`~discretisedfield.Field.array`

        """
        value_array = dfu.as_array(self.mesh, self.dim, self._value)
        if np.array_equal(self.array, value_array):
            return self._value
        else:
            return self.array

    @value.setter
    def value(self, val):
        self._value = val
        self.array = dfu.as_array(self.mesh, self.dim, val)

    @property
    def components(self):
        """Vector components of the field."""
        return self._components

    @components.setter
    def components(self, components):
        if components is not None:
            if len(components) != self.dim:
                raise ValueError('Number of components does not match'
                                 f' {self.dim=}.')
            if len(components) != len(set(components)):
                raise ValueError('Components must be unique.')
            for c in components:
                if hasattr(self, c):
                    # redefining component labels is okay.
                    if self._components is None or c not in self._components:
                        raise ValueError(
                            f'Component name {c} is already '
                            'used by a different method/property.')
            self._components = list(components)
        else:
            if 2 <= self.dim <= 3:
                components = ['x', 'y', 'z'][:self.dim]
            elif self.dim > 3:
                warnings.warn(f'Component labels must be specified for '
                              f'{self.dim=} fields to get access to individual'
                              ' vector components.')
            self._components = components

    @property
    def array(self):
        """Field value as ``numpy.ndarray``.

        The shape of the array is ``(*mesh.n, dim)``.

        Parameters
        ----------
        array : numpy.ndarray

            Array with shape ``(*mesh.n, dim)``.

        Returns
        -------
        numpy.ndarray

            Field values array.

        Raises
        ------
        ValueError

            If unsupported type or shape is passed.

        Examples
        --------
        1. Accessing and setting the field array.

        >>> import discretisedfield as df
        >>> import numpy as np
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (1, 1, 1)
        >>> cell = (0.5, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> value = (0, 0, 1)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=value)
        >>> field.array
        array(...)
        >>> field.average
        (0.0, 0.0, 1.0)
        >>> field.array.shape
        (2, 1, 1, 3)
        >>> field.array = np.ones_like(field.array)
        >>> field.array
        array(...)
        >>> field.average
        (1.0, 1.0, 1.0)

        .. seealso:: :py:func:`~discretisedfield.Field.value`

        """
        return self._array

    @array.setter
    def array(self, val):
        self._array = dfu.as_array(self.mesh, self.dim, val)

    @property
    def norm(self):
        """Norm of the field.

        Computes the norm of the field and returns ``discretisedfield.Field``
        with ``dim=1``. Norm of a scalar field is interpreted as an absolute
        value of the field. Alternatively, ``discretisedfield.Field.__abs__``
        can be called for obtaining the norm of the field.

        The field norm can be set by passing ``numbers.Real``,
        ``numpy.ndarray``, or callable. If the field has ``dim=1`` or it
        contains zero values, norm cannot be set and ``ValueError`` is raised.

        Parameters
        ----------
        numbers.Real, numpy.ndarray, callable

            Norm value.

        Returns
        -------
        discretisedfield.Field

            Norm of the field if ``dim>1`` or absolute value for ``dim=1``.

        Raises
        ------
        ValueError

            If the norm is set with wrong type, shape, or value. In addition,
            if the field is scalar (``dim=1``) or the field contains zero
            values.

        Examples
        --------
        1. Manipulating the field norm.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (1, 1, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(0, 0, 1))
        >>> field.norm
        Field(...)
        >>> field.norm.average
        1.0
        >>> field.norm = 2
        >>> field.average
        (0.0, 0.0, 2.0)
        >>> field.value = (1, 0, 0)
        >>> field.norm.average
        1.0
        >>> # An attempt to set the norm for a zero field.
        >>> field.value = 0
        >>> field.average
        (0.0, 0.0, 0.0)
        >>> field.norm = 1
        Traceback (most recent call last):
        ...
        ValueError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.__abs__`

        """
        if self.dim == 1:
            res = abs(self.value)
        else:
            res = np.linalg.norm(self.array, axis=-1)[..., np.newaxis]

        return self.__class__(self.mesh, dim=1, value=res)

    @norm.setter
    def norm(self, val):
        if val is not None:
            if self.dim == 1:
                msg = f'Cannot set norm for field with dim={self.dim}.'
                raise ValueError(msg)

            if not np.all(self.norm.array):
                msg = 'Cannot normalise field with zero values.'
                raise ValueError(msg)

            self.array /= self.norm.array  # normalise to 1
            self.array *= dfu.as_array(self.mesh, dim=1, val=val)

    def __abs__(self):
        """Field norm.

        This is a convenience operator and it returns
        ``discretisedfield.Field.norm``. For details, please refer to
        ``discretisedfield.Field.norm``.

        Returns
        -------
        discretisedfield.Field

            Norm of the field if ``dim>1`` or absolute value for ``dim=1``.

        Examples
        --------
        1. Computing the absolute value of a scalar field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 10, 13)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=1, value=-5)
        >>> abs(field).average
        5.0

        .. seealso:: :py:func:`~discretisedfield.Field.norm`

        """
        return self.norm

    @property
    def zero(self):
        """Zero field.

        This method returns a zero field defined on the same mesh and with the
        same value dimension.

        Returns
        -------
        discretisedfield.Field

            Zero field.

        Examples
        --------
        1. Getting the zero-field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 10, 13)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(3, -1, 1))
        >>> zero_field = field.zero
        >>> zero_field.average
        (0.0, 0.0, 0.0)

        """
        return self.__class__(self.mesh, dim=self.dim, value=0,
                              components=self.components)

    @property
    def orientation(self):
        """Orientation field.

        This method computes the orientation (direction) of a vector field and
        returns ``discretisedfield.Field`` with the same dimension. More
        precisely, at every mesh discretisation cell, the vector is divided by
        its norm, so that a unit vector is obtained. However, if the vector at
        a discretisation cell is a zero-vector, it remains unchanged. In the
        case of a scalar (``dim=1``) field, ``ValueError`` is raised.

        Returns
        -------
        discretisedfield.Field

            Orientation field.

        Raises
        ------
        ValueError

            If the field is has ``dim=1``.

        Examples
        --------
        1. Computing the orientation field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(6, 0, 8))
        >>> field.orientation
        Field(...)
        >>> field.orientation.norm.average
        1.0

        """
        if self.dim == 1:
            msg = (f'Cannot compute orientation field for a '
                   f'dim={self.dim} field.')
            raise ValueError(msg)

        orientation_array = np.divide(self.array,
                                      self.norm.array,
                                      out=np.zeros_like(self.array),
                                      where=(self.norm.array != 0))
        return self.__class__(self.mesh, dim=self.dim, value=orientation_array,
                              components=self.components)

    @property
    def average(self):
        """Field average.

        It computes the average of the field over the entire volume of the
        mesh. It returns a tuple with the length same as the dimension
        (``dim``) of the field.

        Returns
        -------
        tuple

            Field average tuple, whose length equals to the field's dimension.

        Examples
        --------
        1. Computing the vector field average.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 5, 5)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(0, 0, 1))
        >>> field.average
        (0.0, 0.0, 1.0)

        2. Computing the scalar field average.

        >>> field = df.Field(mesh=mesh, dim=1, value=55)
        >>> field.average
        55.0

        """
        return dfu.array2tuple(self.array.mean(axis=(0, 1, 2)))

    def __repr__(self):
        """Representation string.

        Returns
        -------
        str

            Representation string.

        Example
        -------
        1. Getting representation string.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh, dim=1, value=1)
        >>> repr(field)
        "Field(mesh=..., dim=1)"

        """
        frepr = f"Field(mesh={repr(self.mesh)}, dim={self.dim}"
        if self.components:
            frepr += f", components={self.components})"
        else:
            frepr += ")"
        return frepr

    def __call__(self, point):
        r"""Sample the field value at ``point``.

        It returns the value of the field in the discretisation cell to which
        ``point`` belongs to. It returns a tuple, whose length is the same as
        the dimension (``dim``) of the field.

        Parameters
        ----------
        point : (3,) array_like

            The mesh point coordinate :math:`\\mathbf{p} = (p_{x}, p_{y},
            p_{z})`.

        Returns
        -------
        tuple

            A tuple, whose length is the same as the dimension of the field.

        Example
        -------
        1. Sampling the field value.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (20, 20, 20)
        >>> n = (20, 20, 20)
        >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), n=n)
        ...
        >>> field = df.Field(mesh, dim=3, value=(1, 3, 4))
        >>> point = (10, 2, 3)
        >>> field(point)
        (1.0, 3.0, 4.0)

        """
        return dfu.array2tuple(self.array[self.mesh.point2index(point)])

    def __getattr__(self, attr):
        """Extract the component of the vector field.

        This method provides access to individual field components for fields
        with dimension > 1. Component labels are defined in the ``components``
        attribute. For dimension 2 and 3 default values ``'x'``, ``'y'``, and
        ``'z'`` are used if no custom component labels are provided. For fields
        with ``dim>3`` component labels must be specified manually to get
        access to individual vector components.

        Parameters
        ----------
        attr : str

            Vector field component defined in ``components``.

        Returns
        -------
        discretisedfield.Field

            Scalar field with vector field component values.

        Examples
        --------
        1. Accessing the default vector field components.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 2)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(0, 0, 1))
        >>> field.x
        Field(...)
        >>> field.x.average
        0.0
        >>> field.y
        Field(...)
        >>> field.y.average
        0.0
        >>> field.z
        Field(...)
        >>> field.z.average
        1.0
        >>> field.z.dim
        1

        2. Accessing custom vector field components.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 2)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh=mesh, dim=3, value=(0, 0, 1),
        ...                  components=['mx', 'my', 'mz'])
        >>> field.mx
        Field(...)
        >>> field.mx.average
        0.0
        >>> field.my
        Field(...)
        >>> field.my.average
        0.0
        >>> field.mz
        Field(...)
        >>> field.mz.average
        1.0
        >>> field.mz.dim
        1

        """
        if self.components is not None and attr in self.components:
            attr_array = self.array[..., self.components.index(attr),
                                    np.newaxis]
            return self.__class__(mesh=self.mesh, dim=1, value=attr_array)
        else:
            msg = f'Object has no attribute {attr}.'
            raise AttributeError(msg)

    def __dir__(self):
        """Extension of the ``dir(self)`` list.

        Adds component labels to the ``dir(self)`` list. Similarly, adds or
        removes methods (``grad``, ``div``,...) depending on the dimension of
        the field.

        Returns
        -------
        list

            Avalilable attributes.

        """
        dirlist = dir(self.__class__)

        if self.components is not None:
            dirlist += self.components
        if self.dim == 1:
            need_removing = ['div', 'curl', 'orientation']
        if self.dim == 2:
            need_removing = ['grad', 'curl', 'k3d']
        if self.dim == 3:
            need_removing = ['grad']

        for attr in need_removing:
            dirlist.remove(attr)

        return dirlist

    def __iter__(self):
        r"""Generator yielding coordinates and values of all mesh
        discretisation cells.

        Yields
        ------
        tuple (2,)

            The first value is the mesh cell coordinates :math:`\\mathbf{p} =
            (p_{x}, p_{y}, p_{z})`, whereas the second one is the field value.

        Examples
        --------
        1. Iterating through the field coordinates and values

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> field = df.Field(mesh, dim=3, value=(0, 0, 1))
        >>> for coord, value in field:
        ...     print (coord, value)
        (0.5, 0.5, 0.5) (0.0, 0.0, 1.0)
        (1.5, 0.5, 0.5) (0.0, 0.0, 1.0)
        (0.5, 1.5, 0.5) (0.0, 0.0, 1.0)
        (1.5, 1.5, 0.5) (0.0, 0.0, 1.0)

        .. seealso:: :py:func:`~discretisedfield.Mesh.indices`

        """
        for point in self.mesh:
            yield point, self(point)

    def __eq__(self, other):
        """Relational operator ``==``.

        Two fields are considered to be equal if:

          1. They are defined on the same mesh.

          2. They have the same dimension (``dim``).

          3. They both contain the same values in ``array``.

        Parameters
        ----------
        other : discretisedfield.Field

            Second operand.

        Returns
        -------
        bool

            ``True`` if two fields are equal, ``False`` otherwise.

        Examples
        --------
        1. Check if two fields are (not) equal.

        >>> import discretisedfield as df
        ...
        >>> mesh = df.Mesh(p1=(0, 0, 0), p2=(5, 5, 5), cell=(1, 1, 1))
        ...
        >>> f1 = df.Field(mesh, dim=1, value=3)
        >>> f2 = df.Field(mesh, dim=1, value=4-1)
        >>> f3 = df.Field(mesh, dim=3, value=(1, 4, 3))
        >>> f1 == f2
        True
        >>> f1 != f2
        False
        >>> f1 == f3
        False
        >>> f1 != f3
        True
        >>> f2 == f3
        False
        >>> f1 == 'a'
        False

        """
        if not isinstance(other, self.__class__):
            return False
        elif (self.mesh == other.mesh and self.dim == other.dim and
              np.array_equal(self.array, other.array)):
            return True
        else:
            return False

    # TODO The mesh comparison has no tolerance.
    def allclose(self, other, rtol=1e-5, atol=1e-8):
        """Allclose method.

        This method determines whether two fields are:

          1. Defined on the same mesh.

          2. Have the same dimension (``dim``).

          3. All values in are within relative (``rtol``) and absolute
          (``atol``) tolerances.

        Parameters
        ----------
        other : discretisedfield.Field

            Field to be compared to.

        rtol : numbers.Real

            Relative tolerance. Defaults to 1e-5.

        atol : numbers.Real

            Absolute tolerance. Defaults to 1e-8.

        Returns
        -------
        bool

            ``True`` if two fields are within tolerance, ``False`` otherwise.

        Raises
        ------
        TypeError

            If a non field object is passed.

        Examples
        --------
        1. Check if two fields are within a tolerance.

        >>> import discretisedfield as df
        ...
        >>> mesh = df.Mesh(p1=(0, 0, 0), p2=(5, 5, 5), cell=(1, 1, 1))
        ...
        >>> f1 = df.Field(mesh, dim=1, value=3)
        >>> f2 = df.Field(mesh, dim=1, value=3+1e-9)
        >>> f3 = df.Field(mesh, dim=1, value=3.1)
        >>> f1.allclose(f2)
        True
        >>> f1.allclose(f3)
        False
        >>> f1.allclose(f3, atol=1e-2)
        False

        """
        if not isinstance(other, self.__class__):
            msg = (f'Cannot apply allclose method between '
                   f'{type(self)=} and {type(other)=} objects.')
            raise TypeError(msg)

        if (self.mesh == other.mesh and self.dim == other.dim):
            return np.allclose(self.array, other.array, rtol=rtol, atol=atol)
        else:
            return False

    def __pos__(self):
        """Unary ``+`` operator.

        This method defines the unary operator ``+``. It returns the field
        itself:

        .. math::

            +f(x, y, z) = f(x, y, z)

        Returns
        -------
        discretisedfield.Field

            Field itself.

        Example
        -------
        1. Applying unary ``+`` operator on a field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5e-9, 5e-9, 5e-9)
        >>> n = (10, 10, 10)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        ...
        >>> f = df.Field(mesh, dim=3, value=(0, -1000, -3))
        >>> res = +f
        >>> res.average
        (0.0, -1000.0, -3.0)
        >>> res == f
        True
        >>> +(+f) == f
        True

        """
        return self

    def __neg__(self):
        r"""Unary ``-`` operator.

        This method negates the value of each discretisation cell. It is
        equivalent to multiplication with -1:

        .. math::

            -f(x, y, z) = -1 \\cdot f(x, y, z)

        Returns
        -------
        discretisedfield.Field

            Field multiplied with -1.

        Example
        -------
        1. Applying unary ``-`` operator on a scalar field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5e-9, 3e-9, 1e-9)
        >>> n = (10, 5, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        ...
        >>> f = df.Field(mesh, dim=1, value=3.1)
        >>> res = -f
        >>> res.average
        -3.1
        >>> f == -(-f)
        True

        2. Applying unary negation operator on a vector field.

        >>> f = df.Field(mesh, dim=3, value=(0, -1000, -3))
        >>> res = -f
        >>> res.average
        (0.0, 1000.0, 3.0)

        """
        return -1 * self

    def __pow__(self, other):
        """Unary ``**`` operator.

        This method defines the ``**`` operator for scalar (``dim=1``) fields
        only. This operator is not defined for vector (``dim>1``) fields, and
        ``ValueError`` is raised.

        Parameters
        ----------
        other : numbers.Real

            Value to which the field is raised.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Applying unary ``**`` operator on a scalar field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (-25e-3, -25e-3, -25e-3)
        >>> p2 = (25e-3, 25e-3, 25e-3)
        >>> n = (10, 10, 10)
        >>> mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), n=n)
        ...
        >>> f = df.Field(mesh, dim=1, value=2)
        >>> res = f**(-1)
        >>> res
        Field(...)
        >>> res.average
        0.5
        >>> res = f**2
        >>> res.average
        4.0
        >>> f**f  # the power must be numbers.Real
        Traceback (most recent call last):
        ...
        TypeError: ...

        2. Attempt to apply power operator on a vector field.

        >>> p1 = (0, 0, 0)
        >>> p2 = (5e-9, 5e-9, 5e-9)
        >>> n = (10, 10, 10)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        ...
        >>> f = df.Field(mesh, dim=3, value=(0, -1, -3))
        >>> f**2
        Traceback (most recent call last):
        ...
        ValueError: ...

        """
        if self.dim != 1:
            msg = f'Cannot apply ** operator on {self.dim=} field.'
            raise ValueError(msg)
        if not isinstance(other, numbers.Real):
            msg = (f'Unsupported operand type(s) for **: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        return self.__class__(self.mesh, dim=1,
                              value=np.power(self.array, other),
                              components=self.components)

    def __add__(self, other):
        """Binary ``+`` operator.

        It can be applied between two ``discretisedfield.Field`` objects or
        between a ``discretisedfield.Field`` object and a "constant". For
        instance if the field is a scalar field, a scalar field or
        ``numbers.Real`` can be the second operand. Similarly, for a vector
        field, either vector field or an iterable, such as ``tuple``, ``list``,
        or ``numpy.ndarray``, can be the second operand. If the second operand
        is a ``discretisedfield.Field`` object, both must be defined on the
        same mesh and have the same dimensions.

        Parameters
        ----------
        other : discretisedfield.Field, numbers.Real, tuple, list, np.ndarray

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Add vector fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 3, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=3, value=(0, -1, -3.1))
        >>> f2 = df.Field(mesh, dim=3, value=(0, 1, 3.1))
        >>> res = f1 + f2
        >>> res.average
        (0.0, 0.0, 0.0)
        >>> f1 + f2 == f2 + f1
        True
        >>> res = f1 + (1, 2, 3.1)
        >>> res.average
        (1.0, 1.0, 0.0)
        >>> f1 + 5
        Traceback (most recent call last):
        ...
        TypeError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.__sub__`

        """
        if isinstance(other, self.__class__):
            if self.dim != other.dim:
                msg = (f'Cannot apply operator + on {self.dim=} '
                       f'and {other.dim=} fields.')
                raise ValueError(msg)
            if self.mesh != other.mesh:
                msg = ('Cannot apply operator + on fields '
                       'defined on different meshes.')
                raise ValueError(msg)
        elif self.dim == 1 and isinstance(other, numbers.Real):
            return self + self.__class__(self.mesh, dim=self.dim, value=other)
        elif self.dim == 3 and isinstance(other, (tuple, list, np.ndarray)):
            return self + self.__class__(self.mesh, dim=self.dim, value=other)
        else:
            msg = (f'Unsupported operand type(s) for +: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        return self.__class__(self.mesh, dim=self.dim,
                              value=self.array + other.array,
                              components=self.components)

    def __radd__(self, other):
        return self + other

    def __sub__(self, other):
        """Binary ``-`` operator.

        It can be applied between two ``discretisedfield.Field`` objects or
        between a ``discretisedfield.Field`` object and a "constant". For
        instance if the field is a scalar field, a scalar field or
        ``numbers.Real`` can be the second operand. Similarly, for a vector
        field, either vector field or an iterable, such as ``tuple``, ``list``,
        or ``numpy.ndarray``, can be the second operand. If the second operand
        is a ``discretisedfield.Field`` object, both must be defined on the
        same mesh and have the same dimensions.

        Parameters
        ----------
        other : discretisedfield.Field, numbers.Real, tuple, list, np.ndarray

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Subtract vector fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 3, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=3, value=(0, 1, 6))
        >>> f2 = df.Field(mesh, dim=3, value=(0, 1, 3))
        >>> res = f1 - f2
        >>> res.average
        (0.0, 0.0, 3.0)
        >>> f1 - f2 == -(f2 - f1)
        True
        >>> res = f1 - (0, 1, 0)
        >>> res.average
        (0.0, 0.0, 6.0)

        .. seealso:: :py:func:`~discretisedfield.Field.__add__`

        """
        # Ensure unary '-' can be applied to other.
        if isinstance(other, (list, tuple)):
            other = np.array(other)

        return self + (-other)

    def __rsub__(self, other):
        return -self + other

    def __mul__(self, other):
        """Binary ``*`` operator.

        It can be applied between:

        1. Two scalar (``dim=1``) fields,

        2. A field of any dimension and ``numbers.Real``,

        3. A field of any dimension and a scalar (``dim=1``) field, or

        4. A field and an "abstract" integration variable (e.g. ``df.dV``)

        If both operands are ``discretisedfield.Field`` objects, they must be
        defined on the same mesh.

        Parameters
        ----------
        other : discretisedfield.Field, numbers.Real

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Multiply two scalar fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=1, value=5)
        >>> f2 = df.Field(mesh, dim=1, value=9)
        >>> res = f1 * f2
        >>> res.average
        45.0
        >>> f1 * f2 == f2 * f1
        True

        2. Multiply vector field with a scalar.

        >>> f1 = df.Field(mesh, dim=3, value=(0, 2, 5))
        ...
        >>> res = f1 * 5  # discretisedfield.Field.__mul__ is called
        >>> res.average
        (0.0, 10.0, 25.0)
        >>> res = 10 * f1  # discretisedfield.Field.__rmul__ is called
        >>> res.average
        (0.0, 20.0, 50.0)

        .. seealso:: :py:func:`~discretisedfield.Field.__truediv__`

        """
        if isinstance(other, self.__class__):
            if self.dim == 3 and other.dim == 3:
                msg = (f'Cannot apply operator * on {self.dim=} '
                       f'and {other.dim=} fields.')
                raise ValueError(msg)
            if self.mesh != other.mesh:
                msg = ('Cannot apply operator * on fields '
                       'defined on different meshes.')
                raise ValueError(msg)
        elif isinstance(other, numbers.Complex):
            return self * self.__class__(self.mesh, dim=1, value=other)
        elif self.dim == 1 and isinstance(other, (tuple, list, np.ndarray)):
            return self * self.__class__(self.mesh,
                                         dim=np.array(other).shape[-1],
                                         value=other)
        elif isinstance(other, df.DValue):
            return self * other(self)
        else:
            msg = (f'Unsupported operand type(s) for *: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        res_array = np.multiply(self.array, other.array)
        components = self.components if self.dim == res_array.shape[-1] \
            else None
        return self.__class__(self.mesh, dim=res_array.shape[-1],
                              value=res_array,
                              components=components)

    def __rmul__(self, other):
        return self * other

    def __truediv__(self, other):
        """Binary ``/`` operator.

        It can be applied between:

        1. Two scalar (``dim=1``) fields,

        2. A field of any dimension and ``numbers.Real``, or

        3. A field of any dimension and a scalar (``dim=1``) field.

        If both operands are ``discretisedfield.Field`` objects, they must be
        defined on the same mesh.

        Parameters
        ----------
        other : discretisedfield.Field, numbers.Real

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Divide two scalar fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=1, value=100)
        >>> f2 = df.Field(mesh, dim=1, value=20)
        >>> res = f1 / f2
        >>> res.average
        5.0
        >>> f1 / f2 == (f2 / f1)**(-1)
        True

        2. Divide vector field by a scalar.

        >>> f1 = df.Field(mesh, dim=3, value=(0, 10, 5))
        >>> res = f1 / 5  # discretisedfield.Field.__mul__ is called
        >>> res.average
        (0.0, 2.0, 1.0)
        >>> 10 / f1  # division by a vector is not allowed
        Traceback (most recent call last):
        ...
        ValueError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.__mul__`

        """
        return self * other**(-1)

    def __rtruediv__(self, other):
        return self**(-1) * other

    def __matmul__(self, other):
        """Binary ``@`` operator, defined as dot product.

        This method computes the dot product between two fields. Both fields
        must be three-dimensional (``dim=3``) and defined on the same mesh.

        Parameters
        ----------
        other : discretisedfield.Field

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Compute the dot product of two vector fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 10e-9, 10e-9)
        >>> cell = (2e-9, 2e-9, 2e-9)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=3, value=(1, 3, 6))
        >>> f2 = df.Field(mesh, dim=3, value=(-1, -2, 2))
        >>> (f1@f2).average
        5.0

        """
        if isinstance(other, self.__class__):
            if self.mesh != other.mesh:
                msg = ('Cannot apply operator @ on fields '
                       'defined on different meshes.')
                raise ValueError(msg)
            if self.dim != 3 or other.dim != 3:
                msg = (f'Cannot apply operator @ on {self.dim=} '
                       f'and {other.dim=} fields.')
                raise ValueError(msg)
        elif isinstance(other, (tuple, list, np.ndarray)):
            return self @ self.__class__(self.mesh, dim=3, value=other,
                                         components=self.components)
        elif isinstance(other, df.DValue):
            return self @ other(self)
        else:
            msg = (f'Unsupported operand type(s) for @: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        res_array = np.einsum('ijkl,ijkl->ijk', self.array, other.array)
        return df.Field(self.mesh, dim=1, value=res_array[..., np.newaxis])

    def __rmatmul__(self, other):
        return self @ other

    def __and__(self, other):
        """Binary ``&`` operator, defined as cross product.

        This method computes the cross product between two fields. Both fields
        must be three-dimensional (``dim=3``) and defined on the same mesh.

        Parameters
        ----------
        other : discretisedfield.Field, tuple, list, numpy.ndarray

            Second operand.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError, TypeError

            If the operator cannot be applied.

        Example
        -------
        1. Compute the cross product of two vector fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=3, value=(1, 0, 0))
        >>> f2 = df.Field(mesh, dim=3, value=(0, 1, 0))
        >>> (f1 & f2).average
        (0.0, 0.0, 1.0)
        >>> (f1 & (0, 0, 1)).average
        (0.0, -1.0, 0.0)

        """
        if isinstance(other, self.__class__):
            if self.mesh != other.mesh:
                msg = ('Cannot apply operator & on fields '
                       'defined on different meshes.')
                raise ValueError(msg)
            if self.dim != 3 or other.dim != 3:
                msg = (f'Cannot apply operator & on {self.dim=} '
                       f'and {other.dim=} fields.')
                raise ValueError(msg)
        elif isinstance(other, (tuple, list, np.ndarray)):
            return self & self.__class__(self.mesh, dim=3, value=other,
                                         components=self.components)
        else:
            msg = (f'Unsupported operand type(s) for &: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        res_array = np.cross(self.array, other.array)
        return self.__class__(self.mesh, dim=3, value=res_array,
                              components=self.components)

    def __rand__(self, other):
        return self & other

    def __lshift__(self, other):
        """Stacks multiple scalar fields in a single vector field.

        This method takes a list of scalar (``dim=1``) fields and returns a
        vector field, whose components are defined by the scalar fields passed.
        If any of the fields passed has ``dim!=1`` or they are not defined on
        the same mesh, an exception is raised. The dimension of the resulting
        field is equal to the length of the passed list.

        Parameters
        ----------
        fields : list

            List of ``discretisedfield.Field`` objects with ``dim=1``.

        Returns
        -------
        disrectisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError

            If the dimension of any of the fields is not 1, or the fields
            passed are not defined on the same mesh.

        Example
        -------
        1. Stack 3 scalar fields.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f1 = df.Field(mesh, dim=1, value=1)
        >>> f2 = df.Field(mesh, dim=1, value=5)
        >>> f3 = df.Field(mesh, dim=1, value=-3)
        ...
        >>> f = f1 << f2 << f3
        >>> f.average
        (1.0, 5.0, -3.0)
        >>> f.dim
        3
        >>> f.x == f1
        True
        >>> f.y == f2
        True
        >>> f.z == f3
        True

        """
        if isinstance(other, self.__class__):
            if self.mesh != other.mesh:
                msg = ('Cannot apply operator << on fields '
                       'defined on different meshes.')
                raise ValueError(msg)
        elif isinstance(other, numbers.Real):
            return self << self.__class__(self.mesh, dim=1, value=other)
        elif isinstance(other, (tuple, list, np.ndarray)):
            return self << self.__class__(self.mesh, dim=len(other),
                                          value=other)
        else:
            msg = (f'Unsupported operand type(s) for <<: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

        array_list = [self.array[..., i] for i in range(self.dim)]
        array_list += [other.array[..., i] for i in range(other.dim)]

        if self.components is None or other.components is None:
            components = None
        else:
            components = self.components + other.components
            if len(components) != len(set(components)):
                # Component name duplicated; could happen e.g. for lshift with
                # a number -> choose labels automatically
                components = None

        return self.__class__(self.mesh, dim=len(array_list),
                              value=np.stack(array_list, axis=3),
                              components=components)

    def __rlshift__(self, other):
        if isinstance(other, numbers.Real):
            return self.__class__(self.mesh, dim=1, value=other) << self
        elif isinstance(other, (tuple, list, np.ndarray)):
            return self.__class__(self.mesh, dim=len(other),
                                  value=other) << self
        else:
            msg = (f'Unsupported operand type(s) for <<: '
                   f'{type(self)=} and {type(other)=}.')
            raise TypeError(msg)

    def pad(self, pad_width, mode, **kwargs):
        """Field padding.

        This method pads the field by adding more cells in chosen direction and
        assigning to them the values as specified by the ``mode`` argument.
        The way in which the field is going to padded is defined by passing
        ``pad_width`` dictionary. The keys of the dictionary are the directions
        (axes), e.g. ``'x'``, ``'y'``, or ``'z'``, whereas the values are the
        tuples of length 2. The first integer in the tuple is the number of
        cells added in the negative direction, and the second integer is the
        number of cells added in the positive direction.

        This method accepts any other arguments allowed by ``numpy.pad``
        function.

        Parameters
        ----------
        pad_width : dict

            The keys of the dictionary are the directions (axes), e.g. ``'x'``,
            ``'y'``, or ``'z'``, whereas the values are the tuples of length 2.
            The first integer in the tuple is the number of cells added in the
            negative direction, and the second integer is the number of cells
            added in the positive direction.

        mode: str

            Padding mode as defined in ``numpy.pad``.

        Returns
        -------
        discretisedfield.Field

            Padded field.

        Examples
        --------
        1. Padding a field in the x direction by 1 cell with ``constant`` mode.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 1, 1)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> field = df.Field(mesh, dim=1, value=1)
        ...
        >>> # Two cells with value 1
        >>> pf = field.pad({'x': (1, 1)}, mode='constant')  # zeros padded
        >>> pf.average
        0.5

        """
        d = {}
        for key, value in pad_width.items():
            d[dfu.axesdict[key]] = value
        padding_sequence = dfu.assemble_index((0, 0), len(self.array.shape), d)

        padded_array = np.pad(self.array, padding_sequence,
                              mode=mode, **kwargs)
        padded_mesh = self.mesh.pad(pad_width)

        return self.__class__(padded_mesh, dim=self.dim, value=padded_array,
                              components=self.components)

    def derivative(self, direction, n=1):
        """Directional derivative.

        This method computes a directional derivative of the field and returns
        a field. The direction in which the derivative is computed is passed
        via ``direction`` argument, which can be ``'x'``, ``'y'``, or ``'z'``.
        The order of the computed derivative can be 1 or 2 and it is specified
        using argument ``n`` and it defaults to 1.

        Directional derivative cannot be computed if only one discretisation
        cell exists in a specified direction. In that case, a zero field is
        returned. More precisely, it is assumed that the field does not change
        in that direction. Computing of the directional derivative depends
        strongly on the boundary condition specified in the mesh on which the
        field is defined on. More precisely, the values of the derivatives at
        the boundary are different for periodic, Neumann, or no boundary
        conditions. For details on boundary conditions, please refer to the
        ``disretisedfield.Mesh`` class. The derivatives are computed using
        central differences inside the sample and using forward/backward
        differences at the boundaries.

        Parameters
        ----------
        direction : str

            The direction in which the derivative is computed. It can be
            ``'x'``, ``'y'``, or ``'z'``.

        n : int

            The order of the derivative. It can be 1 or 2 and it defaults to 1.

        Returns
        -------
        discretisedfield.Field

            Directional derivative.

        Raises
        ------
        NotImplementedError

            If order ``n`` higher than 2 is asked for.

        Example
        -------
        1. Compute the first-order directional derivative of a scalar field in
        the y-direction of a spatially varying field. For the field we choose
        :math:`f(x, y, z) = 2x + 3y - 5z`. Accordingly, we expect the
        derivative in the y-direction to be to be a constant scalar field
        :math:`df/dy = 3`.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (100e-9, 100e-9, 10e-9)
        >>> cell = (10e-9, 10e-9, 10e-9)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return 2*x + 3*y + -5*z
        ...
        >>> f = df.Field(mesh, dim=1, value=value_fun)
        >>> f.derivative('y').average  # first-order derivative by default
        3.0

        2. Try to compute the second-order directional derivative of the vector
        field which has only one discretisation cell in the z-direction. For
        the field we choose :math:`f(x, y, z) = (2x, 3y, -5z)`. Accordingly, we
        expect the directional derivatives to be: :math:`df/dx = (2, 0, 0)`,
        :math:`df/dy=(0, 3, 0)`, :math:`df/dz = (0, 0, -5)`. However, because
        there is only one discretisation cell in the z-direction, the
        derivative cannot be computed and a zero field is returned. Similarly,
        second-order derivatives in all directions are expected to be zero.

        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return (2*x, 3*y, -5*z)
        ...
        >>> f = df.Field(mesh, dim=3, value=value_fun)
        >>> f.derivative('x', n=1).average
        (2.0, 0.0, 0.0)
        >>> f.derivative('y', n=1).average
        (0.0, 3.0, 0.0)
        >>> f.derivative('z', n=1).average  # derivative cannot be calculated
        (0.0, 0.0, 0.0)
        >>> # second-order derivatives

        """
        direction = dfu.axesdict[direction]

        # If there are no neighbouring cells in the specified direction, zero
        # field is returned.
        if self.mesh.n[direction] == 1:
            return self.zero

        # Preparation (padding) for computing the derivative, depending on the
        # boundary conditions (PBC, Neumann, or no BC). Depending on the BC,
        # the field array is padded.
        if dfu.raxesdict[direction] in self.mesh.bc:  # PBC
            pad_width = {dfu.raxesdict[direction]: (1, 1)}
            padding_mode = 'wrap'
        elif self.mesh.bc == 'neumann':
            pad_width = {dfu.raxesdict[direction]: (1, 1)}
            padding_mode = 'edge'
        else:  # No BC - no padding
            pad_width = {}
            padding_mode = 'constant'

        padded_array = self.pad(pad_width, mode=padding_mode).array

        if n not in (1, 2):
            msg = f'Derivative of the n={n} order is not implemented.'
            raise NotImplementedError(msg)

        elif n == 1:
            if self.dim == 1:
                derivative_array = np.gradient(padded_array[..., 0],
                                               self.mesh.cell[direction],
                                               axis=direction)[..., np.newaxis]
            else:
                derivative_array = np.gradient(padded_array,
                                               self.mesh.cell[direction],
                                               axis=direction)

        elif n == 2:
            derivative_array = np.zeros_like(padded_array)
            for i in range(padded_array.shape[direction]):
                if i == 0:
                    i1, i2, i3 = i+2, i+1, i
                elif i == padded_array.shape[direction] - 1:
                    i1, i2, i3 = i, i-1, i-2
                else:
                    i1, i2, i3 = i+1, i, i-1
                index1 = dfu.assemble_index(slice(None), 4, {direction: i1})
                index2 = dfu.assemble_index(slice(None), 4, {direction: i2})
                index3 = dfu.assemble_index(slice(None), 4, {direction: i3})
                index = dfu.assemble_index(slice(None), 4, {direction: i})
                derivative_array[index] = ((padded_array[index1] -
                                           2*padded_array[index2] +
                                           padded_array[index3]) /
                                           self.mesh.cell[direction]**2)

        # Remove padded values (if any).
        if derivative_array.shape != self.array.shape:
            derivative_array = np.delete(derivative_array,
                                         (0, self.mesh.n[direction]+1),
                                         axis=direction)

        return self.__class__(self.mesh, dim=self.dim, value=derivative_array,
                              components=self.components)

    @property
    def grad(self):
        r"""Gradient.

        This method computes the gradient of a scalar (``dim=1``) field and
        returns a vector field:

        .. math::

            \\nabla f = (\\frac{\\partial f}{\\partial x},
                         \\frac{\\partial f}{\\partial y},
                         \\frac{\\partial f}{\\partial z})

        Directional derivative cannot be computed if only one discretisation
        cell exists in a certain direction. In that case, a zero field is
        considered to be that directional derivative. More precisely, it is
        assumed that the field does not change in that direction.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError

            If the dimension of the field is not 1.

        Example
        -------
        1. Compute gradient of a contant field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 10e-9, 10e-9)
        >>> cell = (2e-9, 2e-9, 2e-9)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f = df.Field(mesh, dim=1, value=5)
        >>> f.grad.average
        (0.0, 0.0, 0.0)

        2. Compute gradient of a spatially varying field. For a field we choose
        :math:`f(x, y, z) = 2x + 3y - 5z`. Accordingly, we expect the gradient
        to be a constant vector field :math:`\\nabla f = (2, 3, -5)`.

        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return 2*x + 3*y - 5*z
        ...
        >>> f = df.Field(mesh, dim=1, value=value_fun)
        >>> f.grad.average
        (2.0, 3.0, -5.0)

        3. Attempt to compute the gradient of a vector field.

        >>> f = df.Field(mesh, dim=3, value=(1, 2, -3))
        >>> f.grad
        Traceback (most recent call last):
        ...
        ValueError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.derivative`

        """
        if self.dim != 1:
            msg = f'Cannot compute gradient for dim={self.dim} field.'
            raise ValueError(msg)

        return (self.derivative('x') <<
                self.derivative('y') <<
                self.derivative('z'))

    @property
    def div(self):
        r"""Divergence.

        This method computes the divergence of a vector (``dim=2`` or
        ``dim=3``) field and returns a scalar (``dim=1``) field as a result.

        .. math::

            \\nabla\\cdot\\mathbf{v} = \\sum_i\\frac{\\partial v_{i}}
            {\\partial i}

        Directional derivative cannot be computed if only one discretisation
        cell exists in a certain direction. In that case, a zero field is
        considered to be that directional derivative. More precisely, it is
        assumed that the field does not change in that direction.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError

            If the dimension of the field is 1.

        Example
        -------
        1. Compute the divergence of a vector field. For a field we choose
        :math:`\\mathbf{v}(x, y, z) = (2x, -2y, 5z)`. Accordingly, we expect
        the divergence to be to be a constant scalar field :math:`\\nabla\\cdot
        \\mathbf{v} = 5`.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (100e-9, 100e-9, 100e-9)
        >>> cell = (10e-9, 10e-9, 10e-9)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return (2*x, -2*y, 5*z)
        ...
        >>> f = df.Field(mesh, dim=3, value=value_fun)
        >>> f.div.average
        5.0

        2. Attempt to compute the divergence of a scalar field.

        >>> f = df.Field(mesh, dim=1, value=3.14)
        >>> f.div
        Traceback (most recent call last):
        ...
        ValueError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.derivative`

        """
        if self.dim not in [2, 3]:
            msg = f'Cannot compute divergence for dim={self.dim} field.'
            raise ValueError(msg)

        return sum([getattr(self,
                            self.components[i]).derivative(dfu.raxesdict[i])
                    for i in range(self.dim)])

    @property
    def curl(self):
        r"""Curl.

        This method computes the curl of a vector (``dim=3``) field and returns
        a vector (``dim=3``) as a result:

        .. math::

            \\nabla \\times \\mathbf{v} = \\left(\\frac{\\partial
            v_{z}}{\\partial y} - \\frac{\\partial v_{y}}{\\partial z},
            \\frac{\\partial v_{x}}{\\partial z} - \\frac{\\partial
            v_{z}}{\\partial x}, \\frac{\\partial v_{y}}{\\partial x} -
            \\frac{\\partial v_{x}}{\\partial y},\\right)

        Directional derivative cannot be computed if only one discretisation
        cell exists in a certain direction. In that case, a zero field is
        considered to be that directional derivative. More precisely, it is
        assumed that the field does not change in that direction.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Raises
        ------
        ValueError

            If the dimension of the field is not 3.

        Example
        -------
        1. Compute curl of a vector field. For a field we choose
        :math:`\\mathbf{v}(x, y, z) = (2xy, -2y, 5xz)`. Accordingly, we expect
        the curl to be to be a constant vector field :math:`\\nabla\\times
        \\mathbf{v} = (0, -5z, -2x)`.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return (2*x*y, -2*y, 5*x*z)
        ...
        >>> f = df.Field(mesh, dim=3, value=value_fun)
        >>> f.curl((1, 1, 1))
        (0.0, -5.0, -2.0)

        2. Attempt to compute the curl of a scalar field.

        >>> f = df.Field(mesh, dim=1, value=3.14)
        >>> f.curl
        Traceback (most recent call last):
        ...
        ValueError: ...

        .. seealso:: :py:func:`~discretisedfield.Field.derivative`

        """
        if self.dim != 3:
            msg = f'Cannot compute curl for dim={self.dim} field.'
            raise ValueError(msg)

        x, y, z = self.components
        curl_x = (getattr(self, z).derivative('y')
                  - getattr(self, y).derivative('z'))
        curl_y = (getattr(self, x).derivative('z')
                  - getattr(self, z).derivative('x'))
        curl_z = (getattr(self, y).derivative('x')
                  - getattr(self, x).derivative('y'))

        return curl_x << curl_y << curl_z

    @property
    def laplace(self):
        r"""Laplace operator.

        This method computes the laplacian of a scalar (``dim=1``) or a vector
        (``dim=3``) field and returns a resulting field:

        .. math::

            \\nabla^2 f = \\frac{\\partial^{2} f}{\\partial x^{2}} +
                          \\frac{\\partial^{2} f}{\\partial y^{2}} +
                          \\frac{\\partial^{2} f}{\\partial z^{2}}

        .. math::

            \\nabla^2 \\mathbf{f} = (\\nabla^2 f_{x},
                                     \\nabla^2 f_{y},
                                     \\nabla^2 f_{z})

        Directional derivative cannot be computed if only one discretisation
        cell exists in a certain direction. In that case, a zero field is
        considered to be that directional derivative. More precisely, it is
        assumed that the field does not change in that direction.

        Returns
        -------
        discretisedfield.Field

            Resulting field.

        Example
        -------
        1. Compute Laplacian of a contant scalar field.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 10e-9, 10e-9)
        >>> cell = (2e-9, 2e-9, 2e-9)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f = df.Field(mesh, dim=1, value=5)
        >>> f.laplace.average
        0.0

        2. Compute Laplacian of a spatially varying field. For a field we
        choose :math:`f(x, y, z) = 2x^{2} + 3y - 5z`. Accordingly, we expect
        the Laplacian to be a constant vector field :math:`\\nabla f = (4, 0,
        0)`.

        >>> def value_fun(point):
        ...     x, y, z = point
        ...     return 2*x**2 + 3*y - 5*z
        ...
        >>> f = df.Field(mesh, dim=1, value=value_fun)
        >>> assert abs(f.laplace.average - 4) < 1e-3

        .. seealso:: :py:func:`~discretisedfield.Field.derivative`

        """
        if self.dim not in [1, 3]:
            raise ValueError(
                f'Cannot compute laplace for dim={self.dim} field.')
        if self.dim == 1:
            return (self.derivative('x', n=2) +
                    self.derivative('y', n=2) +
                    self.derivative('z', n=2))
        else:
            x, y, z = self.components
            return (getattr(self, x).laplace
                    << getattr(self, y).laplace
                    << getattr(self, z).laplace)

    def integral(self, direction='xyz', improper=False):
        r"""Integral.

        This method integrates the field over the mesh along the specified
        direction(s), which can ce specified using ``direction``. Field must be
        explicitly multiplied by an infinitesimal value (``DValue``) before
        integration. Improper integral can be computed by passing
        ``improper=True`` and by specifying a single direction.

        Parameters
        ----------
        direction : str, optional

            Direction(s) along which the field is integrated. Defaults to
            ``'xyz'``.

        improper : bool, optional

            If ``True``, an improper (cumulative) integral is computed.
            Defaults to ``False``.

        Returns
        -------
        discretisedfield.Field, numbers.Real, or (3,) array_like

            Integration result. If the field is integrated in all directions,
            ``numbers.Real`` or ``array_like`` value is returned depending on
            the dimension of the field.

        Raises
        ------
        ValueError

            If ``improper=True`` and more than one integration direction is
            specified.

        Example
        -------
        1. Volume integral of a scalar field.

        .. math::

            \\int_\\mathrm{V} f(\\mathbf{r}) \\mathrm{d}V

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10, 10, 10)
        >>> cell = (2, 2, 2)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        ...
        >>> f = df.Field(mesh, dim=1, value=5)
        >>> (f * df.dV).integral()
        5000.0

        2. Volume integral of a vector field.

        .. math::

            \\int_\\mathrm{V} \\mathbf{f}(\\mathbf{r}) \\mathrm{d}V

        >>> f = df.Field(mesh, dim=3, value=(-1, -2, -3))
        >>> (f * df.dV).integral()
        (-1000.0, -2000.0, -3000.0)

        3. Surface integral of a scalar field.

        .. math::

            \\int_\\mathrm{S} f(\\mathbf{r}) |\\mathrm{d}\\mathbf{S}|

        >>> f = df.Field(mesh, dim=1, value=5)
        >>> f_plane = f.plane('z')
        >>> (f_plane * abs(df.dS)).integral()
        500.0

        4. Surface integral of a vector field (flux).

        .. math::

            \\int_\\mathrm{S} \\mathbf{f}(\\mathbf{r}) \\cdot
            \\mathrm{d}\\mathbf{S}

        >>> f = df.Field(mesh, dim=3, value=(1, 2, 3))
        >>> f_plane = f.plane('z')
        >>> (f_plane @ df.dS).integral()
        300.0

        5. Integral along x-direction.

        .. math::

            \\int_{x_\\mathrm{min}}^{x_\\mathrm{max}} \\mathbf{f}(\\mathbf{r})
            \\mathrm{d}x

        >>> f = df.Field(mesh, dim=3, value=(1, 2, 3))
        >>> f_plane = f.plane('z')
        >>> (f_plane * df.dx).integral(direction='x').average
        (10.0, 20.0, 30.0)

        6. Improper integral along x-direction.

        .. math::

            \\int_{x_\\mathrm{min}}^{x} \\mathbf{f}(\\mathbf{r})
            \\mathrm{d}x'

        >>> f = df.Field(mesh, dim=3, value=(1, 2, 3))
        >>> f_plane = f.plane('z')
        >>> (f_plane * df.dx).integral(direction='x', improper=True)
        Field(...)

        """
        if improper and len(direction) > 1:
            msg = 'Cannot compute improper integral along multiple directions.'
            raise ValueError(msg)

        mesh = self.mesh

        if not improper:
            for i in direction:
                mesh = mesh.plane(i)
            axes = [dfu.axesdict[i] for i in direction]
            res_array = np.sum(self.array, axis=tuple(axes), keepdims=True)
        else:
            res_array = np.cumsum(self.array, axis=dfu.axesdict[direction])

        res = self.__class__(mesh, dim=self.dim, value=res_array,
                             components=self.components)

        if len(direction) == 3:
            return dfu.array2tuple(res.array.squeeze())
        else:
            return res

    def line(self, p1, p2, n=100):
        r"""Sample the field along the line.

        Given two points :math:`p_{1}` and :math:`p_{2}`, :math:`n` position
        coordinates are generated and the corresponding field values.

        .. math::

           \\mathbf{r}_{i} = i\\frac{\\mathbf{p}_{2} -
           \\mathbf{p}_{1}}{n-1}

        Parameters
        ----------
        p1, p2 : (3,) array_like

            Two points between which the line is generated.

        n : int, optional

            Number of points on the line. Defaults to 100.

        Returns
        -------
        discretisedfield.Line

            Line object.

        Raises
        ------
        ValueError

            If ``p1`` or ``p2`` is outside the mesh domain.

        Examples
        --------
        1. Sampling the field along the line.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 2)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> field = df.Field(mesh, dim=2, value=(0, 3))
        ...
        >>> line = field.line(p1=(0, 0, 0), p2=(2, 0, 0), n=5)

        """
        points = list(self.mesh.line(p1=p1, p2=p2, n=n))
        values = [self(p) for p in points]

        return df.Line(points=points, values=values)

    def plane(self, *args, n=None, **kwargs):
        """Extracts field on the plane mesh.

        If one of the axes (``'x'``, ``'y'``, or ``'z'``) is passed as a
        string, a plane mesh perpendicular to that axis is extracted,
        intersecting the mesh region at its centre, and the field is sampled on
        that mesh. Alternatively, if a keyword argument is passed (e.g.
        ``x=1e-9``), a plane perpendicular to the x-axis (parallel to yz-plane)
        and intersecting it at ``x=1e-9`` is extracted. The number of points in
        two dimensions on the plane can be defined using ``n`` tuple (e.g.
        ``n=(10, 15)``).

        Parameters
        ----------
        n : (2,) tuple

            The number of points on the plane in two dimensions.

        Returns
        ------
        discretisedfield.Field

            An extracted field.

        Examples
        --------
        1. Extracting the field on a plane at a specific point.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (5, 5, 5)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> f = df.Field(mesh, dim=3, value=(0, 0, 1))
        ...
        >>> f.plane(y=1)
        Field(...)

        2. Extracting the field at the mesh region centre.

        >>> f.plane('z')
        Field(...)

        3. Specifying the number of points.

        >>> f.plane('z', n=(10, 10))
        Field(...)

        .. seealso:: :py:func:`~discretisedfield.Mesh.plane`

        """
        plane_mesh = self.mesh.plane(*args, n=n, **kwargs)
        return self.__class__(plane_mesh, dim=self.dim, value=self,
                              components=self.components)

    def __getitem__(self, item):
        """Extracts the field on a subregion.

        If subregions were defined by passing ``subregions`` dictionary when
        the mesh was created, this method returns a field in a subregion
        ``subregions[item]``. Alternatively, a ``discretisedfield.Region``
        object can be passed and a minimum-sized field containing it will be
        returned. The resulting mesh has the same discretisation cell as the
        original field's mesh.

        Parameters
        ----------
        item : str, discretisedfield.Region

            The key of a subregion in ``subregions`` dictionary or a region
            object.

        Returns
        -------
        disretisedfield.Field

            Field on a subregion.

        Example
        -------
        1. Extract field on the subregion by passing a key.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (100, 100, 100)
        >>> cell = (10, 10, 10)
        >>> subregions = {'r1': df.Region(p1=(0, 0, 0), p2=(50, 100, 100)),
        ...               'r2': df.Region(p1=(50, 0, 0), p2=(100, 100, 100))}
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell, subregions=subregions)
        >>> def value_fun(point):
        ...     x, y, z = point
        ...     if x <= 50:
        ...         return (1, 2, 3)
        ...     else:
        ...         return (-1, -2, -3)
        ...
        >>> f = df.Field(mesh, dim=3, value=value_fun)
        >>> f.average
        (0.0, 0.0, 0.0)
        >>> f['r1']
        Field(...)
        >>> f['r1'].average
        (1.0, 2.0, 3.0)
        >>> f['r2'].average
        (-1.0, -2.0, -3.0)

        2. Extracting a subfield by passing a region.

        >>> import discretisedfield as df
        ...
        >>> p1 = (-50e-9, -25e-9, 0)
        >>> p2 = (50e-9, 25e-9, 5e-9)
        >>> cell = (5e-9, 5e-9, 5e-9)
        >>> region = df.Region(p1=p1, p2=p2)
        >>> mesh = df.Mesh(region=region, cell=cell)
        >>> field = df.Field(mesh=mesh, dim=1, value=5)
        ...
        >>> subregion = df.Region(p1=(-9e-9, -1e-9, 1e-9),
        ...                       p2=(9e-9, 14e-9, 4e-9))
        >>> subfield = field[subregion]
        >>> subfield.array.shape
        (4, 4, 1, 1)

        """
        submesh = self.mesh[item]

        index_min = self.mesh.point2index(submesh.index2point((0, 0, 0)))
        index_max = np.add(index_min, submesh.n)
        slices = [slice(i, j) for i, j in zip(index_min, index_max)]
        return self.__class__(submesh, dim=self.dim,
                              value=self.array[tuple(slices)],
                              components=self.components)

    def project(self, direction):
        """Projects the field along one direction and averages it out along
        that direction.

        One of the axes (``'x'``, ``'y'``, or ``'z'``) is passed and the field
        is projected (averaged) along that direction. For example
        ``project('z')`` would average the field in the z-direction and return
        the field which has only one discretisation cell in the z-direction.

        Parameters
        ----------
        direction : str

            Direction along which the field is projected (``'x'``, ``'y'``, or
            ``'z'``).

        Returns
        ------
        discretisedfield.Field

            A projected field.

        Example
        -------
        1. Projecting the field along a certain direction.

        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (2, 2, 2)
        >>> cell = (1, 1, 1)
        >>> mesh = df.Mesh(p1=p1, p2=p2, cell=cell)
        >>> field = df.Field(mesh, dim=3, value=(1, 2, 3))
        ...
        >>> field.project('z')
        Field(...)
        >>> field.project('z').average
        (1.0, 2.0, 3.0)
        >>> field.project('z').array.shape
        (2, 2, 1, 3)

        """
        n_cells = self.mesh.n[dfu.axesdict[direction]]
        return self.integral(direction=direction) / n_cells

    @property
    def angle(self):
        r"""In-plane angle of the vector field.

        This method can be applied only on sliced fields, when a plane is
        defined. This method then returns a scalar field which is an angle
        between the in-plane compoenent of the vector field and the horizontal
        axis. The angle is computed in radians and all values are in :math:`(0,
        2\\pi)` range.

        Returns
        -------
        discretisedfield.Field

            Angle scalar field.

        Raises
        ------
        ValueError

            If the field is not sliced.

        Example
        -------
        1. Computing the angle of the field in yz-plane.

        >>> import discretisedfield as df
        >>> import numpy as np
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (100, 100, 100)
        >>> n = (10, 10, 10)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        >>> field = df.Field(mesh, dim=3, value=(0, 1, 0))
        ...
        >>> abs(field.plane('z').angle.average - np.pi/2) < 1e-3
        True

        """
        if not self.mesh.attributes['isplane']:
            msg = 'The field must be sliced before angle can be computed.'
            raise ValueError(msg)

        angle_array = np.arctan2(
            self.array[..., self.mesh.attributes['axis2']],
            self.array[..., self.mesh.attributes['axis1']])

        # Place all values in [0, 2pi] range
        angle_array[angle_array < 0] += 2 * np.pi

        return self.__class__(self.mesh, dim=1,
                              value=angle_array[..., np.newaxis])

    def write(self, filename, representation='txt', extend_scalar=False):
        """Write the field to OVF, HDF5, or VTK file.

        If the extension of ``filename`` is ``.vtk``, a VTK file is written
        (:py:func:`~discretisedfield.Field._writevtk`).

        For ``.ovf``, ``.omf``, or ``.ohf`` extensions, the field is saved to
        OVF file (:py:func:`~discretisedfield.Field._writeovf`). In that case,
        the representation of data (``'bin4'``, ``'bin8'``, or ``'txt'``) is
        passed as ``representation`` and if ``extend_scalar=True``, a scalar
        field will be saved as a vector field. More precisely, if the value at
        a cell is X, that cell will be saved as (X, 0, 0).

        Finally, if the extension of ``filename`` is ``.hdf5``, HDF5 file will
        be written (:py:func:`~discretisedfield.Field._writehdf5`).

        Parameters
        ----------
        filename : str

            Name of the file written.

        representation : str, optional

            In the case of OVF files (``.ovf``, ``.omf``, or ``.ohf``),
            representation can be specified (``'bin4'``, ``'bin8'``, or
            ``'txt'``). Defaults to ``'txt'``.

        extend_scalar : bool, optional

            If ``True``, a scalar field will be saved as a vector field. More
            precisely, if the value at a cell is 3, that cell will be saved as
            (3, 0, 0). This is valid only for the OVF file formats. Defaults to
            ``False``.

        Example
        -------
        1. Write field to the OVF file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, -5e-9)
        >>> p2 = (5e-9, 15e-9, 15e-9)
        >>> n = (5, 15, 20)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        >>> field = df.Field(mesh, dim=3, value=(5, 6, 7))
        ...
        >>> filename = 'mytestfile.omf'
        >>> field.write(filename)  # write the file
        >>> os.path.isfile(filename)
        True
        >>> field_read = df.Field.fromfile(filename)  # read the file
        >>> field_read == field
        True
        >>> os.remove(filename)  # delete the file

        2. Write field to the VTK file.

        >>> filename = 'mytestfile.vtk'
        >>> field.write(filename)  # write the file
        >>> os.path.isfile(filename)
        True
        >>> os.remove(filename)  # delete the file

        3. Write field to the HDF5 file.

        >>> filename = 'mytestfile.hdf5'
        >>> field.write(filename)  # write the file
        >>> os.path.isfile(filename)
        True
        >>> field_read = df.Field.fromfile(filename)  # read the file
        >>> field_read == field
        True
        >>> os.remove(filename)  # delete the file

        .. seealso:: :py:func:`~discretisedfield.Field.fromfile`

        """
        if any([filename.endswith(ext) for ext in ['.omf', '.ovf', '.ohf']]):
            self._writeovf(filename, representation=representation,
                           extend_scalar=extend_scalar)
        elif any([filename.endswith(ext) for ext in ['.hdf5', '.h5']]):
            self._writehdf5(filename)
        elif filename.endswith('.vtk'):
            self._writevtk(filename)
        else:
            msg = (f'Writing file with extension {filename.split(".")[-1]} '
                   f'not supported.')
            raise ValueError(msg)

    def _writeovf(self, filename, representation='txt', extend_scalar=False):
        """Write the field to an OVF2.0 file.

        Data representation (``'bin4'``, ``'bin8'``, or ``'txt'``) is passed
        using ``representation`` argument. If ``extend_scalar=True``, a scalar
        field will be saved as a vector field. More precisely, if the value at
        a cell is X, that cell will be saved as (X, 0, 0).

        Parameters
        ----------
        filename : str

            Name with an extension of the file written.

        representation : str, optional

            Representation; ``'bin4'``, ``'bin8'``, or ``'txt'``. Defaults to
            ``'txt'``.

        extend_scalar : bool, optional

            If ``True``, a scalar field will be saved as a vector field. More
            precisely, if the value at a cell is 3, that cell will be saved as
            (3, 0, 0). Defaults to ``False``.

        Example
        -------
        1. Write field to the OVF file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 5e-9, 3e-9)
        >>> n = (10, 5, 3)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        >>> value_fun = lambda point: (point[0], point[1], point[2])
        >>> field = df.Field(mesh, dim=3, value=value_fun)
        ...
        >>> filename = 'mytestfile.ohf'
        >>> field._writeovf(filename, representation='bin8')  # write the file
        >>> os.path.isfile(filename)
        True
        >>> field_read = df.Field.fromfile(filename)  # read the file
        >>> field_read == field
        True
        >>> os.remove(filename)  # delete the file

        .. seealso:: :py:func:`~discretisedfield.Field.fromfile`

        """
        if self.dim != 1 and self.dim != 3:
            msg = (f'Cannot write dim={self.dim} field.')
            raise TypeError(msg)

        if extend_scalar and self.dim == 1:
            write_dim = 3
        else:
            write_dim = self.dim

        header = ['OOMMF OVF 2.0',
                  '',
                  'Segment count: 1',
                  '',
                  'Begin: Segment',
                  'Begin: Header',
                  '',
                  'Title: Field',
                  'Desc: File generated by Field class',
                  'meshunit: m',
                  'meshtype: rectangular',
                  f'xbase: {self.mesh.region.pmin[0] + self.mesh.cell[0]/2}',
                  f'ybase: {self.mesh.region.pmin[1] + self.mesh.cell[1]/2}',
                  f'zbase: {self.mesh.region.pmin[2] + self.mesh.cell[2]/2}',
                  f'xnodes: {self.mesh.n[0]}',
                  f'ynodes: {self.mesh.n[1]}',
                  f'znodes: {self.mesh.n[2]}',
                  f'xstepsize: {self.mesh.cell[0]}',
                  f'ystepsize: {self.mesh.cell[1]}',
                  f'zstepsize: {self.mesh.cell[2]}',
                  f'xmin: {self.mesh.region.pmin[0]}',
                  f'ymin: {self.mesh.region.pmin[1]}',
                  f'zmin: {self.mesh.region.pmin[2]}',
                  f'xmax: {self.mesh.region.pmax[0]}',
                  f'ymax: {self.mesh.region.pmax[1]}',
                  f'zmax: {self.mesh.region.pmax[2]}',
                  f'valuedim: {write_dim}',
                  'valuelabels: field_x field_y field_z',
                  'valueunits: None None None',
                  '',
                  'End: Header',
                  '']

        if representation == 'bin4':
            header.append('Begin: Data Binary 4')
            footer = ['End: Data Binary 4',
                      'End: Segment']
        elif representation == 'bin8':
            header.append('Begin: Data Binary 8')
            footer = ['End: Data Binary 8',
                      'End: Segment']
        elif representation == 'txt':
            header.append('Begin: Data Text')
            footer = ['End: Data Text',
                      'End: Segment']

        # Write header lines
        with open(filename, 'w') as f:
            f.write(''.join(map(lambda line: f'# {line}\n', header)))

        binary_reps = {'bin4': (1234567.0, 'f'),
                       'bin8': (123456789012345.0, 'd')}

        # Write data.
        if representation in binary_reps:
            # Reopen with binary write, appending to the end of the file.
            with open(filename, 'ab') as f:
                # Add the binary checksum.
                packarray = [binary_reps[representation][0]]

                # Write data to the ovf file.
                for point, value in self:
                    if self.dim == 3:
                        v = value
                    else:
                        if extend_scalar:
                            v = [value, 0.0, 0.0]
                        else:
                            v = [value]
                    for vi in v:
                        packarray.append(vi)

                format = binary_reps[representation][1]*len(packarray)
                f.write(struct.pack(format, *packarray))
        else:
            # Reopen with txt representation, appending to the end of the file.
            with open(filename, 'a') as f:
                for point, value in self:
                    if self.dim == 3:
                        v = value
                    else:
                        if extend_scalar:
                            v = [value, 0.0, 0.0]
                        else:
                            v = [value]
                    for vi in v:
                        f.write(f' {str(vi)}')
                    f.write('\n')

        # Write footer lines to OOMMF file.
        with open(filename, 'a') as f:
            f.write(''.join(map(lambda line: f'# {line}\n', footer)))

    def _writevtk(self, filename):
        """Write the field to a VTK file.

        The data is saved as a ``RECTILINEAR_GRID`` dataset. Scalar field
        (``dim=1``) is saved as ``SCALARS``. On the other hand, vector field
        (``dim=3``) is saved as both ``VECTORS`` as well as ``SCALARS`` for all
        three components to enable easy colouring of vectors in some
        visualisation packages.

        The saved VTK file can be opened with `Paraview
        <https://www.paraview.org/>`_ or `Mayavi
        <https://docs.enthought.com/mayavi/mayavi/>`_.

        Parameters
        ----------
        filename : str

            File name with an extension.

        Example
        -------
        1. Write field to a VTK file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 5e-9, 3e-9)
        >>> n = (10, 5, 3)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        >>> value_fun = lambda point: (point[0], point[1], point[2])
        >>> field = df.Field(mesh, dim=3, value=value_fun)
        ...
        >>> filename = 'mytestfile.vtk'
        >>> field._writevtk(filename)  # write the file
        >>> os.path.isfile(filename)
        True
        >>> os.remove(filename)  # delete the file

        """
        header = ['# vtk DataFile Version 3.0',
                  'Field',
                  'ASCII',
                  'DATASET RECTILINEAR_GRID',
                  'DIMENSIONS {} {} {}'.format(*self.mesh.n),
                  f'X_COORDINATES {self.mesh.n[0]} float',
                  ' '.join(map(str, self.mesh.axis_points('x'))),
                  f'Y_COORDINATES {self.mesh.n[1]} float',
                  ' '.join(map(str, self.mesh.axis_points('y'))),
                  f'Z_COORDINATES {self.mesh.n[2]} float',
                  ' '.join(map(str, self.mesh.axis_points('z'))),
                  f'POINT_DATA {len(self.mesh)}']

        if self.dim == 1:
            data = dfu.vtk_scalar_data(self, 'field')
        elif self.dim == 3:
            x, y, z = self.components
            data = dfu.vtk_scalar_data(getattr(self, x), 'x-component')
            data += dfu.vtk_scalar_data(getattr(self, y), 'y-component')
            data += dfu.vtk_scalar_data(getattr(self, z), 'z-component')
            data += dfu.vtk_vector_data(self, 'field')

        with open(filename, 'w') as f:
            f.write('\n'.join(header+data))

    def _writehdf5(self, filename):
        """Write the field to an HDF5 file.

        Parameters
        ----------
        filename : str

            Name with an extension of the file written.

        Example
        -------
        1. Write field to an HDF5 file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> p1 = (0, 0, 0)
        >>> p2 = (10e-9, 5e-9, 3e-9)
        >>> n = (10, 5, 3)
        >>> mesh = df.Mesh(p1=p1, p2=p2, n=n)
        >>> value_fun = lambda point: (point[0], point[1], point[2])
        >>> field = df.Field(mesh, dim=3, value=value_fun)
        ...
        >>> filename = 'mytestfile.h5'
        >>> field._writehdf5(filename)  # write the file
        >>> os.path.isfile(filename)
        True
        >>> field_read = df.Field.fromfile(filename)  # read the file
        >>> field_read == field
        True
        >>> os.remove(filename)  # delete the file

        .. seealso:: :py:func:`~discretisedfield.Field.fromfile`

        """
        with h5py.File(filename, 'w') as f:
            # Set up the file structure
            gfield = f.create_group('field')
            gmesh = gfield.create_group('mesh')
            gregion = gmesh.create_group('region')

            # Save everything as datasets
            gregion.create_dataset('p1', data=self.mesh.region.p1)
            gregion.create_dataset('p2', data=self.mesh.region.p2)
            gmesh.create_dataset('n', dtype='i4', data=self.mesh.n)
            gfield.create_dataset('dim', dtype='i4', data=self.dim)
            gfield.create_dataset('array', data=self.array)

    @classmethod
    def fromfile(cls, filename):
        """Read the field from an OVF (1.0 or 2.0), VTK, or HDF5 file.

        The extension of the ``filename`` should correspond to either:
            - OVF (``.ovf``, ``.omf``, ``.ohf``, ``.oef``)
            - VTK (``.vtk``), or
            - HDF5 (``.hdf5`` or ``.h5``).

        This is a ``classmethod`` and should be called as, for instance,
        ``discretisedfield.Field.fromfile('myfile.omf')``.

        Parameters
        ----------
        filename : str

            Name of the file to be read.

        Returns
        -------
        discretisedfield.Field

            Field read from the file.

        Example
        -------
        1. Read the field from an OVF file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> dirname = os.path.join(os.path.dirname(__file__),
        ...                        'tests', 'test_sample')
        >>> filename = os.path.join(dirname, 'oommf-ovf2-bin4.omf')
        >>> field = df.Field.fromfile(filename)
        >>> field
        Field(mesh=...)

        2. Read a field from the VTK file.

        >>> filename = os.path.join(dirname, 'vtk-file.vtk')
        >>> field = df.Field.fromfile(filename)
        >>> field
        Field(mesh=...)

        3. Read a field from the HDF5 file.

        >>> filename = os.path.join(dirname, 'hdf5-file.hdf5')
        >>> field = df.Field.fromfile(filename)
        >>> field
        Field(mesh=...)

        .. seealso:: :py:func:`~discretisedfield.Field._fromovf`
        .. seealso:: :py:func:`~discretisedfield.Field._fromhdf5`
        .. seealso:: :py:func:`~discretisedfield.Field._fromhdf5`
        .. seealso:: :py:func:`~discretisedfield.Field.write`

        """
        if any([filename.endswith(ext) for ext in ['.omf', '.ovf',
                                                   '.ohf', '.oef']]):
            return cls._fromovf(filename)
        elif any([filename.endswith(ext) for ext in ['.vtk']]):
            return cls._fromvtk(filename)
        elif any([filename.endswith(ext) for ext in ['.hdf5', '.h5']]):
            return cls._fromhdf5(filename)
        else:
            msg = (f'Reading file with extension {filename.split(".")[-1]} '
                   f'not supported.')
            raise ValueError(msg)

    @classmethod
    def _fromovf(cls, filename):
        """Read the field from an OVF file.

        Data representation (``txt``, ``bin4``, or ``bin8``) as well as the OVF
        version (OVF1.0 or OVF2.0) are extracted from the file itself.

        This is a ``classmethod`` and should be called as, for instance,
        ``discretisedfield.Field._fromovf('myfile.omf')``.

        Parameters
        ----------
        filename : str

            Name of the file to be read.

        Returns
        -------
        discretisedfield.Field

            Field read from the file.

        Example
        -------
        1. Read a field from the OVF file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> dirname = os.path.join(os.path.dirname(__file__),
        ...                        'tests', 'test_sample')
        >>> filename = os.path.join(dirname, 'oommf-ovf2-bin8.omf')
        >>> field = df.Field._fromovf(filename)
        >>> field
        Field(mesh=...)

        .. seealso:: :py:func:`~discretisedfield.Field._writeovf`

        """
        # valuedim is not in OVF1 metadata and has to be extracted
        # from the data itself.
        mdatalist = ['xmin', 'ymin', 'zmin',
                     'xmax', 'ymax', 'zmax',
                     'xstepsize', 'ystepsize', 'zstepsize',
                     'valuedim']
        mdatadict = dict()

        try:
            # Encoding in open is important on Windows.
            with open(filename, 'r', encoding='utf-8') as ovffile:
                lines = ovffile.readlines()

            mdatalines = list(filter(lambda s: s.startswith('#'), lines))
            datalines = np.genfromtxt(filter(lambda s: not s.startswith('#'),
                                             lines), dtype=None)

            if '1.0' in mdatalines[0]:
                # valuedim is not in OVF1 file.
                mdatadict['valuedim'] = datalines.shape[-1]

            for line in mdatalines:
                for mdatum in mdatalist:
                    if mdatum in line:
                        mdatadict[mdatum] = float(line.split()[-1])
                        break

        except UnicodeDecodeError:
            with open(filename, 'rb') as ovffile:
                f = ovffile.read()
                lines = f.split(b'\n')

            mdatalines = list(filter(lambda s: s.startswith(bytes('#',
                                                                  'utf-8')),
                                     lines))

            if bytes('2.0', 'utf-8') in mdatalines[0]:
                endian = '<'  # little-endian
            elif bytes('1.0', 'utf-8') in mdatalines[0]:
                endian = '>'  # big-endian

            for line in mdatalines:
                for mdatum in mdatalist:
                    if bytes(mdatum, 'utf-8') in line:
                        mdatadict[mdatum] = float(line.split()[-1])
                        break

            header = b'# Begin: Data Binary '
            data_start = f.find(header)
            header = f[data_start:(data_start + len(header) + 1)]

            data_start += len(header)
            data_end = f.find(b'# End: Data Binary ')

            if b'4' in header:
                nbytes = 4
                formatstr = endian + 'f'
                checkvalue = 1234567.0
            elif b'8' in header:
                nbytes = 8
                formatstr = endian + 'd'
                checkvalue = 123456789012345.0

            newlines = [b'\n\r', b'\r\n', b'\n']  # ordered by length
            for nl in newlines:
                if f.startswith(nl, data_start):
                    data_start += len(nl)
                    # There is a difference between files written by OOMMF and
                    # mumax3. OOMMF has a newline character before the "end
                    # metadata line', whereas mumax3 does not. Therefore if the
                    # length of the data stream is not a multiple of nbytes, we
                    # have to subtract the length of newline character from
                    # data_end.
                    if (data_end - data_start) % nbytes != 0:
                        data_end -= len(nl)
                    break

            listdata = list(struct.iter_unpack(formatstr,
                                               f[data_start:data_end]))
            datalines = np.array(listdata)

            if datalines[0] != checkvalue:
                # These two lines cannot be accessed via tests. Therefore, they
                # are excluded from coverage.
                msg = 'Error in checksum comparison.'  # pragma: no cover
                raise AssertionError(msg)  # pragma: no cover

            datalines = datalines[1:]

        p1 = (mdatadict[key] for key in ['xmin', 'ymin', 'zmin'])
        p2 = (mdatadict[key] for key in ['xmax', 'ymax', 'zmax'])
        cell = (mdatadict[key] for key in ['xstepsize', 'ystepsize',
                                           'zstepsize'])

        mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), cell=cell)

        # valuedim is not in OVF1 file and for binary data it has to be
        # extracted here.
        if 'valuedim' not in mdatadict.keys():
            mdatadict['valuedim'] = len(datalines) / len(mesh)

        r_tuple = (*tuple(reversed(mesh.n)), int(mdatadict['valuedim']))
        t_tuple = (*tuple(reversed(range(3))), 3)

        return cls(mesh, dim=int(mdatadict['valuedim']),
                   value=datalines.reshape(r_tuple).transpose(t_tuple))

    @classmethod
    def _fromvtk(cls, filename):
        """Read the field from a VTK file.

        This method reads the field from a VTK file defined on
        STRUCTURED_POINTS written by ``discretisedfield._writevtk``.

        This is a ``classmethod`` and should be called as, for instance,
        ``discretisedfield.Field._fromvtk('myfile.vtk')``.

        Parameters
        ----------
        filename : str

            Name of the file to be read.

        Returns
        -------
        discretisedfield.Field

            Field read from the file.

        Example
        -------
        1. Read a field from the VTK file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> dirname = os.path.join(os.path.dirname(__file__),
        ...                        'tests', 'test_sample')
        >>> filename = os.path.join(dirname, 'vtk-file.vtk')
        >>> field = df.Field._fromvtk(filename)
        >>> field
        Field(mesh=...)

        .. seealso:: :py:func:`~discretisedfield.Field._writevtk`

        """
        with open(filename, 'r') as f:
            content = f.read()
        lines = content.split('\n')

        # Determine the dimension of the field.
        if 'VECTORS' in content:
            dim = 3
            data_marker = 'VECTORS'
            skip = 0  # after how many lines data starts after marker
        else:
            dim = 1
            data_marker = 'SCALARS'
            skip = 1

        # Extract the metadata
        mdatalist = ['X_COORDINATES', 'Y_COORDINATES', 'Z_COORDINATES']
        n = []
        cell = []
        origin = []
        for i, line in enumerate(lines):
            for mdatum in mdatalist:
                if mdatum in line:
                    n.append(int(line.split()[1]))
                    coordinates = list(map(float, lines[i+1].split()))
                    origin.append(coordinates[0])
                    if len(coordinates) > 1:
                        cell.append(coordinates[1] - coordinates[0])
                    else:
                        # If only one cell exists, 1nm cell is used by default.
                        cell.append(1e-9)

        # Create objects from metadata info
        p1 = np.subtract(origin, np.multiply(cell, 0.5))
        p2 = np.add(p1, np.multiply(n, cell))
        mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), n=n)
        field = cls(mesh, dim=dim)

        # Find where data starts.
        for i, line in enumerate(lines):
            if line.startswith(data_marker):
                start_index = i
                break

        # Extract data.
        for i, line in zip(mesh.indices, lines[start_index+skip+1:]):
            if not line[0].isalpha():
                field.array[i] = list(map(float, line.split()))

        return field

    @classmethod
    def _fromhdf5(cls, filename):
        """Read the field from an HDF5 file.

        This method reads the field from an HDF5 file defined on written by
        ``discretisedfield._writevtk``.

        This is a ``classmethod`` and should be called as, for instance,
        ``discretisedfield.Field._fromhdf5('myfile.h5')``.

        Parameters
        ----------
        filename : str

            Name of the file to be read.

        Returns
        -------
        discretisedfield.Field

            Field read from the file.

        Example
        -------
        1. Read a field from the HDF5 file.

        >>> import os
        >>> import discretisedfield as df
        ...
        >>> dirname = os.path.join(os.path.dirname(__file__),
        ...                        'tests', 'test_sample')
        >>> filename = os.path.join(dirname, 'hdf5-file.hdf5')
        >>> field = df.Field._fromhdf5(filename)
        >>> field
        Field(mesh=...)

        .. seealso:: :py:func:`~discretisedfield.Field._writehdf5`

        """
        with h5py.File(filename, 'r') as f:
            # Read data from the file.
            p1 = f['field/mesh/region/p1']
            p2 = f['field/mesh/region/p2']
            n = np.array(f['field/mesh/n']).tolist()
            dim = np.array(f['field/dim']).tolist()
            array = f['field/array']

            # Create field.
            mesh = df.Mesh(region=df.Region(p1=p1, p2=p2), n=n)
            return cls(mesh, dim=dim, value=array[:])

    @property
    def mpl(self):
        """Plot interface, matplotlib based."""
        return dfp.Mpl(self)

    @property
    def k3d(self):
        """Plot interface, k3d based."""
        return dfp.K3d(self)

    @property
    def fftn(self):
        """Fourier transform.

        Computes 3D FFT for "normal" fields, 2D FFT if the field is sliced.

        Returns
        -------
        discretisedfield.Field
        """
        mesh = self._fft_mesh()

        values = []
        for idx in range(self.dim):
            ft = np.fft.fftshift(np.fft.fftn(self.array[..., idx].squeeze()))
            values.append(ft.reshape(mesh.n))

        return self.__class__(mesh, dim=len(values),
                              value=np.stack(values, axis=3),
                              components=self.components)

    @property
    def ifftn(self):
        """Inverse Fourier transform.

        Returns
        -------
        discretisedfield.Field
        """
        mesh = self.mesh.attributes['realspace_mesh']
        if self.mesh.attributes['isplane'] and not mesh.attributes['isplane']:
            mesh = mesh.plane(dfu.raxesdict[self.mesh.attributes['planeaxis']])

        values = []
        for idx in range(self.dim):
            ft = np.fft.ifftn(np.fft.ifftshift(self.array[..., idx].squeeze()))
            values.append(ft.reshape(mesh.n))

        return self.__class__(mesh, dim=len(values),
                              value=np.stack(values, axis=3),
                              components=self.components)

    @property
    def rfftn(self):
        """Real Fourier transform.

        Returns
        -------
        discretisedfield.Field
        """
        mesh = self._fft_mesh(rfft=True)

        values = []
        for idx in range(self.dim):
            array = self.array[..., idx].squeeze()
            # no shifting for the last axis
            ft = np.fft.fftshift(np.fft.rfftn(array),
                                 axes=range(len(array.shape) - 1))
            values.append(ft.reshape(mesh.n))

        return self.__class__(mesh, dim=len(values),
                              value=np.stack(values, axis=3),
                              components=self.components)

    @property
    def irfftn(self):
        """Inverse real Fourier transform.

        Returns
        -------
        discretisedfield.Field
        """
        mesh = self.mesh.attributes['realspace_mesh']
        if self.mesh.attributes['isplane'] and not mesh.attributes['isplane']:
            mesh = mesh.plane(dfu.raxesdict[self.mesh.attributes['planeaxis']])

        values = []
        for idx in range(self.dim):
            array = self.array[..., idx].squeeze()
            ft = np.fft.irfftn(
                np.fft.ifftshift(array, axes=range(len(array.shape) - 1)),
                s=[i for i in mesh.n if i > 1])
            values.append(ft.reshape(mesh.n))

        return self.__class__(mesh, dim=len(values),
                              value=np.stack(values, axis=3),
                              components=self.components)

    def _fft_mesh(self, rfft=False):
        """FFT can be one of fftfreq, rfftfreq."""
        p1 = []
        p2 = []
        n = []
        for i in range(3):
            if self.mesh.n[i] == 1:
                p1.append(0)
                p2.append(1 / self.mesh.cell[i])
                n.append(1)
            else:
                freqs = np.fft.fftshift(np.fft.fftfreq(self.mesh.n[i],
                                                       self.mesh.cell[i]))
                # Shift the region boundaries to get the correct coordinates of
                # mesh cells.
                dfreq = (freqs[1] - freqs[0]) / 2
                p1.append(min(freqs) - dfreq)
                p2.append(max(freqs) + dfreq)
                n.append(len(freqs))

        if rfft:
            # last frequency is different for rfft
            for i in [2, 1, 0]:
                if self.mesh.n[i] > 1:
                    freqs = np.fft.rfftfreq(self.mesh.n[i], self.mesh.cell[i])
                    dfreq = (freqs[1] - freqs[0]) / 2
                    p1[i] = (min(freqs) - dfreq)
                    p2[i] = (max(freqs) + dfreq)
                    n[i] = len(freqs)
                    break

        # TODO: Using PlaneMesh will simplify the code a lot here.
        mesh = df.Mesh(p1=p1, p2=p2, n=n)
        if self.mesh.attributes['isplane']:
            mesh = mesh.plane(dfu.raxesdict[self.mesh.attributes['planeaxis']])

        mesh.attributes['realspace_mesh'] = self.mesh
        mesh.attributes['fourierspace'] = True
        mesh.attributes['unit'] = rf'({mesh.attributes["unit"]})$^{{-1}}$'
        return mesh

    @property
    def real(self):
        """Real part of complex field."""
        return self.__class__(self.mesh, dim=self.dim, value=self.array.real,
                              components=self.components)

    @property
    def imag(self):
        """Imaginary part of complex field."""
        return self.__class__(self.mesh, dim=self.dim, value=self.array.imag,
                              components=self.components)

    @property
    def phase(self):
        """Phase of complex field."""
        return self.__class__(self.mesh, dim=self.dim,
                              value=np.angle(self.array),
                              components=self.components)

    @property
    def conjugate(self):
        """Complex conjugate of complex field."""
        return self.__class__(self.mesh, dim=self.dim,
                              value=self.array.conjugate(),
                              components=self.components)

    # TODO check and write tests
    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
        """Field class support for numpy ``ufuncs``."""
        # See reference implementation at:
        # https://numpy.org/doc/stable/reference/generated/numpy.lib.mixins.NDArrayOperatorsMixin.html#numpy.lib.mixins.NDArrayOperatorsMixin
        for x in inputs:
            if not isinstance(x, (Field, np.ndarray, numbers.Number)):
                return NotImplemented
        out = kwargs.get('out', ())
        if out:
            for x in out:
                if not isinstance(x, Field):
                    return NotImplemented

        mesh = [x.mesh for x in inputs if isinstance(x, Field)]
        inputs = tuple(x.array if isinstance(x, Field) else x
                       for x in inputs)
        if out:
            kwargs['out'] = tuple(x.array for x in out)

        result = getattr(ufunc, method)(*inputs, **kwargs)
        if isinstance(result, tuple):
            if len(result) != len(mesh):
                raise ValueError('wrong number of Field objects')
            return tuple(self.__class__(m, dim=x.shape[-1], value=x,
                                        components=self.components)
                         for x, m in zip(result, mesh))
        elif method == 'at':
            return None
        else:
            return self.__class__(mesh[0], dim=result.shape[-1], value=result,
                                  components=self.components)
