Python scipy.special.hankel2() Examples

def theo_fun(k):
    r"""Returns the value of Theodorsen's function at a reduced frequency :math:`k`.

    .. math:: \mathcal{C}(jk) = \frac{H_1^{(2)}(k)}{H_1^{(2)}(k) + jH_0^{(2)}(k)}

    where :math:`H_0^{(2)}(k)` and :math:`H_1^{(2)}(k)` are Hankel functions of the second kind.

    Args:
        k (np.array): Reduced frequency/frequencies at which to evaluate the function.

    Returns:
        np.array: Value of Theodorsen's function evaluated at the desired reduced frequencies.

    """

    H1 = scsp.hankel2(1, k)
    H0 = scsp.hankel2(0, k)

    C = H1 / (H1 + j * H0)

    return C 

def sphankel2(n, kr):
    """Spherical Hankel (second kind) of order n at kr

    Parameters
    ----------
    n : array_like
       Order
    kr: array_like
       Argument

    Returns
    -------
    hn2 : complex float
       Spherical Hankel function hn (second kind)
    """
    n, kr = scalar_broadcast_match(n, kr)
    hn2 = _np.full(n.shape, _np.nan, dtype=_np.complex_)
    kr_nonzero = kr != 0
    hn2[kr_nonzero] = _np.sqrt(_np.pi / 2) / _np.lib.scimath.sqrt(kr[kr_nonzero]) * hankel2(n[kr_nonzero] + 0.5,
                                                                                            kr[kr_nonzero])
    return hn2 

def hankel2(n, z):
    """Bessel function of third kind (Hankel function) of order n at kr.
    Wraps scipy.special.hankel2(n, z)

    Parameters
    ----------
    n : array_like
       Order
    z: array_like
       Argument

    Returns
    -------
    H2 : array_like
       Values of Hankel function of order n at position z
    """
    return scy.hankel2(n, z) 

def test_hankel2(self):
        hank2 = special.hankel2(1,.1)
        hankrl2 = (special.jv(1,.1) - special.yv(1,.1)*1j)
        assert_almost_equal(hank2,hankrl2,8) 

def line_dipole(omega, x0, n0, grid, *, c=None):
    r"""Line source with dipole characteristics parallel to the z-axis.

    Parameters
    ----------
    omega : float
        Frequency of source.
    x0 : (3,) array_like
        Position of source. Note: third component of x0 is ignored.
    x0 : (3,) array_like
        Normal vector of the source.
    grid : triple of array_like
        The grid that is used for the sound field calculations.
        See `sfs.util.xyz_grid()`.
    c : float, optional
        Speed of sound.

    Notes
    -----
    .. math::

        G(\x-\x_0,\w) = \frac{\i k}{4} \Hankel{2}{1}{\wc|\x-\x_0|} \cos{\phi}

    """
    k = _util.wavenumber(omega, c)
    x0 = _util.asarray_1d(x0)[:2]  # ignore z-components
    n0 = _util.asarray_1d(n0)[:2]
    grid = _util.as_xyz_components(grid)
    dx = grid[:2] - x0

    r = _np.linalg.norm(dx)
    p = 1j*k/4 * _special.hankel2(1, k * r) * _np.inner(dx, n0) / r
    return _duplicate_zdirection(p, grid) 

def _print_bessel_functions(n, k):
    print(' '.join(('bessel:', str(besselj(n, k)))))
    print(' '.join(('hankel2:', str(hankel2(n, k)))))
    print(' '.join(('spbessel:', str(spbessel(n, k)))))
    print(' '.join(('dspbessel:', str(dspbessel(n, k)))))
    print(' '.join(('spneumann:', str(spneumann(n, k)))))
    print(' '.join(('dspneumann:', str(dspneumann(n, k)))))
    print(' '.join(('sphankel2:', str(sphankel2(n, k)))))
    print(' '.join(('dsphankel2:', str(dsphankel2(n, k))))) 

def test_hankel2(self):
        assert_mpmath_equal(sc.hankel2,
                            exception_to_nan(lambda v, x: mpmath.hankel2(v, x, **HYPERKW)),
                            [Arg(-1e20, 1e20), Arg()]) 

def circular_ls(N, k, r, rs, setup):
    r"""Radial coefficients for a line source.

    Computes the radial component of the circular harmonics expansion of a
    line source impinging on a circular array.

    .. math::

        \mathring{P}_n(k) = \frac{-i}{4} H_n^{(2)}(k r_s) b_n(kr)

    Parameters
    ----------
    N : int
        Maximum order.
    k : (M,) array_like
        Wavenumber.
    r : float
        Radius of microphone array.
    rs : float
        Distance of source.
    setup : {'open', 'card', 'rigid'}
        Array configuration (open, cardioids, rigid).

    Returns
    -------
    bn : (M, 2*N+1) numpy.ndarray
        Radial weights for all orders up to N and the given wavenumbers.
    """
    k = util.asarray_1d(k)
    krs = k*rs
    n = np.roll(np.arange(-N, N+1), -N)

    bn = circ_radial_weights(N, k*r, setup)
    if len(k) == 1:
        bn = bn[np.newaxis, :]
    for i, x in enumerate(krs):
        Hn = special.hankel2(n, x)
        bn[i, :] = bn[i, :] * Hn
    return -1j/4 * np.squeeze(bn) 

def test_negv2(self):
        assert_almost_equal(special.hankel2(-3,2), -special.hankel2(3,2), 14) 

def test_hankel2(self):
        cephes.hankel2(1,1) 

def test_hankel2(self):
        assert_mpmath_equal(sc.hankel2,
                            _exception_to_nan(lambda v, x: mpmath.hankel2(v, x, **HYPERKW)),
                            [Arg(-1e20, 1e20), Arg()]) 

def test_hankel2(self):
        hank2 = special.hankel2(1,.1)
        hankrl2 = (special.jv(1,.1) - special.yv(1,.1)*1j)
        assert_almost_equal(hank2,hankrl2,8) 

def test_negv2(self):
        assert_almost_equal(special.hankel2(-3,2), -special.hankel2(3,2), 14) 

def test_hankel2(self):
        cephes.hankel2(1,1) 

def test_ticket_853(self):
        """Negative-order Bessels"""
        # cephes
        assert_allclose(special.jv(-1, 1), -0.4400505857449335)
        assert_allclose(special.jv(-2, 1), 0.1149034849319005)
        assert_allclose(special.yv(-1, 1), 0.7812128213002887)
        assert_allclose(special.yv(-2, 1), -1.650682606816255)
        assert_allclose(special.iv(-1, 1), 0.5651591039924851)
        assert_allclose(special.iv(-2, 1), 0.1357476697670383)
        assert_allclose(special.kv(-1, 1), 0.6019072301972347)
        assert_allclose(special.kv(-2, 1), 1.624838898635178)
        assert_allclose(special.jv(-0.5, 1), 0.43109886801837607952)
        assert_allclose(special.yv(-0.5, 1), 0.6713967071418031)
        assert_allclose(special.iv(-0.5, 1), 1.231200214592967)
        assert_allclose(special.kv(-0.5, 1), 0.4610685044478945)
        # amos
        assert_allclose(special.jv(-1, 1+0j), -0.4400505857449335)
        assert_allclose(special.jv(-2, 1+0j), 0.1149034849319005)
        assert_allclose(special.yv(-1, 1+0j), 0.7812128213002887)
        assert_allclose(special.yv(-2, 1+0j), -1.650682606816255)

        assert_allclose(special.iv(-1, 1+0j), 0.5651591039924851)
        assert_allclose(special.iv(-2, 1+0j), 0.1357476697670383)
        assert_allclose(special.kv(-1, 1+0j), 0.6019072301972347)
        assert_allclose(special.kv(-2, 1+0j), 1.624838898635178)

        assert_allclose(special.jv(-0.5, 1+0j), 0.43109886801837607952)
        assert_allclose(special.jv(-0.5, 1+1j), 0.2628946385649065-0.827050182040562j)
        assert_allclose(special.yv(-0.5, 1+0j), 0.6713967071418031)
        assert_allclose(special.yv(-0.5, 1+1j), 0.967901282890131+0.0602046062142816j)

        assert_allclose(special.iv(-0.5, 1+0j), 1.231200214592967)
        assert_allclose(special.iv(-0.5, 1+1j), 0.77070737376928+0.39891821043561j)
        assert_allclose(special.kv(-0.5, 1+0j), 0.4610685044478945)
        assert_allclose(special.kv(-0.5, 1+1j), 0.06868578341999-0.38157825981268j)

        assert_allclose(special.jve(-0.5,1+0.3j), special.jv(-0.5, 1+0.3j)*exp(-0.3))
        assert_allclose(special.yve(-0.5,1+0.3j), special.yv(-0.5, 1+0.3j)*exp(-0.3))
        assert_allclose(special.ive(-0.5,0.3+1j), special.iv(-0.5, 0.3+1j)*exp(-0.3))
        assert_allclose(special.kve(-0.5,0.3+1j), special.kv(-0.5, 0.3+1j)*exp(0.3+1j))

        assert_allclose(special.hankel1(-0.5, 1+1j), special.jv(-0.5, 1+1j) + 1j*special.yv(-0.5,1+1j))
        assert_allclose(special.hankel2(-0.5, 1+1j), special.jv(-0.5, 1+1j) - 1j*special.yv(-0.5,1+1j)) 

def line_2d(omega, x0, n0, xs, *, c=None):
    r"""Driving function for 2-dimensional SDM for a virtual line source.

    Parameters
    ----------
    omega : float
        Angular frequency of line source.
    x0 : (N, 3) array_like
        Sequence of secondary source positions.
    n0 : (N, 3) array_like
        Sequence of normal vectors of secondary sources.
    xs : (3,) array_like
        Position of line source.
    c : float, optional
        Speed of sound.

    Returns
    -------
    d : (N,) numpy.ndarray
        Complex weights of secondary sources.
    selection : (N,) numpy.ndarray
        Boolean array containing ``True`` or ``False`` depending on
        whether the corresponding secondary source is "active" or not.
    secondary_source_function : callable
        A function that can be used to create the sound field of a
        single secondary source.  See `sfs.fd.synthesize()`.

    Notes
    -----
    The secondary sources have to be located on the x-axis (y0=0).
    Derived from :cite:`Spors2009`, Eq.(9), Eq.(4).

    Examples
    --------
    .. plot::
        :context: close-figs

        d, selection, secondary_source = sfs.fd.sdm.line_2d(
            omega, array.x, array.n, xs)
        plot(d, selection, secondary_source)

    """
    x0 = _util.asarray_of_rows(x0)
    n0 = _util.asarray_of_rows(n0)
    xs = _util.asarray_1d(xs)
    k = _util.wavenumber(omega, c)
    ds = x0 - xs
    r = _np.linalg.norm(ds, axis=1)
    d = - 1j/2 * k * xs[1] / r * _hankel2(1, k * r)
    selection = _util.source_selection_all(len(x0))
    return d, selection, _secondary_source_line(omega, c) 

def line_2d(omega, x0, n0, xs, *, c=None):
    r"""Driving function for 2-dimensional WFS for a virtual line source.

    Parameters
    ----------
    omega : float
        Angular frequency of line source.
    x0 : (N, 3) array_like
        Sequence of secondary source positions.
    n0 : (N, 3) array_like
        Sequence of normal vectors of secondary sources.
    xs : (3,) array_like
        Position of virtual line source.
    c : float, optional
        Speed of sound.

    Returns
    -------
    d : (N,) numpy.ndarray
        Complex weights of secondary sources.
    selection : (N,) numpy.ndarray
        Boolean array containing ``True`` or ``False`` depending on
        whether the corresponding secondary source is "active" or not.
    secondary_source_function : callable
        A function that can be used to create the sound field of a
        single secondary source.  See `sfs.fd.synthesize()`.

    Notes
    -----
    .. math::

        D(\x_0,\w) = \frac{\i}{2} \wc
            \frac{\scalarprod{\x-\x_0}{\n_0}}{|\x-\x_\text{s}|}
            \Hankel{2}{1}{\wc|\x-\x_\text{s}|}

    Examples
    --------
    .. plot::
        :context: close-figs

        d, selection, secondary_source = sfs.fd.wfs.line_2d(
            omega, array.x, array.n, xs)
        plot(d, selection, secondary_source)

    """
    x0 = _util.asarray_of_rows(x0)
    n0 = _util.asarray_of_rows(n0)
    xs = _util.asarray_1d(xs)
    k = _util.wavenumber(omega, c)
    ds = x0 - xs
    r = _np.linalg.norm(ds, axis=1)
    d = -1j/2 * k * _inner1d(ds, n0) / r * _hankel2(1, k * r)
    selection = _util.source_selection_line(n0, x0, xs)
    return d, selection, _secondary_source_line(omega, c) 

def plane_25d(omega, x0, n0, n=[0, 1, 0], *, xref=[0, 0, 0], c=None):
    r"""Driving function for 2.5-dimensional SDM for a virtual plane wave.

    Parameters
    ----------
    omega : float
        Angular frequency of plane wave.
    x0 : (N, 3) array_like
        Sequence of secondary source positions.
    n0 : (N, 3) array_like
        Sequence of normal vectors of secondary sources.
    n: (3,) array_like, optional
        Normal vector (traveling direction) of plane wave.
    xref : (3,) array_like, optional
        Reference point for synthesized sound field.
    c : float, optional
        Speed of sound.

    Returns
    -------
    d : (N,) numpy.ndarray
        Complex weights of secondary sources.
    selection : (N,) numpy.ndarray
        Boolean array containing ``True`` or ``False`` depending on
        whether the corresponding secondary source is "active" or not.
    secondary_source_function : callable
        A function that can be used to create the sound field of a
        single secondary source.  See `sfs.fd.synthesize()`.

    Notes
    -----
    The secondary sources have to be located on the x-axis (y0=0).
    Eq.(3.79) from :cite:`Ahrens2012`.

    Examples
    --------
    .. plot::
        :context: close-figs

        d, selection, secondary_source = sfs.fd.sdm.plane_25d(
            omega, array.x, array.n, npw, xref=[0, -1, 0])
        plot(d, selection, secondary_source)

    """
    x0 = _util.asarray_of_rows(x0)
    n0 = _util.asarray_of_rows(n0)
    n = _util.normalize_vector(n)
    xref = _util.asarray_1d(xref)
    k = _util.wavenumber(omega, c)
    d = 4j * _np.exp(-1j*k*n[1]*xref[1]) / _hankel2(0, k*n[1]*xref[1]) * \
        _np.exp(-1j*k*n[0]*x0[:, 0])
    selection = _util.source_selection_all(len(x0))
    return d, selection, _secondary_source_point(omega, c) 

def line_velocity(omega, x0, grid, *, c=None, rho0=None):
    """Velocity of line source parallel to the z-axis.

    Parameters
    ----------
    omega : float
        Frequency of source.
    x0 : (3,) array_like
        Position of source. Note: third component of x0 is ignored.
    grid : triple of array_like
        The grid that is used for the sound field calculations.
        See `sfs.util.xyz_grid()`.
    c : float, optional
        Speed of sound.

    Returns
    -------
    `XyzComponents`
        Particle velocity at positions given by *grid*.

    Examples
    --------
    The particle velocity can be plotted on top of the sound pressure:

    .. plot::
        :context: close-figs

        v = sfs.fd.source.line_velocity(omega, x0, vgrid)
        sfs.plot2d.amplitude(p * normalization_line, grid)
        sfs.plot2d.vectors(v * normalization_line, vgrid)
        plt.title("Sound Pressure and Particle Velocity")

    """
    if c is None:
        c = _default.c
    if rho0 is None:
        rho0 = _default.rho0
    k = _util.wavenumber(omega, c)
    x0 = _util.asarray_1d(x0)[:2]  # ignore z-component
    grid = _util.as_xyz_components(grid)

    offset = grid[:2] - x0
    r = _np.linalg.norm(offset)
    v = -1/(4 * c * rho0) * _special.hankel2(1, k * r)
    v = [v * o / r for o in offset]

    assert v[0].shape == v[1].shape

    if len(grid) > 2:
        v.append(_np.zeros_like(v[0]))

    return _util.XyzComponents([_duplicate_zdirection(vi, grid) for vi in v]) 

def test_ticket_853(self):
        """Negative-order Bessels"""
        # cephes
        assert_tol_equal(special.jv(-1, 1), -0.4400505857449335)
        assert_tol_equal(special.jv(-2, 1), 0.1149034849319005)
        assert_tol_equal(special.yv(-1, 1), 0.7812128213002887)
        assert_tol_equal(special.yv(-2, 1), -1.650682606816255)
        assert_tol_equal(special.iv(-1, 1), 0.5651591039924851)
        assert_tol_equal(special.iv(-2, 1), 0.1357476697670383)
        assert_tol_equal(special.kv(-1, 1), 0.6019072301972347)
        assert_tol_equal(special.kv(-2, 1), 1.624838898635178)
        assert_tol_equal(special.jv(-0.5, 1), 0.43109886801837607952)
        assert_tol_equal(special.yv(-0.5, 1), 0.6713967071418031)
        assert_tol_equal(special.iv(-0.5, 1), 1.231200214592967)
        assert_tol_equal(special.kv(-0.5, 1), 0.4610685044478945)
        # amos
        assert_tol_equal(special.jv(-1, 1+0j), -0.4400505857449335)
        assert_tol_equal(special.jv(-2, 1+0j), 0.1149034849319005)
        assert_tol_equal(special.yv(-1, 1+0j), 0.7812128213002887)
        assert_tol_equal(special.yv(-2, 1+0j), -1.650682606816255)

        assert_tol_equal(special.iv(-1, 1+0j), 0.5651591039924851)
        assert_tol_equal(special.iv(-2, 1+0j), 0.1357476697670383)
        assert_tol_equal(special.kv(-1, 1+0j), 0.6019072301972347)
        assert_tol_equal(special.kv(-2, 1+0j), 1.624838898635178)

        assert_tol_equal(special.jv(-0.5, 1+0j), 0.43109886801837607952)
        assert_tol_equal(special.jv(-0.5, 1+1j), 0.2628946385649065-0.827050182040562j)
        assert_tol_equal(special.yv(-0.5, 1+0j), 0.6713967071418031)
        assert_tol_equal(special.yv(-0.5, 1+1j), 0.967901282890131+0.0602046062142816j)

        assert_tol_equal(special.iv(-0.5, 1+0j), 1.231200214592967)
        assert_tol_equal(special.iv(-0.5, 1+1j), 0.77070737376928+0.39891821043561j)
        assert_tol_equal(special.kv(-0.5, 1+0j), 0.4610685044478945)
        assert_tol_equal(special.kv(-0.5, 1+1j), 0.06868578341999-0.38157825981268j)

        assert_tol_equal(special.jve(-0.5,1+0.3j), special.jv(-0.5, 1+0.3j)*exp(-0.3))
        assert_tol_equal(special.yve(-0.5,1+0.3j), special.yv(-0.5, 1+0.3j)*exp(-0.3))
        assert_tol_equal(special.ive(-0.5,0.3+1j), special.iv(-0.5, 0.3+1j)*exp(-0.3))
        assert_tol_equal(special.kve(-0.5,0.3+1j), special.kv(-0.5, 0.3+1j)*exp(0.3+1j))

        assert_tol_equal(special.hankel1(-0.5, 1+1j), special.jv(-0.5, 1+1j) + 1j*special.yv(-0.5,1+1j))
        assert_tol_equal(special.hankel2(-0.5, 1+1j), special.jv(-0.5, 1+1j) - 1j*special.yv(-0.5,1+1j)) 

def circ_radial_weights(N, kr, setup):
    r"""Radial weighing functions.

    Computes the radial weighting functions for diferent array types

    For instance for an rigid array

    .. math::

        b_n(kr) = J_n(kr) - \frac{J_n^\prime(kr)}{H_n^{(2)\prime}(kr)}H_n^{(2)}(kr)

    Parameters
    ----------
    N : int
        Maximum order.
    kr : (M,) array_like
        Wavenumber * radius.
    setup : {'open', 'card', 'rigid'}
        Array configuration (open, cardioids, rigid).

    Returns
    -------
    bn : (M, 2*N+1) numpy.ndarray
        Radial weights for all orders up to N and the given wavenumbers.

    """
    kr = util.asarray_1d(kr)
    n = np.arange(N+1)
    Bns = np.zeros((len(kr), N+1), dtype=complex)
    for i, x in enumerate(kr):
        Jn = special.jv(n, x)
        if setup == 'open':
            bn = Jn
        elif setup == 'card':
            bn = Jn - 1j * special.jvp(n, x, n=1)
        elif setup == 'rigid':
            if x == 0:
                # Hn(x)/Hn'(x) -> 0 for x -> 0
                bn = Jn
            else:
                Jnd = special.jvp(n, x, n=1)
                Hn = special.hankel2(n, x)
                Hnd = special.h2vp(n, x)
                bn = Jn - Jnd/Hnd*Hn
        else:
            raise ValueError('setup must be either: open, card or rigid')
        Bns[i, :] = bn
    Bns = np.concatenate((Bns, (Bns*(-1)**np.arange(N+1))[:, :0:-1]), axis=-1)
    return np.squeeze(Bns)