Python scipy.special.spherical_jn() Examples

def spherical_hn2(n, z):
    r"""Spherical Hankel function of 2nd kind.

    Defined as https://dlmf.nist.gov/10.47.E6,

    .. math::

        \hankel{2}{n}{z} = \sqrt{\frac{\pi}{2z}}
        \Hankel{2}{n + \frac{1}{2}}{z},

    where :math:`\Hankel{2}{n}{\cdot}` is the Hankel function of the
    second kind and n-th order, and :math:`z` its complex argument.

    Parameters
    ----------
    n : array_like
        Order of the spherical Hankel function (n >= 0).
    z : array_like
        Argument of the spherical Hankel function.

    """
    return spherical_jn(n, z) - 1j * spherical_yn(n, z) 

def sphere_attenuation(self, q, tau, diameter):
        """Implements the finite time Callaghan model for planes."""
        radius = diameter / 2.0
        q_argument = 2 * np.pi * q * radius
        q_argument_2 = q_argument ** 2
        res = np.zeros_like(q)

        # J = special.spherical_jn(q_argument)
        Jder = special.spherical_jn(q_argument, derivative=True)
        for k in range(0, self.alpha.shape[0]):
            for n in range(0, self.alpha.shape[1]):
                a_nk2 = self.alpha[k, n] ** 2
                update = np.exp(-a_nk2 * self.Dintra * tau / radius ** 2)
                update *= ((2 * n + 1) * a_nk2) / \
                    (a_nk2 - (n - 0.5) ** 2 + 0.25)
                update *= q_argument * Jder
                update /= (q_argument_2 - a_nk2) ** 2
                res += update
        return res 

def simple_Qext(self, m, x):
        max_n = self.get_num_iters(m, x)
        all_psi, all_psi_derivs = riccati_jn(max_n, x)

        jn = spherical_jn(range(max_n + 1), m*x)
        all_mx_vals = m * x * jn
        all_mx_derivs = jn + m * x * spherical_jn(range(max_n + 1), m * x, derivative=True)
        all_D = all_mx_derivs/all_mx_vals

        all_xi = all_psi - 1j * riccati_yn(max_n, x)[0]

        all_n = np.arange(1, max_n+1)

        all_a = ((all_D[1:]/m + all_n/x)*all_psi[1:] - all_psi[0:-1])/((all_D[1:]/m + all_n/x)*all_xi[1:] - all_xi[0:-1])
        all_b = ((m*all_D[1:] + all_n/x)*all_psi[1:] - all_psi[0:-1])/((m*all_D[1:] + all_n/x)*all_xi[1:] - all_xi[0:-1])

        all_terms = 2.0/x**2 * (2*all_n + 1) * (all_a + all_b).real
        Qext = np.sum(all_terms[~np.isnan(all_terms)])
        return Qext 

def test_sph_jn(self):
        s1 = np.empty((2,3))
        x = 0.2

        s1[0][0] = spherical_jn(0, x)
        s1[0][1] = spherical_jn(1, x)
        s1[0][2] = spherical_jn(2, x)
        s1[1][0] = spherical_jn(0, x, derivative=True)
        s1[1][1] = spherical_jn(1, x, derivative=True)
        s1[1][2] = spherical_jn(2, x, derivative=True)

        s10 = -s1[0][1]
        s11 = s1[0][0]-2.0/0.2*s1[0][1]
        s12 = s1[0][1]-3.0/0.2*s1[0][2]
        assert_array_almost_equal(s1[0],[0.99334665397530607731,
                                      0.066400380670322230863,
                                      0.0026590560795273856680],12)
        assert_array_almost_equal(s1[1],[s10,s11,s12],12) 

def test_spherical_jn_yn_cross_product_1(self):
        # http://dlmf.nist.gov/10.50.E3
        n = np.array([1, 5, 8])
        x = np.array([0.1, 1, 10])
        left = (spherical_jn(n + 1, x) * spherical_yn(n, x) -
                spherical_jn(n, x) * spherical_yn(n + 1, x))
        right = 1/x**2
        assert_allclose(left, right) 

def spbessel(n, kr):
    """Spherical Bessel function (first kind) of order n at kr

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

    Returns
    -------
    J : complex float
       Spherical Bessel
    """
    n, kr = scalar_broadcast_match(n, kr)

    if _np.any(n < 0) | _np.any(kr < 0) | _np.any(_np.mod(n, 1) != 0):
        J = _np.zeros(kr.shape, dtype=_np.complex_)

        kr_non_zero = kr != 0
        J[kr_non_zero] = _np.lib.scimath.sqrt(_np.pi / 2 / kr[kr_non_zero]) * besselj(n[kr_non_zero] + 0.5,
                                                                                      kr[kr_non_zero])
        J[_np.logical_and(kr == 0, n == 0)] = 1
    else:
        J = scy.spherical_jn(n.astype(_np.int), kr)
    return _np.squeeze(J) 

def test_spherical_jn_complex(self):
        def mp_spherical_jn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besselj(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_jn(int(n.real), z),
                            exception_to_nan(mp_spherical_jn),
                            [IntArg(0, 200), ComplexArg()]) 

def test_spherical_jn(self):
        def mp_spherical_jn(n, z):
            arg = mpmath.mpmathify(z)
            out = (mpmath.besselj(n + mpmath.mpf(1)/2, arg) /
                   mpmath.sqrt(2*arg/mpmath.pi))
            if arg.imag == 0:
                return out.real
            else:
                return out

        assert_mpmath_equal(lambda n, z: sc.spherical_jn(int(n), z),
                            exception_to_nan(mp_spherical_jn),
                            [IntArg(0, 200), Arg(-1e8, 1e8)],
                            dps=300) 

def test_spherical_jn_d_zero(self):
        n = np.array([0, 1, 2, 3, 7, 15])
        assert_allclose(spherical_jn(n, 0, derivative=True),
                        np.array([0, 1/3, 0, 0, 0, 0])) 

def df(self, n, z):
        return spherical_jn(n, z, derivative=True) 

def test_spherical_jn_yn_cross_product_2(self):
        # http://dlmf.nist.gov/10.50.E3
        n = np.array([1, 5, 8])
        x = np.array([0.1, 1, 10])
        left = (spherical_jn(n + 2, x) * spherical_yn(n, x) -
                spherical_jn(n, x) * spherical_yn(n + 2, x))
        right = (2*n + 3)/x**3
        assert_allclose(left, right) 

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

    Computes the radial component of the spherical harmonics expansion of a
    point source impinging on a spherical array.

    .. math::

        \mathring{P}_n(k) = 4 \pi (-i) k 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, 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.arange(N+1)

    bn = weights(N, k*r, setup)
    if len(k) == 1:
        bn = bn[np.newaxis, :]

    for i, x in enumerate(krs):
        hn = special.spherical_jn(n, x) - 1j * special.spherical_yn(n, x)
        bn[i, :] = bn[i, :] * 4*np.pi * (-1j) * hn * k[i]

    return np.squeeze(bn) 

def test_spherical_jn_at_zero(self):
        # http://dlmf.nist.gov/10.52.E1
        # But note that n = 0 is a special case: j0 = sin(x)/x -> 1
        n = np.array([0, 1, 2, 5, 10, 100])
        x = 0
        assert_allclose(spherical_jn(n, x), np.array([1, 0, 0, 0, 0, 0])) 

def test_spherical_jn_large_arg_2(self):
        # https://github.com/scipy/scipy/issues/1641
        # Reference value computed using mpmath, via
        # besselj(n + mpf(1)/2, z)*sqrt(pi/(2*z))
        assert_allclose(spherical_jn(2, 10000), 3.0590002633029811e-05) 

def test_spherical_jn_large_arg_1(self):
        # https://github.com/scipy/scipy/issues/2165
        # Reference value computed using mpmath, via
        # besselj(n + mpf(1)/2, z)*sqrt(pi/(2*z))
        assert_allclose(spherical_jn(2, 3350.507), -0.00029846226538040747) 

def test_spherical_jn_inf_real(self):
        # http://dlmf.nist.gov/10.52.E3
        n = 6
        x = np.array([-inf, inf])
        assert_allclose(spherical_jn(n, x), np.array([0, 0])) 

def test_spherical_jn_recurrence_real(self):
        # http://dlmf.nist.gov/10.51.E1
        n = np.array([1, 2, 3, 7, 12])
        x = 0.12
        assert_allclose(spherical_jn(n - 1, x) + spherical_jn(n + 1,x),
                        (2*n + 1)/x*spherical_jn(n, x)) 

def test_spherical_jn_recurrence_complex(self):
        # http://dlmf.nist.gov/10.51.E1
        n = np.array([1, 2, 3, 7, 12])
        x = 1.1 + 1.5j
        assert_allclose(spherical_jn(n - 1, x) + spherical_jn(n + 1, x),
                        (2*n + 1)/x*spherical_jn(n, x)) 

def test_spherical_jn_exact(self):
        # http://dlmf.nist.gov/10.49.E3
        # Note: exact expression is numerically stable only for small
        # n or z >> n.
        x = np.array([0.12, 1.23, 12.34, 123.45, 1234.5])
        assert_allclose(spherical_jn(2, x),
                        (-1/x + 3/x**3)*sin(x) - 3/x**2*cos(x)) 

def test_riccati_jn(self):
        N, x = 2, 0.2
        S = np.empty((N, N))
        for n in range(N):
            j = special.spherical_jn(n, x)
            jp = special.spherical_jn(n, x, derivative=True)
            S[0,n] = x*j
            S[1,n] = x*jp + j
        assert_array_almost_equal(S, special.riccati_jn(n, x), 8) 

def spherical_head_filter(max_order, full_order, kr, is_tapering=False):
    """Generate coloration compensation filter of specified maximum SH order.

    Parameters
    ----------
    max_order : int
        Maximum order
    full_order : int
        Full order necessary to expand sound field in entire modal range
    kr : array_like
        Vector of corresponding wave numbers
    is_tapering : bool, optional
        If set, spherical head filter will be adapted applying a Hann window, according to [2]

    Returns
    -------
    G_SHF : array_like
       Vector of frequency domain filter of shape [NFFT / 2 + 1]

    References
    ----------
    .. [1] Ben-Hur, Z., Brinkmann, F., Sheaffer, J., et al. (2017). "Spectral equalization in binaural signals
           represented by order-truncated spherical harmonics. The Journal of the Acoustical Society of America".

    .. [2] Hold, Christoph, Hannes Gamper, Ville Pulkki, Nikunj Raghuvanshi, and Ivan J. Tashev (2019). “Improving
           Binaural Ambisonics Decoding by Spherical Harmonics Domain Tapering and Coloration Compensation.”
    """

    # noinspection PyShadowingNames
    def pressure_on_sphere(max_order, kr, taper_weights=None):
        """
        Calculate the diffuse field pressure frequency response of a spherical scatterer, up to the specified SH order.
        If tapering weights are specified, pressure on the sphere function will be adapted.
        """
        if taper_weights is None:
            taper_weights = _np.ones(max_order + 1)  # no weighting

        p = _np.zeros_like(kr)
        for order in range(max_order + 1):
            # Calculate mode strength b_n(kr) for an incident plane wave on sphere according to [1, Eq.(9)]
            b_n = 4 * _np.pi * 1j ** order * (spherical_jn(order, kr) -
                                              (spherical_jn(order, kr, True) /
                                               dsphankel2(order, kr)) * sphankel2(order, kr))
            p += (2 * order + 1) * _np.abs(b_n) ** 2 * taper_weights[order]

        # according to [1, Eq.(11)]
        return 1 / (4 * _np.pi) * _np.sqrt(p)

    # according to [1, Eq.(12)].
    taper_weights = tapering_window(max_order) if is_tapering else None
    G_SHF = pressure_on_sphere(full_order, kr) / pressure_on_sphere(max_order, kr, taper_weights)

    G_SHF[G_SHF == 0] = 1e-12  # catch zeros
    G_SHF[_np.isnan(G_SHF)] = 1  # catch NaNs
    return G_SHF 

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

    Computes the radial weighting functions for diferent array types
    (cf. eq.(2.62), Rafaely 2015).

    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, 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.spherical_jn(n, x)
        if setup == 'open':
            bn = jn
        elif setup == 'card':
            bn = jn - 1j * special.spherical_jn(n, x, derivative=True)
        elif setup == 'rigid':
            if x == 0:
                # hn(x)/hn'(x) -> 0 for x -> 0
                bn = jn
            else:
                jnd = special.spherical_jn(n, x, derivative=True)
                hn = jn - 1j * special.spherical_yn(n, x)
                hnd = jnd - 1j * special.spherical_yn(n, x, derivative=True)
                bn = jn - jnd/hnd*hn
        else:
            raise ValueError('setup must be either: open, card or rigid')
        bns[i, :] = bn
    return np.squeeze(bns)