def kernel(r, h):
    """
    Compute the cubic spline softened kernel

    Arguments:
      -r : Radial distance [FLOAT]
      -h : Gravitational softening length [FLOAT]

    Returns:
      Kernel value
    """

    if h == 0.0:
        return -1.0 / r

    hinv = 1.0 / h
    q = r * hinv
    if q <= 0.5:
        return (-2.8 + q * q * (5.33333333333333333 + q * q * (6.4 * q - 9.6))) * hinv
    elif q <= 1:
        return (
            -3.2
            + 0.066666666666666666666 / q
            + q
            * q
            * (
                10.666666666666666666666
                + q * (-16.0 + q * (9.6 - 2.1333333333333333333333 * q))
            )
        ) * hinv
    else:
        return -1.0 / r
    return