Compute the projection of a disc using the Fourier slice theorem
    and the Bessel function of the first kind of order 1.

    Parameters
    ----------
    radius: float
        Radius of the sphere [m]
    sphere_index: float
        Refractive index of the sphere
    medium_index: float
        Refractive index of the surrounding medium
    wavelength: float
        Vacuum wavelength of the imaging light [m]
    pixel_size: float
        Pixel size [m]
    grid_size: tuple of floats
        Resulting image size in x and y [px]
    center: tuple of floats
        Center position in image coordinates [px]
    lD: float
        The axial distance [m] from the center of the sphere
        at which the field is computed.
    approx: str
        Which approximation to use (either "born" or "rytov")
    zeropad: int
        Zero-padding factor
    oversample: int
        Oversampling factor

def sphere_prop_fslice_bessel(radius, sphere_index, medium_index,
                              wavelength=550e-9, pixel_size=1e-7,
                              grid_size=(80, 80), lD=0, approx="rytov",
                              zeropad=5, oversample=1
                              ):
    """Compute the projection of a disc using the Fourier slice theorem
    and the Bessel function of the first kind of order 1.

    Parameters
    ----------
    radius: float
        Radius of the sphere [m]
    sphere_index: float
        Refractive index of the sphere
    medium_index: float
        Refractive index of the surrounding medium
    wavelength: float
        Vacuum wavelength of the imaging light [m]
    pixel_size: float
        Pixel size [m]
    grid_size: tuple of floats
        Resulting image size in x and y [px]
    center: tuple of floats
        Center position in image coordinates [px]
    lD: float
        The axial distance [m] from the center of the sphere
        at which the field is computed.
    approx: str
        Which approximation to use (either "born" or "rytov")
    zeropad: int
        Zero-padding factor
    oversample: int
        Oversampling factor
    """
    assert approx in ["born", "rytov"]
    assert oversample > 0
    assert zeropad > 0
    assert isinstance(oversample, int)
    assert isinstance(zeropad, int)

    # convert everything to pixels
    radius /= pixel_size
    wavelength /= pixel_size
    lD /= pixel_size

    # apply over-sampling and zero-padding
    wavelength *= oversample
    opad_size = np.array(grid_size) * zeropad * oversample

    assert (int(s) == s for s in opad_size), "grid_size must be integer type"
    opad_size = np.array(np.round(opad_size), dtype=int)
    grid_size = np.array(np.round(grid_size), dtype=int)

    kx = 2 * np.pi * \
        np.fft.ifftshift(np.fft.fftfreq(opad_size[0])).reshape(-1, 1)
    ky = 2 * np.pi * \
        np.fft.ifftshift(np.fft.fftfreq(opad_size[1])).reshape(1, -1)
    km = 2 * np.pi * medium_index / wavelength

    filter_klp = (kx**2 + ky**2 < km**2)

    kz = np.sqrt((km**2 - kx**2 - ky**2) * filter_klp) - km

    r = np.sqrt((kx**2 + ky**2 + kz**2) * filter_klp) / (2 * np.pi)

    comp_id = r != 0
    F = np.zeros_like(r)
    F[comp_id] = spspec.spherical_jn(1, r[comp_id] * radius * np.pi * 2) \
        * radius**2 / r[comp_id] * 2
    # center has analytical value
    center_fft = np.where(np.abs(kx) + np.abs(ky) + np.abs(kz) == 0)
    F[center_fft] = 4 / 3 * np.pi * radius**3

    # object amplitude
    F *= km**2 * ((sphere_index / medium_index)**2 - 1)

    # prefactor A
    M = 1. / km * np.sqrt((km**2 - kx**2 - ky**2) * filter_klp)

    # division factor
    A = -2j * km * M * np.exp(-1j * km * M * lD)

    # rotate phase by half a pixel so the ifft is centered in real space
    if grid_size[0] % 2:
        doffx = 0
    else:
        doffx = .5
    if grid_size[1] % 2:
        doffy = 0
    else:
        doffy = .5
    transl = np.exp(1j * ((doffx) * kx + (doffy) * ky))

    valid = F != 0
    Fconv = np.zeros((opad_size[0], opad_size[1]), dtype=complex)
    Fconv[valid] = F[valid] / A[valid] * transl[valid]

    p = np.fft.ifftn(np.fft.fftshift(Fconv))

    p = np.fft.ifftshift(p)

    if oversample > 1:
        p = p[::oversample, ::oversample]

    if zeropad > 1:
        # Slice
        a0, a1 = np.array(np.floor(opad_size / 2), dtype=int) // oversample
        b0, b1 = np.array(np.floor(grid_size / 2), dtype=int)
        if grid_size[0] % 2 != 0:
            of0 = 1
            a0 += 1
        else:
            of0 = 0
        if grid_size[1] % 2 != 0:
            of1 = 1
            a1 += 1
        else:
            of1 = 0
        # remove zero-padding
        p = p[a0 - b0:a0 + b0 + of0, a1 - b1:a1 + b1 + of1]

    if approx == "born":
        # norm = (u0 + ub)/u0
        # norm = 1 + ub/u0
        return 1 + p / np.exp(1j * km * lD)
    elif approx == "rytov":
        # norm = (u0 + ur)/u0
        # ur = u0 ( exp(ub/u0) -1 )
        # norm = ( u0 + u0 *(exp(ub/u0)-1) )/u0
        # norm = exp(ub/u0)
        return np.exp(p / np.exp(1j * km * lD))