Python special.jacobi方法代码示例

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import jacobi [as 别名]
def test_jacobi(self):
        a = 5*rand() - 1
        b = 5*rand() - 1
        P0 = special.jacobi(0,a,b)
        P1 = special.jacobi(1,a,b)
        P2 = special.jacobi(2,a,b)
        P3 = special.jacobi(3,a,b)

        assert_array_almost_equal(P0.c,[1],13)
        assert_array_almost_equal(P1.c,array([a+b+2,a-b])/2.0,13)
        cp = [(a+b+3)*(a+b+4), 4*(a+b+3)*(a+2), 4*(a+1)*(a+2)]
        p2c = [cp[0],cp[1]-2*cp[0],cp[2]-cp[1]+cp[0]]
        assert_array_almost_equal(P2.c,array(p2c)/8.0,13)
        cp = [(a+b+4)*(a+b+5)*(a+b+6),6*(a+b+4)*(a+b+5)*(a+3),
              12*(a+b+4)*(a+2)*(a+3),8*(a+1)*(a+2)*(a+3)]
        p3c = [cp[0],cp[1]-3*cp[0],cp[2]-2*cp[1]+3*cp[0],cp[3]-cp[2]+cp[1]-cp[0]]
        assert_array_almost_equal(P3.c,array(p3c)/48.0,13) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import jacobi [as 别名]
def test_jacobi(self):
        a = 5*np.random.random() - 1
        b = 5*np.random.random() - 1
        P0 = special.jacobi(0,a,b)
        P1 = special.jacobi(1,a,b)
        P2 = special.jacobi(2,a,b)
        P3 = special.jacobi(3,a,b)

        assert_array_almost_equal(P0.c,[1],13)
        assert_array_almost_equal(P1.c,array([a+b+2,a-b])/2.0,13)
        cp = [(a+b+3)*(a+b+4), 4*(a+b+3)*(a+2), 4*(a+1)*(a+2)]
        p2c = [cp[0],cp[1]-2*cp[0],cp[2]-cp[1]+cp[0]]
        assert_array_almost_equal(P2.c,array(p2c)/8.0,13)
        cp = [(a+b+4)*(a+b+5)*(a+b+6),6*(a+b+4)*(a+b+5)*(a+3),
              12*(a+b+4)*(a+2)*(a+3),8*(a+1)*(a+2)*(a+3)]
        p3c = [cp[0],cp[1]-3*cp[0],cp[2]-2*cp[1]+3*cp[0],cp[3]-cp[2]+cp[1]-cp[0]]
        assert_array_almost_equal(P3.c,array(p3c)/48.0,13) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import jacobi [as 别名]
def wigner_d_naive_v2(l, m, n, beta):
    """
    Wigner d functions as defined in the SOFT 2.0 documentation.
    When approx_lim is set to a high value, this function appears to give
    identical results to Johann Goetz' wignerd() function.

    However, integration fails: does not satisfy orthogonality relations everywhere...
    """
    from scipy.special import jacobi

    if n >= m:
        xi = 1
    else:
        xi = (-1)**(n - m)

    mu = np.abs(m - n)
    nu = np.abs(n + m)
    s = l - (mu + nu) * 0.5

    sq = np.sqrt((np.math.factorial(s) * np.math.factorial(s + mu + nu))
                 / (np.math.factorial(s + mu) * np.math.factorial(s + nu)))
    sinb = np.sin(beta * 0.5) ** mu
    cosb = np.cos(beta * 0.5) ** nu
    P = jacobi(s, mu, nu)(np.cos(beta))
    return xi * sq * sinb * cosb * P 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import jacobi [as 别名]
def test_jacobi_1_4_match_scipy(n, alpha, beta):
    x = np.linspace(-1, 1, 32)
    prysm_ = pjac.jacobi(n=n, alpha=alpha, beta=beta, x=x)
    scipy_ = sps_jac(n=n, alpha=alpha, beta=beta)(x)
    assert np.allclose(prysm_, scipy_) 

# 需要导入模块: from scipy import special [as 别名]
# 或者: from scipy.special import jacobi [as 别名]
def wigner_d_naive(l, m, n, beta):
    """
    Numerically naive implementation of the Wigner-d function.
    This is useful for checking the correctness of other implementations.

    :param l: the degree of the Wigner-d function. l >= 0
    :param m: the order of the Wigner-d function. -l <= m <= l
    :param n: the order of the Wigner-d function. -l <= n <= l
    :param beta: the argument. 0 <= beta <= pi
    :return: d^l_mn(beta) in the TODO: what basis? complex, quantum(?), centered, cs(?)
    """
    from scipy.special import eval_jacobi
    try:
        from scipy.misc import factorial
    except:
        from scipy.special import factorial

    from sympy.functions.special.polynomials import jacobi, jacobi_normalized
    from sympy.abc import j, a, b, x
    from sympy import N
    #jfun = jacobi_normalized(j, a, b, x)
    jfun = jacobi(j, a, b, x)
    # eval_jacobi = lambda q, r, p, o: float(jfun.eval(int(q), int(r), int(p), float(o)))
    # eval_jacobi = lambda q, r, p, o: float(N(jfun, int(q), int(r), int(p), float(o)))
    eval_jacobi = lambda q, r, p, o: float(jfun.subs({j:int(q), a:int(r), b:int(p), x:float(o)}))

    mu = np.abs(m - n)
    nu = np.abs(m + n)
    s = l - (mu + nu) / 2
    xi = 1 if n >= m else (-1) ** (n - m)

    # print(s, mu, nu, np.cos(beta), type(s), type(mu), type(nu), type(np.cos(beta)))
    jac = eval_jacobi(s, mu, nu, np.cos(beta))
    z = np.sqrt((factorial(s) * factorial(s + mu + nu)) / (factorial(s + mu) * factorial(s + nu)))

    # print(l, m, n, beta, np.isfinite(mu), np.isfinite(nu), np.isfinite(s), np.isfinite(xi), np.isfinite(jac), np.isfinite(z))
    assert np.isfinite(mu) and np.isfinite(nu) and np.isfinite(s) and np.isfinite(xi) and np.isfinite(jac) and np.isfinite(z)
    assert np.isfinite(xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac)
    return xi * z * np.sin(beta / 2) ** mu * np.cos(beta / 2) ** nu * jac