def METHOD_NAME(
    self, R, R_sersic, n_sersic, max_R_frac=1000.0, alpha=1.0, R_break=0.0
):
    """
    :param R: radius (array or float)
    :param R_sersic: Sersic radius (half-light radius)
    :param n_sersic: Sersic index (float)
    :param max_R_frac: maximum window outside which the mass is zeroed, in units of R_sersic (float)
    :return: kernel of the Sersic surface brightness at R
    """
    R_ = self._R_stable(R)
    R_sersic_ = self._R_stable(R_sersic)
    bn = self.b_n(n_sersic)
    R_frac = R_ / R_sersic_
    if isinstance(R_, int) or isinstance(R_, float):
        if R_frac > max_R_frac:
            result = 0
        else:
            exponent = -bn * (R_frac ** (1.0 / n_sersic) - 1.0)
            result = np.exp(exponent)
    else:
        R_frac_real = R_frac[R_frac <= max_R_frac]
        exponent = -bn * (R_frac_real ** (1.0 / n_sersic) - 1.0)
        result = np.zeros_like(R_)
        result[R_frac <= max_R_frac] = np.exp(exponent)
    return np.nan_to_num(result)