Goodwin SCN circadian network (20 oscillators)

def _light(t, i):
    """Photic drive: half-sinusoid with 24h period, only for VL neurons."""
    if i >= _N_VL:
        return 0.0
    phase = np.sin(2.0 * np.pi * t / 24.0)
    return _L0 * max(phase, 0.0)

def goodwin_scn_20_rhs(t, y):
    Y = y.reshape(_N_OSC, 3)
    X = Y[:, 0]
    F = np.mean(X)               # mean-field coupling

    K1n = _K1 ** _HILL_N
    dy = np.empty(3 * _N_OSC)

    for i in range(_N_OSC):
        xi, yi, zi = Y[i, 0], Y[i, 1], Y[i, 2]
        zi_safe = max(zi, 0.0)
        xi_safe = max(xi, 0.0)
        yi_safe = max(yi, 0.0)

        # dX/dt
        repression = K1n / (K1n + zi_safe ** _HILL_N)
        deg_x = _V2 * xi_safe / (_K2 + xi_safe)
        dx = _V1[i] * repression + _light(t, i) - deg_x + _G_COUPLING * (F - xi_safe)

        # dY/dt
        deg_y = _V4 * yi_safe / (_K4 + yi_safe)
        dyi = _K3 * xi_safe - deg_y

        # dZ/dt
        deg_z = _V6 * zi_safe / (_K6 + zi_safe)
        dzi = _K5 * yi_safe - deg_z

        base = 3 * i
        dy[base]     = dx
        dy[base + 1] = dyi
        dy[base + 2] = dzi

    return dy