# 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


### Explanation

$$\dot{z_i} = z_i \Bigg( a f_i + b_1 f_i |z_i|^2 + b_2 f_i \frac{\varepsilon |z_i|^4}{1-\varepsilon |z_i|^2} \Bigg) + x$$

where

$$a = \alpha + j2\pi \\ b_1 = \beta_1 + j\delta_1 \\ b_2 = \beta_2 + j\delta_2$$