def add_cubic_gather_for_mode(m, Fr, Ft, Fz, exptheta_m, Fr_grid, Ft_grid,
    Fz_grid, ir_lowest, iz_lowest, Sr_arr, Sz_arr, Nr, Nz):
    """
    Add the contribution of the gathered field from azimuthal mode `m` to the
    fields felt by one macroparticle (`Fr`, `Ft`, `Fz`), using cubic weights.

    Parameters:
    -----------
    m: int
        The index of the azimuthal mode that is added.

    Fr, Ft, Fz: floats
        The fields felt by one macroparticle, which represent either E or B
        (before the contribution of mode `m` has been added)

    exptheta_m: complex
        The complex azimuthal factor $e^{-i m 	heta}$ where $	heta$ is
        the azimuthal position of the macroparticle considered.

    Fr_grid, Ft_grid, Fz_grid: 2darrays of complexs
        The fields on the interpolation grid for mode `m`

    ir_lowest, iz_lowest: ints
        The lowest indices in r and z from which the macroparticle
        considered should gather the fields (in the arrays F*_grid)
        These indices can in fact be negative and out-of-bound
        (e.g. for particles close to the axis) but get corrected within
        this function.

    Sr_arr, Sz_arr: 1darrays containing 4 floats
        The weights in r and z with which the macroparticle
        considered should gather the fields (in the arrays F*_grid)

    Nr, Nz: ints
        Dimensions of the field arrays.

    Returns:
    --------
    Fr, Ft, Fz: floats
        The fields felt by one macroparticle, which represent either E or B
        (after the contribution of mode `m` has been added)
    """
    Fr_m = 0.0j
    Ft_m = 0.0j
    Fz_m = 0.0j
    for index_r in range(4):
        ir = ir_lowest + index_r
        Sr_long = Sr_arr[index_r]
        Sr_perp = Sr_long
        if ir < 0:
            Sr_long *= (-1) ** m
            Sr_perp *= -(-1) ** m
        if ir < 0:
            ir = abs(ir) - 1
        elif ir > Nr - 1:
            ir = Nr - 1
        for index_z in range(4):
            iz = iz_lowest + index_z
            Sz = Sz_arr[index_z]
            if iz < 0:
                iz += Nz
            elif iz > Nz - 1:
                iz -= Nz
            Fr_m += Sz * Sr_perp * Fr_grid[iz, ir]
            Ft_m += Sz * Sr_perp * Ft_grid[iz, ir]
            Fz_m += Sz * Sr_long * Fz_grid[iz, ir]
    if m == 0:
        factor = 1.0
    else:
        factor = 2.0
    Fr += factor * (Fr_m * exptheta_m).real
    Ft += factor * (Ft_m * exptheta_m).real
    Fz += factor * (Fz_m * exptheta_m).real
    return Fr, Ft, Fz
