def time_average(rho, t):
    r"""Return a time-averaged density matrix (using trapezium rule).

    We test a basic two-level system.

    >>> import numpy as np
    >>> from scipy.constants import physical_constants
    >>> from sympy import Matrix, symbols
    >>> from fast.electric_field import electric_field_amplitude_top
    >>> from fast.symbolic import (define_laser_variables,
    ...                            polarization_vector)

    >>> Ne = 2
    >>> Nl = 1
    >>> a0 = physical_constants["Bohr radius"][0]
    >>> rm = [np.array([[0, 0], [a0, 0]]),
    ...       np.array([[0, 0], [0, 0]]),
    ...       np.array([[0, 0], [0, 0]])]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> omega_level = [0, 1.0e9]
    >>> gamma21 = 2*np.pi*6e6
    >>> gamma = np.array([[0, -gamma21], [gamma21, 0]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)

    We define symbolic variables to be used as token arguments.
    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]
    >>> detuning_knob = [symbols("delta1", real=True)]

    An map to unfold the density matrix.
    >>> unfolding = Unfolding(Ne, True, True, True)

    We obtain a function to calculate Hamiltonian terms.
    >>> aux = (Ep, epsilonps, detuning_knob, gamma,
    ...        omega_level, rm, xi, theta)
    >>> time_evolution = fast_time_evolution(*aux)

    We specify values for the variables
    >>> detuning_knobs = [100e6]
    >>> Eps = electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")
    >>> Eps *= np.exp(1j*np.pi)
    >>> Eps = [Eps]

    >>> t = np.linspace(0, 1e-6, 11)
    >>> rho0 = np.array([[1, 0], [0, 0]])
    >>> rho0 = unfolding(rho0)

    >>> rho = time_evolution(t, rho0, Eps, detuning_knobs)
    >>> print(time_average(rho, t))
    [ 0.0175  0.1244 -0.0222]

    """
    T = t[-1]-t[0]
    dt = t[1]-t[0]
    rhoav = np.sum(rho[1:-1], axis=0) + 0.5*(rho[0]+rho[-1])
    rhoav = dt/T*rhoav
    return rhoav
r"""Return a time-averaged density matrix (using trapezium rule).

    We test a basic two-level system.

    >>> import numpy as np
    >>> from scipy.constants import physical_constants
    >>> from sympy import Matrix, symbols
    >>> from fast.electric_field import electric_field_amplitude_top
    >>> from fast.symbolic import (define_laser_variables,
    ...                            polarization_vector)

    >>> Ne = 2
    >>> Nl = 1
    >>> a0 = physical_constants["Bohr radius"][0]
    >>> rm = [np.array([[0, 0], [a0, 0]]),
    ...       np.array([[0, 0], [0, 0]]),
    ...       np.array([[0, 0], [0, 0]])]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> omega_level = [0, 1.0e9]
    >>> gamma21 = 2*np.pi*6e6
    >>> gamma = np.array([[0, -gamma21], [gamma21, 0]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)

    We define symbolic variables to be used as token arguments.
    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]
    >>> detuning_knob = [symbols("delta1", real=True)]

    An map to unfold the density matrix.
    >>> unfolding = Unfolding(Ne, True, True, True)

    We obtain a function to calculate Hamiltonian terms.
    >>> aux = (Ep, epsilonps, detuning_knob, gamma,
    ...        omega_level, rm, xi, theta)
    >>> time_evolution = fast_time_evolution(*aux)

    We specify values for the variables
    >>> detuning_knobs = [100e6]
    >>> Eps = electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")
    >>> Eps *= np.exp(1j*np.pi)
    >>> Eps = [Eps]

    >>> t = np.linspace(0, 1e-6, 11)
    >>> rho0 = np.array([[1, 0], [0, 0]])
    >>> rho0 = unfolding(rho0)

    >>> rho = time_evolution(t, rho0, Eps, detuning_knobs)
    >>> print(time_average(rho, t))
    [ 0.0175  0.1244 -0.0222]