def Theta_Equation(self, l="undef", m="undef", theta=theta):
        """ l - orbital quantum number
            m - magnetic quantum number """

        if l == "undef":
            l = self.l_val
        if m == "undef":
            m = self.m_val

        # Prevents integer division and float l,m
        l = int(l)
        m = int(m)

        mode = self.select_exec_mode(theta)

        gL = self.Generalized_Legendre(l, m, theta).subs(theta, sympy.cos(theta))

        if gL:
            rat_part = sympy.sqrt(
                Rational((2 * l + 1), 2) * Rational(factorial(l - sympy.abs(m)), factorial(l + sympy.abs(m)))
            )
            if mode == "numer":
                return (rat_part * gL).evalf()
            elif mode == "analit":
                return (rat_part * gL).expand()

        return False

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"$\operatorname{beta}\left(x\right)$"
    assert latex(sin(x)) == r"$\operatorname{sin}\left(x\right)$"
    assert latex(sin(x), fold_func_brackets=True) == r"$\operatorname{sin}x$"
    assert latex(sin(2 * x ** 2), fold_func_brackets=True) == r"$\operatorname{sin}2 x^{2}$"
    assert latex(sin(x ** 2), fold_func_brackets=True) == r"$\operatorname{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"$\operatorname{arcsin}^{2}\left(x\right)$"
    assert latex(asin(x) ** 2, inv_trig_style="power") == r"$\operatorname{sin}^{-1}\left(x\right)^{2}$"
    assert latex(asin(x ** 2), inv_trig_style="power", fold_func_brackets=True) == r"$\operatorname{sin}^{-1}x^{2}$"

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

    assert latex(floor(x)) == r"$\lfloor{x}\rfloor$"
    assert latex(ceiling(x)) == r"$\lceil{x}\rceil$"
    assert latex(abs(x)) == r"$\lvert{x}\rvert$"
    assert latex(re(x)) == r"$\Re{x}$"
    assert latex(im(x)) == r"$\Im{x}$"
    assert latex(conjugate(x)) == r"$\overline{x}$"
    assert latex(gamma(x)) == r"$\operatorname{\Gamma}\left(x\right)$"
    assert latex(Order(x)) == r"$\operatorname{\mathcal{O}}\left(x\right)$"

def g(yieldCurve, zeroRates,n, verbose):
    '''
        generates recursively the zero curve 
        expressions eval('(0.06/1.05)+(1.06/(1+x)**2)-1')
        solves these expressions to get the new rate
        for that period
    
    '''
    if len(zeroRates) >= len(yieldCurve):
        print "\n\n\t+zero curve boot strapped [%d iterations]" % (n)
        return
    else:
        legn = ''
        for i in range(0,len(zeroRates),1):
            if i == 0:
                legn = '%2.6f/(1+%2.6f)**%d'%(yieldCurve[n], zeroRates[i],i+1)
            else:
                legn = legn + ' +%2.6f/(1+%2.6f)**%d'%(yieldCurve[n], zeroRates[i],i+1)
        legn = legn + '+ (1+%2.6f)/(1+x)**%d-1'%(yieldCurve[n], n+1)
        # solve the expression for this iteration
        if verbose:
            print "-[%d] %s" % (n, legn.strip())
        rate1 = solve(eval(legn), x)
        # Abs here since some solutions can be complex
        rate1 = min([Real(abs(r)) for r in rate1])
        if verbose:
            print "-[%d] solution %2.6f" % (n, float(rate1))
        # stuff the new rate in the results, will be 
        # used by the next iteration
        zeroRates.append(rate1)
        g(yieldCurve, zeroRates,n+1, verbose)

def run_Lagrange_interp_abs_conv(N=[3, 6, 12, 24]):
    f = sm.abs(1-2*x)
    f = sm.sin(2*sm.pi*x)
    fn = sm.lambdify([x], f, modules='numpy')
    resolution = 50001
    import numpy as np
    xcoor = np.linspace(0, 1, resolution)
    fcoor = fn(xcoor)
    Einf = []
    E2 = []
    h = []
    for _N in N:
        psi, points = Lagrange_polynomials_01(x, _N)
        u = interpolation(f, psi, points)
        un = sm.lambdify([x], u, modules='numpy')
        ucoor = un(xcoor)
        e = fcoor - ucoor
        Einf.append(e.max())
        E2.append(np.sqrt(np.sum(e*e/e.size)))
        h.append(1./_N)
    print Einf
    print E2
    print h
    print N
    # Assumption: error = CN**(-N)
    print 'convergence rates:'
    for i in range(len(E2)):
        C1 = E2[i]/(N[i]**(-N[i]/2))
        C2 = Einf[i]/(N[i]**(-N[i]/2))
        print N[i], C1, C2

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"$\operatorname{beta}\left(x\right)$"
    assert latex(sin(x)) == r"$\operatorname{sin}\left(x\right)$"
    assert latex(sin(x), fold_func_brackets=True) == r"$\operatorname{sin}x$"
    assert latex(sin(2*x**2), fold_func_brackets=True) == \
    r"$\operatorname{sin}2 x^{2}$"

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

    assert latex(floor(x)) == r"$\lfloor{x}\rfloor$"
    assert latex(ceiling(x)) == r"$\lceil{x}\rceil$"
    assert latex(abs(x)) == r"$\lvert{x}\rvert$"
    assert latex(re(x)) == r"$\Re{x}$"
    assert latex(im(x)) == r"$\Im{x}$"
    assert latex(conjugate(x)) == r"$\overline{x}$"
    assert latex(gamma(x)) == r"$\operatorname{\Gamma}\left(x\right)$"
    assert latex(Order(x)) == r"$\operatorname{\mathcal{O}}\left(x\right)$"

def sp_derive():

    import sympy as sp

    vars = 'G_s, s_n, s_p_n, w_n, dw_n, ds_n, G_s, G_w, c, phi'

    syms = sp.symbols(vars)

    for var, sym in zip(vars.split(','), syms):
        globals()[var.strip()] = sym

    s_n1 = s_n + ds_n
    w_n1 = w_n + dw_n

    tau_trial = G_s * (s_n1 - s_p_n)

    print 'diff', sp.diff(tau_trial, ds_n)

    print tau_trial

    sig_n1 = G_w * w_n1

    print sig_n1

    tau_fr = (c + sig_n1 * sp.tan(phi)) * sp.Heaviside(sig_n1 - c / sp.tan(phi))

    print tau_fr

    d_tau_fr = sp.diff(tau_fr, dw_n)

    print d_tau_fr

    f_trial = sp.abs(tau_trial) - tau_fr

    print f_trial

    d_gamma = f_trial / G_s

    print 'd_gamma'
    sp.pretty_print(d_gamma)

    print 'd_gamma_s'
    sp.pretty_print(sp.diff(d_gamma, ds_n))

    print 'tau_n1'
    tau_n1 = sp.simplify(tau_trial - d_gamma * G_s * sp.sign(tau_trial))
    sp.pretty_print(tau_n1)

    print 'dtau_n1_w'
    dtau_n1_w = sp.diff(tau_n1, dw_n)
    sp.pretty_print(dtau_n1_w)

    print 'dtau_n1_s'
    dtau_n1_s = sp.diff(d_gamma * sp.sign(tau_trial), ds_n)
    print dtau_n1_s

    s_p_n1 = s_p_n + d_gamma * sp.sign(tau_trial)

    print s_p_n1

def chop(expr, tol = 1e-10):
    """suppress small numerical values"""

    expr = sp.expand(sp.sympify(expr))
    if expr.is_Symbol: return expr

    assert expr.is_Add

    return sp.Add(*[term for term in expr.as_Add() if sp.abs(term.as_coeff_terms()[0]) >= tol])

def run_Lagrange_leastsq_abs(N):
    """Least-squares with of Lagrange polynomials for |1-2x|."""
    f = sm.abs(1-2*x)
    # This f will lead to failure of sympy integrate, fallback on numerical int.
    psi, points = Lagrange_polynomials_01(x, N)
    Omega=[0, 1]
    u = least_squares(f, psi, Omega)
    comparison_plot(f, u, Omega, filename='Lagrange_ls_abs_%d.pdf' % (N+1),
                    plot_title='Least squares approximation by '\
                    'Lagrange polynomials of degree %d' % N)

def abs(*args, **kw):
    arg0 = args[0]
    if isinstance(arg0, (int, float, long)):
        return _math.abs(*args,**kw)
    elif isinstance(arg0, complex):
        return _cmath.abs(*args,**kw)
    elif isinstance(arg0, _sympy.Basic):
        return _sympy.abs(*args,**kw)
    else:
        return _numpy.abs(*args,**kw)

    def Generalized_Legendre(self, n=0, m=0, x=x):
        mode = self.select_exec_mode(x)

        if mode != "undef":
            legendre_polinomial = self.Legendre(n=n, x=x)
            rat_part = (1 - (x ** 2)) ** (sympy.abs(m) / 2.0)
            diff_part = sympy.diff(legendre_polinomial, x, abs(m))

            return (rat_part * diff_part).expand()
        else:
            return False

    def Phi_Equation(self, m="undef", phi=phi):
        if m == "undef":
            m = self.m_val

        m = int(m)
        mode = self.select_exec_mode(phi)

        if mode == "numer":
            phi_val = self.phi
            self.phi = Symbol("phi")
        elif mode == "custom_var":
            self.phi = Symbol(phi)
        elif mode == "undef":
            return False

        left_part = 1 / sympy.sqrt(2 * sympy.pi)
        right_part = sympy.cos(sympy.abs(m) * self.phi) if m >= 0 else sympy.sin(sympy.abs(m) * self.phi)

        if mode == "numer":
            return (left_part * right_part).evalf()
        elif mode == "analit" or mode == "custom_var":
            return (left_part * right_part).expand()

def l2_norm(f, lim):
    """
    Calculates L2 norm of the function "f", over the domain lim=(x, a, b).

    x ...... the independent variable in f over which to integrate
    a, b ... the limits of the interval

    Example:

    >>> l2_norm(1, (x, -1, 1))
    2**(1/2)
    >>> l2_norm(x, (x, -1, 1))
    1/3*6**(1/2)
    """
    return sqrt(integrate(abs(f)**2, lim))

def test_maxima_functions():
    assert parse_maxima('expand( (x+1)^2)') == x**2 + 2*x + 1
    assert parse_maxima('factor( x**2 + 2*x + 1)') == (x+1)**2
    assert parse_maxima('trigsimp(2*cos(x)^2 + sin(x)^2)') == 2 - sin(x)**2
    assert parse_maxima('trigexpand(sin(2*x)+cos(2*x))') == (-1) + 2*cos(x)**2 + 2*cos(x)*sin(x)
    assert parse_maxima('solve(x^2-4,x)') == [2, -2]
    assert parse_maxima('limit((1+1/x)^x,x,inf)') == E
    assert parse_maxima('limit(sqrt(-x)/x,x,0,minus)') == -oo
    assert parse_maxima('diff(x^x, x)') == x**x*(1 + log(x))
    assert parse_maxima('sum(k, k, 1, n)') == (n**2 +n)/2
    assert parse_maxima('product(k, k, 1, n)',
            name_dict=dict(
                n=Symbol('n',integer=True),
                k=Symbol('k',integer=True)
                )
            )  == factorial(n)
    assert parse_maxima('ratsimp((x^2-1)/(x+1))') == x-1
    assert abs( parse_maxima('float(sec(%pi/3) + csc(%pi/3))') - 3.154700538379252) < 10**(-5)

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"\operatorname{beta}\left(x\right)"
    assert latex(sin(x)) == r"\operatorname{sin}\left(x\right)"
    assert latex(sin(x), fold_func_brackets=True) == r"\operatorname{sin}x"
    assert latex(sin(2*x**2), fold_func_brackets=True) == \
    r"\operatorname{sin}2 x^{2}"
    assert latex(sin(x**2), fold_func_brackets=True) == \
    r"\operatorname{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"\operatorname{arcsin}^{2}\left(x\right)"
    assert latex(asin(x)**2,inv_trig_style="power") == \
        r"\operatorname{sin}^{-1}\left(x\right)^{2}"
    assert latex(asin(x**2),inv_trig_style="power",fold_func_brackets=True) == \
        r"\operatorname{sin}^{-1}x^{2}"

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

    assert latex(floor(x)) == r"\lfloor{x}\rfloor"
    assert latex(ceiling(x)) == r"\lceil{x}\rceil"
    assert latex(abs(x)) == r"\lvert{x}\rvert"
    assert latex(re(x)) == r"\Re{x}"
    assert latex(im(x)) == r"\Im{x}"
    assert latex(conjugate(x)) == r"\overline{x}"
    assert latex(gamma(x)) == r"\operatorname{\Gamma}\left(x\right)"
    assert latex(Order(x)) == r"\operatorname{\mathcal{O}}\left(x\right)"

def test_parser():
    assert abs(parse_maxima('float(1/3)') - 0.333333333) < 10**(-5)
    assert parse_maxima('13^26') == 91733330193268616658399616009
    assert parse_maxima('sin(%pi/2) + cos(%pi/3)') == Rational(3,2)
    assert parse_maxima('log(%e)') == 1

def single_cross_section_calc(theta, phi, rad, N_atoms_uc, atom_list, csection, tau, eig_list, wt,
                       temperature, eief = True, efixed = 14.7):

    temperature = temperature

    front_constant = ((GAMMA_R0_VALUE**2)*DEBYE_WALLER_VALUE/(2*pi*HBAR_VALUE)).evalf()
 
    kx = rad*np.sin(theta)*np.cos(phi)
    ky = rad*np.sin(theta)*np.sin(phi)
    kz = rad*np.cos(theta)
    kap = np.array([kx,ky,kz])

    kapunit = kap.copy()
    kapunit[0]=kap[0]/rad
    kapunit[1]=kap[1]/rad
    kapunit[2]=kap[2]/rad    

    ws=[]

    #csectempp = deepcopy(csection)  #good
    #csectempm = deepcopy(csection)  #good             
    
    csectempp = copy(csection)  #good
    csectempm = copy(csection)  #good
    
    
    cs_expr = sp.lambdify((KAP_SYM,Q_SYM,TAU_SYM,W_SYM,WQ_SYM,KAPXHAT_SYM,KAPYHAT_SYM,KAPZHAT_SYM), expr2, modules = "numpy")

    #res_np = expr_np(Svals[:,newaxis,newaxis,newaxis],xvals[newaxis,:,newaxis,newaxis],yvals[newaxis,newaxis,:,newaxis],zvals[newaxis,newaxis,newaxis,:])
    res_np2 = expr_np2(xvals[:,newaxis,newaxis],yvals[newaxis,:,newaxis],zvals[newaxis,newaxis,:])
    
    
    csectempp = csectempp.subs(sp.DiracDelta(KAP_SYM - Q_SYM - TAU_SYM),sp.S(1))
    csectempp = csectempp.subs(sp.DiracDelta(KAP_SYM + Q_SYM - TAU_SYM),sp.S(0))
    csectempm = csectempm.subs(sp.DiracDelta(KAP_SYM - Q_SYM - TAU_SYM),sp.S(0))
    csectempm = csectempm.subs(sp.DiracDelta(KAP_SYM + Q_SYM - TAU_SYM),sp.S(1))

    csectempp = csectempp.subs(KAPXHAT_SYM,kapunit[0])
    csectempp = csectempp.subs(KAPYHAT_SYM,kapunit[1])
    csectempp = csectempp.subs(KAPZHAT_SYM,kapunit[2])
    csectempm = csectempm.subs(KAPXHAT_SYM,kapunit[0])
    csectempm = csectempm.subs(KAPYHAT_SYM,kapunit[1])
    csectempm = csectempm.subs(KAPZHAT_SYM,kapunit[2])

    for eigi in range(len(eig_list)):
        #eigcsecp=deepcopy(csectempp)  #good
        #eigcsecm=deepcopy(csectempm)
        #eigtemp = deepcopy(eig_list[0][eigi])
        eigcsecp=copy(csectempp)  #good
        eigcsecm=copy(csectempm)
        eigtemp = copy(eig_list[0][eigi])
        eig_expr = sp.lambdify((S_SYM,KX_SYM,KY_SYM,KZ_SYM), sp.abs(eigtemp))#, modules = "numpy")
        eigval = eig_expr(S(1.0),kap[0],kap[1],kap[2]).evalf(chop=True)

        eigtemp = eigtemp.subs(S_SYM, sp.S(1.0))
        eigtemp = eigtemp.subs(KX_SYM, kap[0])
        eigtemp = eigtemp.subs(KY_SYM, kap[1])
        eigtemp = eigtemp.subs(KZ_SYM, kap[2])
        eigtemp = sp.abs(eigtemp.evalf(chop=True)) #works
        #eigtemp = chop(np.abs(eigtemp))

        nval = sp.Pow(sp.exp(eigtemp/(BOLTZ_VALUE*temperature))-1,-1).evalf() #works
        #nval = np.power(np.exp(np.abs(eigtemp)/(BOLTZ*temperature))-1,-1)
        for i in range(N_atoms_uc):
            nq = sp.Symbol('n%i'%(i,), real = True)
            eigcsecp = eigcsecp.subs(nq,nval)
            eigcsecm = eigcsecm.subs(nq,nval) 

        wvalp = eigtemp - wt
        wvalm = eigtemp + wt

        eigcsecp = eigcsecp.subs((W_SYM - WQ_SYM),wvalp)
        eigcsecp = eigcsecp.subs((W_SYM + WQ_SYM),wvalm)
        eigcsecp = sp.re(eigcsecp.evalf(chop = True)) # works
        #eigcsecp = chop(eigcsecp)
        eigcsecm = eigcsecm.subs((W_SYM - WQ_SYM),wvalp)
        eigcsecm = eigcsecm.subs((W_SYM + WQ_SYM),wvalm)
        eigcsecm = sp.re(eigcsecm.evalf(chop = True)) #works
        #eigcsecm = chop(eigcsecm)

        # eief == True => ei=efixed
        # eief == False => ef=efixed
        kpk = 0.
        if eief == True:
            ei = efixed
            ef = ei - eigtemp
            kpk = ef/ei
        else:
            ef = efixed
            ei = ef + eigtemp
            kpk = ef/ei

        ws.append(kpk*(eigcsecp+eigcsecm))
    
    csdata = sp.re(sum(ws))

    #Multiply data by front constants
    csdata = front_constant*csdata

    #Multiply by Form Factor
#.........这里部分代码省略.........

def test_abs():
    assert str(abs(x)) == "abs(x)"
    assert str(abs(Rational(1,6))) == "1/6"
    assert str(abs(Rational(-1,6))) == "1/6"

def SGN_S(x):
    return x/abs(x)

def test_ccode_exceptions():
    assert ccode(ceiling(x)) == "ceil(x)"
    assert ccode(abs(x)) == "fabs(x)"

    lf = lambdify_t(f)
    tt = np.arange(dt,t+dt,dt)
    return lf(tt).sum() * dt 


_gexpr = gamma_expr(5.4, 5.2) - 0.35 * gamma_expr(10.8, 7.35)
_gexpr = _gexpr / _getint(_gexpr)
_glover = lambdify_t(_gexpr)
glover = implemented_function('glover', _glover)
glovert = lambdify_t(glover(T))

# Derivative of Glover HRF

_dgexpr = _gexpr.diff(T)
dpos = sympy.Derivative((T >= 0), T)
_dgexpr = _dgexpr.subs(dpos, 0)
_dgexpr = _dgexpr / _getint(sympy.abs(_dgexpr))
_dglover = lambdify_t(_dgexpr)
dglover = implemented_function('dglover', _dglover)
dglovert = lambdify_t(dglover(T))

del(_glover); del(_gexpr); del(dpos); del(_dgexpr); del(_dglover)

# AFNI's HRF

_aexpr = ((T >= 0) * T)**8.6 * sympy.exp(-T/0.547)
_aexpr = _aexpr / _getint(_aexpr)
_afni = lambdify_t(_aexpr)
afni = implemented_function('afni', _afni)
afnit = lambdify_t(afni(T))