Oscillation

def zdot(x, z, f, zp):
    """Dynamics of a neural oscillator.

    Implements the dynamical system described in equation 20 of the paper
    referenced above.

    .. math::

        \\tau_i\\dot{z} = z \\Bigg(\\alpha + j\\omega + b_1 |z|^2 +
            \\frac{b_2\\varepsilon |z|^4}{1-\\varepsilon |z|^2} \\Bigg) +
            \\frac{x}{1-\\sqrt{\\varepsilon} x} \\frac{1}{1-\\sqrt{\\varepsilon}
            \\bar{z}}

    where

    .. math::

        b_1 &= \\beta_1 + j \\delta_1 \\\\
        b_2 &= \\beta_2 + j \\delta_2 \\\\

    Args:
        x (:class:`numpy.array`): input signal
        z (:class:`numpy.array`): oscillator state
        f (:class:`numpy.array`): natural frequency of the oscillator
        zp (:class:`.Zparam`): oscillator parameters: :math:`\\alpha,
            \\beta_1, \\delta_1, \\beta_2, \\delta_2` and :math:`\\varepsilon`

    Returns:
         (:class:`numpy.array`): The evaluated time derivative (:math:`\\dot{z}`)

    Note:
        Can work with vectors (i.e. simultaneously compute multiple oscillators).
        In that case, ``x``, ``z``, and ``f`` must have the same shape.

    """

    lin = (zp.a) * f
    nl1 = zp.b1 * f * np.abs(z) ** 2
    nl2 = zp.b2 * f * zp.e * (np.abs(z) ** 4) / (1 - zp.e * np.abs(z) ** 2)
    return z * (lin + nl1 + nl2) + x