def rapsd(
    field, fft_method=None, return_freq=False, d=1.0, normalize=False, **fft_kwargs
):
    """
    Compute radially averaged power spectral density (RAPSD) from the given
    2D input field.

    Parameters
    ----------
    field: array_like
        A 2d array of shape (m, n) containing the input field.
    fft_method: object
        A module or object implementing the same methods as numpy.fft and
        scipy.fftpack. If set to None, field is assumed to represent the
        shifted discrete Fourier transform of the input field, where the
        origin is at the center of the array
        (see numpy.fft.fftshift or scipy.fftpack.fftshift).
    return_freq: bool
        Whether to also return the Fourier frequencies.
    d: scalar
        Sample spacing (inverse of the sampling rate). Defaults to 1.
        Applicable if return_freq is 'True'.
    normalize: bool
        If True, normalize the power spectrum so that it sums to one.

    Returns
    -------
    out: ndarray
      One-dimensional array containing the RAPSD. The length of the array is
      int(l/2) (if l is even) or int(l/2)+1 (if l is odd), where l=max(m,n).
    freq: ndarray
      One-dimensional array containing the Fourier frequencies.

    References
    ----------
    :cite:`RC2011`
    """

    if len(field.shape) != 2:
        raise ValueError(
            f"{len(field.shape)} dimensions are found, but the number "
            "of dimensions should be 2"
        )

    if np.sum(np.isnan(field)) > 0:
        raise ValueError("input field should not contain nans")

    m, n = field.shape

    yc, xc = arrays.compute_centred_coord_array(m, n)
    r_grid = np.sqrt(xc * xc + yc * yc).round()
    l = max(field.shape[0], field.shape[1])

    if l % 2 == 1:
        r_range = np.arange(0, int(l / 2) + 1)
    else:
        r_range = np.arange(0, int(l / 2))

    if fft_method is not None:
        psd = fft_method.fftshift(fft_method.fft2(field, **fft_kwargs))
        psd = np.abs(psd) ** 2 / psd.size
    else:
        psd = field

    result = []
    for r in r_range:
        mask = r_grid == r
        psd_vals = psd[mask]
        result.append(np.mean(psd_vals))

    result = np.array(result)

    if normalize:
        result /= np.sum(result)

    if return_freq:
        freq = np.fft.fftfreq(l, d=d)
        freq = freq[r_range]
        return result, freq
    else:
        return result
Compute radially averaged power spectral density (RAPSD) from the given
2D input field.

Parameters
----------
field: array_like
A 2d array of shape (m, n) containing the input field.
fft_method: object
A module or object implementing the same methods as numpy.fft and
scipy.fftpack. If set to None, field is assumed to represent the
shifted discrete Fourier transform of the input field, where the
origin is at the center of the array
(see numpy.fft.fftshift or scipy.fftpack.fftshift).
return_freq: bool
Whether to also return the Fourier frequencies.
d: scalar
Sample spacing (inverse of the sampling rate). Defaults to 1.
Applicable if return_freq is 'True'.
normalize: bool
If True, normalize the power spectrum so that it sums to one.

Returns
-------
out: ndarray
One-dimensional array containing the RAPSD. The length of the array is
int(l/2) (if l is even) or int(l/2)+1 (if l is odd), where l=max(m,n).
freq: ndarray
One-dimensional array containing the Fourier frequencies.

References
----------
:cite:`RC2011`
pysteps/utils/spectral.py