Python scipy.special.jacobi() Examples

The following are 5 code examples of scipy.special.jacobi(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.special , or try the search function .
Example #1
Source File: test_basic.py    From Computable with MIT License 6 votes vote down vote up 
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)
Example #2
Source File: test_basic.py    From GraphicDesignPatternByPython with MIT License 6 votes vote down vote up 
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)
Example #3
Source File: wigner_d.py    From lie_learn with MIT License 6 votes vote down vote up 
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
Example #4
Source File: test_jacobi.py    From prysm with MIT License 5 votes vote down vote up 
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_)
Example #5
Source File: wigner_d.py    From lie_learn with MIT License 5 votes vote down vote up 
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