def decoding3dMaxRe(l) :
	if l is 0 : return numpy.array([1])
	import sympy
	x=sympy.Symbol("x")
	maxRoot = max((float(root) for root in sympy.solve(sympy.legendre(l+1,x),x)))
	return numpy.array([1.]+ [float(sympy.legendre(m, maxRoot))
		for m in xrange(1,l+1)])

def test_assoc_legendre():
    Plm = assoc_legendre
    Q = sqrt(1 - x**2)

    assert Plm(0, 0, x) == 1
    assert Plm(1, 0, x) == x
    assert Plm(1, 1, x) == -Q
    assert Plm(2, 0, x) == (3*x**2 - 1)/2
    assert Plm(2, 1, x) == -3*x*Q
    assert Plm(2, 2, x) == 3*Q**2
    assert Plm(3, 0, x) == (5*x**3 - 3*x)/2
    assert Plm(3, 1, x).expand() == (( 3*(1 - 5*x**2)/2 ).expand() * Q).expand()
    assert Plm(3, 2, x) == 15*x * Q**2
    assert Plm(3, 3, x) == -15 * Q**3

    # negative m
    assert Plm(1, -1, x) == -Plm(1, 1, x)/2
    assert Plm(2, -2, x) == Plm(2, 2, x)/24
    assert Plm(2, -1, x) == -Plm(2, 1, x)/6
    assert Plm(3, -3, x) == -Plm(3, 3, x)/720
    assert Plm(3, -2, x) == Plm(3, 2, x)/120
    assert Plm(3, -1, x) == -Plm(3, 1, x)/12

    n = Symbol("n")
    m = Symbol("m")

    X = Plm(n, m, x)
    assert isinstance(X, assoc_legendre)

    assert Plm(n, 0, x) == legendre(n, x)

    raises(ValueError, lambda: Plm(-1, 0, x))
    raises(ValueError, lambda: Plm(0, 1, x))

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    print "creating"
    e = legendre(p, x).diff(x)
    e = Poly(e, x)
    print "polydone"
    if e == 0:
        return []
    print "roots"
    if e == 1:
        r = []
    else:
        with workdps(40):
            r, err = polyroots(e.all_coeffs(), error=True)
        #if err > 1e-40:
        #    raise Exception("Internal Error: Root is not precise")
    print "done"
    p = []
    p.append("-1.0")
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    p.append("1.0")
    return p

def test_jacobi():
    n = Symbol("n")
    a = Symbol("a")
    b = Symbol("b")

    assert jacobi(0, a, b, x) == 1
    assert jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)

    assert jacobi(n, a, a, x) == RisingFactorial(a + 1, n)*gegenbauer(n, a + S(1)/2, x)/RisingFactorial(2*a + 1, n)
    assert jacobi(n, a, -a, x) == ((-1)**a*(-x + 1)**(-a/2)*(x + 1)**(a/2)*assoc_legendre(n, a, x)*
                                   factorial(-a + n)*gamma(a + n + 1)/(factorial(a + n)*gamma(n + 1)))
    assert jacobi(n, -b, b, x) == ((-x + 1)**(b/2)*(x + 1)**(-b/2)*assoc_legendre(n, b, x)*
                                   gamma(-b + n + 1)/gamma(n + 1))
    assert jacobi(n, 0, 0, x) == legendre(n, x)
    assert jacobi(n, S.Half, S.Half, x) == RisingFactorial(S(3)/2, n)*chebyshevu(n, x)/factorial(n + 1)
    assert jacobi(n, -S.Half, -S.Half, x) == RisingFactorial(S(1)/2, n)*chebyshevt(n, x)/factorial(n)

    X = jacobi(n, a, b, x)
    assert isinstance(X, jacobi)

    assert jacobi(n, a, b, -x) == (-1)**n*jacobi(n, b, a, x)
    assert jacobi(n, a, b, 0) == 2**(-n)*gamma(a + n + 1)*hyper((-b - n, -n), (a + 1,), -1)/(factorial(n)*gamma(a + 1))
    assert jacobi(n, a, b, 1) == RisingFactorial(a + 1, n)/factorial(n)

    m = Symbol("m", positive=True)
    assert jacobi(m, a, b, oo) == oo*RisingFactorial(a + b + m + 1, m)

    assert conjugate(jacobi(m, a, b, x)) == jacobi(m, conjugate(a), conjugate(b), conjugate(x))

    assert diff(jacobi(n,a,b,x), n) == Derivative(jacobi(n, a, b, x), n)
    assert diff(jacobi(n,a,b,x), x) == (a/2 + b/2 + n/2 + S(1)/2)*jacobi(n - 1, a + 1, b + 1, x)

def legendre_norm(i, x):
    """
    Returns the normalized integrated Legendre polynomial.
    """
    f = legendre(i, x)
    n = sqrt(integrate(f**2, (x, -1, 1)))
    return f/n

    def polyrs(self, uhat):
        "Compute the polynomial for the response surface."
        import sympy

        dim = len(self.params)
        var = []
        for i, p in enumerate(self.params):
            var.append(sympy.Symbol(str(p.name)))

        chaos = chaos_sequence(dim, self.level)
        chaos = np.int_(chaos)

        for d in range(0, dim):
            poly = np.array(map(lambda x: sympy.legendre(x, var[d]), chaos[:, d]))
            if d == 0:
                s = poly
            else:
                s *= poly

        eqn = np.sum(uhat * s)
        for i, p in enumerate(self.params):
            pmin, pmax = p.pdf.srange
            c = (pmax + pmin) / 2.0
            s = (pmax - pmin) / 2.0
            eqn = eqn.subs(var[i], (var[i] - c)/s)
        return sympy.sstr(eqn.expand().evalf())

def test_gegenbauer():
    n = Symbol("n")
    a = Symbol("a")

    assert gegenbauer(0, a, x) == 1
    assert gegenbauer(1, a, x) == 2*a*x
    assert gegenbauer(2, a, x) == -a + x**2*(2*a**2 + 2*a)
    assert gegenbauer(3, a, x) == \
        x**3*(4*a**3/3 + 4*a**2 + 8*a/3) + x*(-2*a**2 - 2*a)

    assert gegenbauer(-1, a, x) == 0
    assert gegenbauer(n, S(1)/2, x) == legendre(n, x)
    assert gegenbauer(n, 1, x) == chebyshevu(n, x)
    assert gegenbauer(n, -1, x) == 0

    X = gegenbauer(n, a, x)
    assert isinstance(X, gegenbauer)

    assert gegenbauer(n, a, -x) == (-1)**n*gegenbauer(n, a, x)
    assert gegenbauer(n, a, 0) == 2**n*sqrt(pi) * \
        gamma(a + n/2)/(gamma(a)*gamma(-n/2 + S(1)/2)*gamma(n + 1))
    assert gegenbauer(n, a, 1) == gamma(2*a + n)/(gamma(2*a)*gamma(n + 1))

    assert gegenbauer(n, Rational(3, 4), -1) == zoo

    m = Symbol("m", positive=True)
    assert gegenbauer(m, a, oo) == oo*RisingFactorial(a, m)

    assert conjugate(gegenbauer(n, a, x)) == gegenbauer(n, conjugate(a), conjugate(x))

    assert diff(gegenbauer(n, a, x), n) == Derivative(gegenbauer(n, a, x), n)
    assert diff(gegenbauer(n, a, x), x) == 2*a*gegenbauer(n - 1, a + 1, x)

def test_manualintegrate_orthogonal_poly():
    n = symbols('n')
    a, b = 7, S(5)/3
    polys = [jacobi(n, a, b, x), gegenbauer(n, a, x), chebyshevt(n, x),
        chebyshevu(n, x), legendre(n, x), hermite(n, x), laguerre(n, x),
        assoc_laguerre(n, a, x)]
    for p in polys:
        integral = manualintegrate(p, x)
        for deg in [-2, -1, 0, 1, 3, 5, 8]:
            # some accept negative "degree", some do not
            try:
                p_subbed = p.subs(n, deg)
            except ValueError:
                continue
            assert (integral.subs(n, deg).diff(x) - p_subbed).expand() == 0

        # can also integrate simple expressions with these polynomials
        q = x*p.subs(x, 2*x + 1)
        integral = manualintegrate(q, x)
        for deg in [2, 4, 7]:
            assert (integral.subs(n, deg).diff(x) - q.subs(n, deg)).expand() == 0

        # cannot integrate with respect to any other parameter
        t = symbols('t')
        for i in range(len(p.args) - 1):
            new_args = list(p.args)
            new_args[i] = t
            assert isinstance(manualintegrate(p.func(*new_args), t), Integral)

def legendre_int(i, x):
    """
    Returns the normalized integrated Legendre polynomial.
    """
    y = Symbol("y", dummy=True)
    f = legendre(i, y)
    n = sqrt(integrate(f**2, (y, -1, 1)))
    return integrate(f, (y, -1, x))/n

def acceptance( cosTheta1, phi, cosTheta2, i_max, j_max, k_max, numEvents, c):
	returnValue = 0.
        for i in range(i_max+1):
		for k in range(3):
			for j in range(0, 5, 2): # must have l >= k
 				if (j < k): continue
 				P_i  = N(legendre(i,    cosTheta2))
				Y_jk = abs(N(Zlm(j, k, acos(cosTheta1), phi)))
				#print returnValue, P_i, Y_jk, c[i][j][k]
 				returnValue += c[i][k][j]*(P_i * Y_jk);
	return returnValue

def gauss_lobatto_jacobi_quadruature_points_weights(n=5):
    """
    compute Gauss-Lobatto-Jacobi (GLJ) quadrature weights and points [-1, 1]
    using the sympy symbolic package and analytical definitions of these points

    somehow seems not to work for even numbers n >= 10, I assume this is a bug
    in sympy

    :param n: number of integration points (order + 1)
    :type n: integer

    :returns: tuple of two numpy arrays of floats containing the points and
        weights
    """
    x = sp.symbols('x')

    # GLL points are defined as the roots of the derivative of a sum of
    # Legendre Polynomials and include the boundaries
    Pn_bar = (sp.legendre(n-1, x) + sp.legendre(n, x)) / (1 + x)
    dPn_bar = sp.cancel((1 - x ** 2) * sp.diff(Pn_bar, x))

    # evaluate and sort
    pointsl = []
    for i in np.arange(n):
        p = sp.RootOf(dPn_bar, i)
        pointsl.append(p.evalf())

    pointsl.sort()

    # weights
    weightsl = []
    for p in pointsl:
        if p == -1.:
            wi = 8. / (n ** 2 * (n + 1) * (n - 1))
        else:
            wi = 4. / ((n + 1) * (n - 1) * Pn_bar.subs(x, p) ** 2)
        weightsl.append(wi)

    return np.array(pointsl, dtype='float'), np.array(weightsl, dtype='float')

    def test_five_points(self):
        import numpy as np
        import sympy

        t_symb = sympy.Symbol('t')
        num_points = 5
        p4 = sympy.legendre(num_points - 1, t_symb)
        p4_prime = p4.diff(t_symb)
        start = -1.0
        stop = 1.0
        all_nodes = self._call_func_under_test(start, stop, num_points)
        self.assertEqual(all_nodes[0], start)
        self.assertEqual(all_nodes[-1], stop)
        inner_nodes = all_nodes[1:-1]
        # Make sure the computed roots are actually roots.
        self.assertTrue(np.allclose(
            [float(p4_prime.subs({t_symb: root}))
             for root in inner_nodes], 0))

def gauss_lobatto_points(p):
    """
    Returns the list of Gauss-Lobatto points of the order 'p'.
    """
    x = Symbol("x")
    e = (1-x**2)*legendre(p, x).diff(x)
    e = Poly(e, x)
    if e == 0:
        return []
    with workdps(40):
        r, err = polyroots(e.all_coeffs(), error=True)
    if err > 1e-40:
        raise Exception("Internal Error: Root is not precise")
    p = []
    for x in r:
        if abs(x) < 1e-40: x = 0
        p.append(str(x))
    return p

def gauss_quadruature_points_weights(n=5):
    """
    compute Gaussian quadrature weights and points [-1, 1] using the sympy
    symbolic package and analytical definitions of these points

    somehow seems not to work for even numbers n >= 10, I assume this is a bug
    in sympy

    :param n: number of integration points (order + 1)
    :type n: integer

    :returns: tuple of two numpy arrays of floats containing the points and
        weights
    """
    x = sp.symbols('x')

    # Gauss points are defined as the roots of the Legendre Polynomials
    Pn = sp.legendre(n, x)

    # evaluate and sort
    pointsl = []
    for i in np.arange(n):
        p = sp.RootOf(Pn, i)
        pointsl.append(p.evalf())

    pointsl.sort()

    # weights are found through the derivative of the Legendre Polynomials
    dPn = sp.diff(Pn)

    # evaluate
    weightsl = []
    for p in pointsl:
        wi = 2. / ((1 - p ** 2) * dPn.subs(x, p) ** 2)
        weightsl.append(wi)

    return np.array(pointsl, dtype='float'), np.array(weightsl, dtype='float')

    def __init__(self, n=None, verbose=False):
        self.n = n
        prec = dps = Float.getdps()
        self.prec = prec

        # Reuse old nodes
        if n in _gausscache:
            cacheddps, xs, ws = _gausscache[n]
            if cacheddps >= dps:
                self.x = [x for x in xs]
                self.w = [w for w in ws]
                return

        if verbose:
            print ("calculating nodes for degree-%i Gauss-Legendre "
                "quadrature..." % n)

        Float.store()
        Float.setdps(2*prec + 5)

        self.x = [None] * n
        self.w = [None] * n

        pf = polyfunc(legendre(n, 'x'), True)

        for k in xrange(n//2 + 1):
            if verbose and k % 4 == 0:
                print "  node", k, "of", n//2
            x, w = _gaussnode(pf, k, n)
            self.x[k] = x
            self.x[n-k-1] = -x
            self.w[k] = self.w[n-k-1] = w

        _gausscache[n] = (dps, self.x, self.w)

        Float.revert()

def test_latex_functions():
    assert latex(exp(x)) == "e^{x}"
    assert latex(exp(1) + exp(2)) == "e + e^{2}"

    f = Function("f")
    assert latex(f(x)) == "\\operatorname{f}{\\left (x \\right )}"

    beta = Function("beta")

    assert latex(beta(x)) == r"\beta{\left (x \right )}"
    assert latex(sin(x)) == r"\sin{\left (x \right )}"
    assert latex(sin(x), fold_func_brackets=True) == r"\sin {x}"
    assert latex(sin(2 * x ** 2), fold_func_brackets=True) == r"\sin {2 x^{2}}"
    assert latex(sin(x ** 2), fold_func_brackets=True) == r"\sin {x^{2}}"

    assert latex(asin(x) ** 2) == r"\operatorname{asin}^{2}{\left (x \right )}"
    assert latex(asin(x) ** 2, inv_trig_style="full") == r"\arcsin^{2}{\left (x \right )}"
    assert latex(asin(x) ** 2, inv_trig_style="power") == r"\sin^{-1}{\left (x \right )}^{2}"
    assert latex(asin(x ** 2), inv_trig_style="power", fold_func_brackets=True) == r"\sin^{-1} {x^{2}}"

    assert latex(factorial(k)) == r"k!"
    assert latex(factorial(-k)) == r"\left(- k\right)!"

    assert latex(subfactorial(k)) == r"!k"
    assert latex(subfactorial(-k)) == r"!\left(- k\right)"

    assert latex(factorial2(k)) == r"k!!"
    assert latex(factorial2(-k)) == r"\left(- k\right)!!"

    assert latex(binomial(2, k)) == r"{\binom{2}{k}}"

    assert latex(FallingFactorial(3, k)) == r"{\left(3\right)}_{\left(k\right)}"
    assert latex(RisingFactorial(3, k)) == r"{\left(3\right)}^{\left(k\right)}"

    assert latex(floor(x)) == r"\lfloor{x}\rfloor"
    assert latex(ceiling(x)) == r"\lceil{x}\rceil"
    assert latex(Min(x, 2, x ** 3)) == r"\min\left(2, x, x^{3}\right)"
    assert latex(Min(x, y) ** 2) == r"\min\left(x, y\right)^{2}"
    assert latex(Max(x, 2, x ** 3)) == r"\max\left(2, x, x^{3}\right)"
    assert latex(Max(x, y) ** 2) == r"\max\left(x, y\right)^{2}"
    assert latex(Abs(x)) == r"\lvert{x}\rvert"
    assert latex(re(x)) == r"\Re{x}"
    assert latex(re(x + y)) == r"\Re{x} + \Re{y}"
    assert latex(im(x)) == r"\Im{x}"
    assert latex(conjugate(x)) == r"\overline{x}"
    assert latex(gamma(x)) == r"\Gamma\left(x\right)"
    assert latex(Order(x)) == r"\mathcal{O}\left(x\right)"
    assert latex(lowergamma(x, y)) == r"\gamma\left(x, y\right)"
    assert latex(uppergamma(x, y)) == r"\Gamma\left(x, y\right)"

    assert latex(cot(x)) == r"\cot{\left (x \right )}"
    assert latex(coth(x)) == r"\coth{\left (x \right )}"
    assert latex(re(x)) == r"\Re{x}"
    assert latex(im(x)) == r"\Im{x}"
    assert latex(root(x, y)) == r"x^{\frac{1}{y}}"
    assert latex(arg(x)) == r"\arg{\left (x \right )}"
    assert latex(zeta(x)) == r"\zeta\left(x\right)"

    assert latex(zeta(x)) == r"\zeta\left(x\right)"
    assert latex(zeta(x) ** 2) == r"\zeta^{2}\left(x\right)"
    assert latex(zeta(x, y)) == r"\zeta\left(x, y\right)"
    assert latex(zeta(x, y) ** 2) == r"\zeta^{2}\left(x, y\right)"
    assert latex(dirichlet_eta(x)) == r"\eta\left(x\right)"
    assert latex(dirichlet_eta(x) ** 2) == r"\eta^{2}\left(x\right)"
    assert latex(polylog(x, y)) == r"\operatorname{Li}_{x}\left(y\right)"
    assert latex(polylog(x, y) ** 2) == r"\operatorname{Li}_{x}^{2}\left(y\right)"
    assert latex(lerchphi(x, y, n)) == r"\Phi\left(x, y, n\right)"
    assert latex(lerchphi(x, y, n) ** 2) == r"\Phi^{2}\left(x, y, n\right)"

    assert latex(Ei(x)) == r"\operatorname{Ei}{\left (x \right )}"
    assert latex(Ei(x) ** 2) == r"\operatorname{Ei}^{2}{\left (x \right )}"
    assert latex(expint(x, y) ** 2) == r"\operatorname{E}_{x}^{2}\left(y\right)"
    assert latex(Shi(x) ** 2) == r"\operatorname{Shi}^{2}{\left (x \right )}"
    assert latex(Si(x) ** 2) == r"\operatorname{Si}^{2}{\left (x \right )}"
    assert latex(Ci(x) ** 2) == r"\operatorname{Ci}^{2}{\left (x \right )}"
    assert latex(Chi(x) ** 2) == r"\operatorname{Chi}^{2}{\left (x \right )}"

    assert latex(jacobi(n, a, b, x)) == r"P_{n}^{\left(a,b\right)}\left(x\right)"
    assert latex(jacobi(n, a, b, x) ** 2) == r"\left(P_{n}^{\left(a,b\right)}\left(x\right)\right)^{2}"
    assert latex(gegenbauer(n, a, x)) == r"C_{n}^{\left(a\right)}\left(x\right)"
    assert latex(gegenbauer(n, a, x) ** 2) == r"\left(C_{n}^{\left(a\right)}\left(x\right)\right)^{2}"
    assert latex(chebyshevt(n, x)) == r"T_{n}\left(x\right)"
    assert latex(chebyshevt(n, x) ** 2) == r"\left(T_{n}\left(x\right)\right)^{2}"
    assert latex(chebyshevu(n, x)) == r"U_{n}\left(x\right)"
    assert latex(chebyshevu(n, x) ** 2) == r"\left(U_{n}\left(x\right)\right)^{2}"
    assert latex(legendre(n, x)) == r"P_{n}\left(x\right)"
    assert latex(legendre(n, x) ** 2) == r"\left(P_{n}\left(x\right)\right)^{2}"
    assert latex(assoc_legendre(n, a, x)) == r"P_{n}^{\left(a\right)}\left(x\right)"
    assert latex(assoc_legendre(n, a, x) ** 2) == r"\left(P_{n}^{\left(a\right)}\left(x\right)\right)^{2}"
    assert latex(laguerre(n, x)) == r"L_{n}\left(x\right)"
    assert latex(laguerre(n, x) ** 2) == r"\left(L_{n}\left(x\right)\right)^{2}"
    assert latex(assoc_laguerre(n, a, x)) == r"L_{n}^{\left(a\right)}\left(x\right)"
    assert latex(assoc_laguerre(n, a, x) ** 2) == r"\left(L_{n}^{\left(a\right)}\left(x\right)\right)^{2}"
    assert latex(hermite(n, x)) == r"H_{n}\left(x\right)"
    assert latex(hermite(n, x) ** 2) == r"\left(H_{n}\left(x\right)\right)^{2}"

    # Test latex printing of function names with "_"
    assert latex(polar_lift(0)) == r"\operatorname{polar\_lift}{\left (0 \right )}"
    assert latex(polar_lift(0) ** 3) == r"\operatorname{polar\_lift}^{3}{\left (0 \right )}"

def gen_lobatto(max_order):
    assert max_order > 2

    x = sm.symbols('x')

    lobs = [0, 1]
    lobs[0] = (1 - x) / 2
    lobs[1] = (1 + x) / 2

    dlobs = [lob.diff('x') for lob in lobs]

    legs = [sm.legendre(0, 'y')]
    clegs = [sm.ccode(legs[0])]
    dlegs = [sm.legendre(0, 'y').diff('y')]
    cdlegs = [sm.ccode(dlegs[0])]

    clobs = [sm.ccode(lob) for lob in lobs]
    cdlobs = [sm.ccode(dlob) for dlob in dlobs]

    denoms = [] # for lobs.

    for ii in range(2, max_order + 1):
        coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % ii)
        leg = sm.legendre(ii - 1, 'y')

        pleg = leg.as_poly()
        coefs = pleg.all_coeffs()
        denom = max(sm.denom(val) for val in coefs)

        cleg = sm.ccode(sm.horner(leg*denom)/denom)

        dleg = leg.diff('y')
        cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

        lob = sm.simplify(coef * sm.integrate(leg, ('y', -1, x)))
        lobnc = sm.simplify(sm.integrate(leg, ('y', -1, x)))

        plobnc = lobnc.as_poly()
        coefs = plobnc.all_coeffs()
        denom = sm.denom(coef) * max(sm.denom(val) for val in coefs)

        clob = sm.ccode(sm.horner(lob*denom)/denom)

        dlob = lob.diff('x')
        cdlob = sm.ccode(sm.horner(dlob*denom)/denom)

        legs.append(leg)
        clegs.append(cleg)
        dlegs.append(dleg)
        cdlegs.append(cdleg)
        lobs.append(lob)
        clobs.append(clob)
        dlobs.append(dlob)
        cdlobs.append(cdlob)
        denoms.append(denom)

    coef = sm.sympify('sqrt(2 * (2 * %s - 1)) / 2' % (max_order + 1))
    leg = sm.legendre(max_order, 'y')

    pleg = leg.as_poly()
    coefs = pleg.all_coeffs()
    denom = max(sm.denom(val) for val in coefs)

    cleg = sm.ccode(sm.horner(leg*denom)/denom)

    dleg = leg.diff('y')
    cdleg = sm.ccode(sm.horner(dleg*denom)/denom)

    legs.append(leg)
    clegs.append(cleg)
    dlegs.append(dleg)
    cdlegs.append(cdleg)

    kerns = []
    ckerns = []
    dkerns = []
    cdkerns = []
    for ii, lob in enumerate(lobs[2:]):
        kern = sm.simplify(lob / (lobs[0] * lobs[1]))
        dkern = kern.diff('x')

        denom = denoms[ii] / 4
        ckern = sm.ccode(sm.horner(kern*denom)/denom)
        cdkern = sm.ccode(sm.horner(dkern*denom)/denom)

        kerns.append(kern)
        ckerns.append(ckern)
        dkerns.append(dkern)
        cdkerns.append(cdkern)

    return (legs, clegs, dlegs, cdlegs,
            lobs, clobs, dlobs, cdlobs,
            kerns, ckerns, dkerns, cdkerns,
            denoms)

	def __init__(self, decoding) :
		self.x = sympy.symbols("x")
		self.pattern = sum( (decoding[l] * sh_normalization2(l,0) * sympy.legendre(l,self.x) for l in xrange(0,len(decoding)) ))
		print self.pattern

a, z = sympy.symbols("a z")
a0, z0 = 0, numpy.pi/2
epsilon = a*1e-15 + z*1e-15
order = 5
nSamples = 70
angles = numpy.arange(0,nSamples+1)*2*numpy.pi/nSamples
w = 2*math.pi*spectralRange/spectrumBins * numpy.arange(spectrumBins)

if False:
	print "Checking that the azimuthal formula is equivalent to the zenital formula"
	for l in xrange(0,order+1) :
	#	pattern1 = sum([sh_normalization2(l,m)*sympy.assoc_legendre(l,m,sympy.cos(z))*sympy.assoc_legendre(l,m,1) for m in xrange(0,l+1)])
	#	pattern1 = sh_normalization2(l,0)*sympy.assoc_legendre(l,0,sympy.cos(z))
		pattern1 = sh_normalization2(l,0)*sympy.legendre(l,sympy.cos(z))
		pylab.polar(angles, [ pattern1.subs(z,angle)/(2*l+1) for angle in angles ], label="Zenital %s"%l)

		pattern2 = sum([sh_normalization2(l,m)*sympy.cos(m*z)*sympy.assoc_legendre(l,m,0)**2 for m in xrange(0,l+1)])
		pylab.polar(angles, [ pattern2.subs(z,angle)/(2*l+1) for angle in angles ], label="Azimuthal %s"%l)
	pylab.title("Zenital vs Azimuthal variation",horizontalalignment='center', verticalalignment='baseline', position=(.5,-.1))
	pylab.rgrids(numpy.arange(.4,1,.2),angle=220)
	pylab.legend(loc=2)
	pylab.savefig(figurePath(__file__,"pdf"))
	pylab.show()
# We take the azimuthal formula which is faster
print "Computing component patterns"
patternComponents = [
	sh_normalization2(l,0)*sympy.legendre(l,sympy.cos(z))
	for l in xrange(0,order+1)
	]

def test_legendre():
    raises(ValueError, lambda: legendre(-1, x))
    assert legendre(0, x) == 1
    assert legendre(1, x) == x
    assert legendre(2, x) == ((3*x**2 - 1)/2).expand()
    assert legendre(3, x) == ((5*x**3 - 3*x)/2).expand()
    assert legendre(4, x) == ((35*x**4 - 30*x**2 + 3)/8).expand()
    assert legendre(5, x) == ((63*x**5 - 70*x**3 + 15*x)/8).expand()
    assert legendre(6, x) == ((231*x**6 - 315*x**4 + 105*x**2 - 5)/16).expand()

    assert legendre(10, -1) == 1
    assert legendre(11, -1) == -1
    assert legendre(10, 1) == 1
    assert legendre(11, 1) == 1
    assert legendre(10, 0) != 0
    assert legendre(11, 0) == 0

    assert roots(legendre(4, x), x) == {
        sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) - Rational(2, 35)*sqrt(30)): 1,
        sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
        -sqrt(Rational(3, 7) + Rational(2, 35)*sqrt(30)): 1,
    }

    n = Symbol("n")

    X = legendre(n, x)
    assert isinstance(X, legendre)

    assert legendre(-n, x) == legendre(n - 1, x)
    assert legendre(n, -x) == (-1)**n*legendre(n, x)

    assert conjugate(legendre(n, x)) == legendre(n, conjugate(x))

    assert diff(legendre(n, x), x) == \
        n*(x*legendre(n, x) - legendre(n - 1, x))/(x**2 - 1)
    assert diff(legendre(n, x), n) == Derivative(legendre(n, x), n)