Python special.spherical_jn方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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)) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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])) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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])) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import spherical_jn [as 别名]
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)