def test_ode_solutions():
    # only a few examples here, the rest will be tested in the actual dsolve tests
    assert ode_renumber(constantsimp(C1*exp(2*x)+exp(x)*(C2+C3), x, 3), 'C', 1, 3) == \
        ode_renumber(C1*exp(x)+C2*exp(2*x), 'C', 1, 2)
    assert ode_renumber(constantsimp(Eq(f(x),I*C1*sinh(x/3) + C2*cosh(x/3)), x, 2),
        'C', 1, 2) == ode_renumber(Eq(f(x), C1*sinh(x/3) + C2*cosh(x/3)), 'C', 1, 2)
    assert ode_renumber(constantsimp(Eq(f(x),acos((-C1)/cos(x))), x, 1), 'C', 1, 1) == \
        Eq(f(x),acos(C1/cos(x)))
    assert ode_renumber(constantsimp(Eq(log(f(x)/C1) + 2*exp(x/f(x)), 0), x, 1),
        'C', 1, 1) ==  Eq(log(C1*f(x)) + 2*exp(x/f(x)), 0)
    assert ode_renumber(constantsimp(Eq(log(x*2**Rational(1,2)*(1/x)**Rational(1,2)*f(x)\
        **Rational(1,2)/C1) + x**2/(2*f(x)**2), 0), x, 1), 'C', 1, 1) == \
        Eq(log(C1*x*(1/x)**Rational(1,2)*f(x)**Rational(1,2)) + x**2/(2*f(x)**2), 0)
    assert ode_renumber(constantsimp(Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(x/C1) - \
        cos(f(x)/x)*exp(-f(x)/x)/2, 0), x, 1), 'C', 1, 1) == \
        Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(C1*x) - cos(f(x)/x)*exp(-f(x)/x)/2, 0)
    u2 = Symbol('u2')
    _a = Symbol('_a')
    assert ode_renumber(constantsimp(Eq(-Integral(-1/((1 - u2**2)**Rational(1,2)*u2), \
        (u2, _a, x/f(x))) + log(f(x)/C1), 0), x, 1), 'C', 1, 1) == \
        Eq(-Integral(-1/(u2*(1 - u2**2)**Rational(1,2)), (u2, _a, x/f(x))) + \
        log(C1*f(x)), 0)
    assert map(lambda i: ode_renumber(constantsimp(i, x, 1), 'C', 1, 1),
        [Eq(f(x), (-C1*x + x**2)**Rational(1,2)), Eq(f(x), -(-C1*x +
        x**2)**Rational(1,2))]) == [Eq(f(x), (C1*x + x**2)**Rational(1,2)),
        Eq(f(x), -(C1*x + x**2)**Rational(1,2))]

def test_acos_series():
    assert acos(x).series(x, 0, 8) == \
        pi/2 - x - x**3/6 - 3*x**5/40 - 5*x**7/112 + O(x**8)
    assert acos(x).series(x, 0, 8) == pi/2 - asin(x).series(x, 0, 8)
    t5 = acos(x).taylor_term(5, x)
    assert t5 == -3*x**5/40
    assert acos(x).taylor_term(7, x, t5, 0) == -5*x**7/112

def test_quaternion_conversions():
    q1 = Quaternion(1, 2, 3, 4)

    assert q1.to_axis_angle() == ((2 * sqrt(29)/29,
                                   3 * sqrt(29)/29,
                                   4 * sqrt(29)/29),
                                   2 * acos(sqrt(30)/30))

    assert q1.to_rotation_matrix() == Matrix([[-S(2)/3, S(2)/15, S(11)/15],
                                     [S(2)/3, -S(1)/3, S(14)/15],
                                     [S(1)/3, S(14)/15, S(2)/15]])

    assert q1.to_rotation_matrix((1, 1, 1)) == Matrix([[-S(2)/3, S(2)/15, S(11)/15, S(4)/5],
                                                  [S(2)/3, -S(1)/3, S(14)/15, -S(4)/15],
                                                  [S(1)/3, S(14)/15, S(2)/15, -S(2)/5],
                                                  [S(0), S(0), S(0), S(1)]])

    theta = symbols("theta", real=True)
    q2 = Quaternion(cos(theta/2), 0, 0, sin(theta/2))

    assert trigsimp(q2.to_rotation_matrix()) == Matrix([
                                               [cos(theta), -sin(theta), 0],
                                               [sin(theta),  cos(theta), 0],
                                               [0,           0,          1]])

    assert q2.to_axis_angle() == ((0, 0, sin(theta/2)/Abs(sin(theta/2))),
                                   2*acos(cos(theta/2)))

    assert trigsimp(q2.to_rotation_matrix((1, 1, 1))) == Matrix([
               [cos(theta), -sin(theta), 0, sin(theta) - cos(theta) + 1],
               [sin(theta),  cos(theta), 0, -sin(theta) - cos(theta) + 1],
               [0,           0,          1,  0],
               [0,           0,          0,  1]])

def test_ode_solutions():
    # only a few examples here, the rest will be tested in the actual dsolve tests
    assert constant_renumber(constantsimp(C1*exp(2*x) + exp(x)*(C2 + C3), [C1, C2, C3])) == \
        constant_renumber((C1*exp(x) + C2*exp(2*x)))
    assert constant_renumber(
        constantsimp(Eq(f(x), I*C1*sinh(x/3) + C2*cosh(x/3)), [C1, C2])
        ) == constant_renumber(Eq(f(x), C1*sinh(x/3) + C2*cosh(x/3)))
    assert constant_renumber(constantsimp(Eq(f(x), acos((-C1)/cos(x))), [C1])) == \
        Eq(f(x), acos(C1/cos(x)))
    assert constant_renumber(
        constantsimp(Eq(log(f(x)/C1) + 2*exp(x/f(x)), 0), [C1])
        ) == Eq(log(C1*f(x)) + 2*exp(x/f(x)), 0)
    assert constant_renumber(constantsimp(Eq(log(x*sqrt(2)*sqrt(1/x)*sqrt(f(x))
        /C1) + x**2/(2*f(x)**2), 0), [C1])) == \
        Eq(log(C1*sqrt(x)*sqrt(f(x))) + x**2/(2*f(x)**2), 0)
    assert constant_renumber(constantsimp(Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(x/C1) -
        cos(f(x)/x)*exp(-f(x)/x)/2, 0), [C1])) == \
        Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(C1*x) - cos(f(x)/x)*
           exp(-f(x)/x)/2, 0)
    assert constant_renumber(constantsimp(Eq(-Integral(-1/(sqrt(1 - u2**2)*u2),
        (u2, _a, x/f(x))) + log(f(x)/C1), 0), [C1])) == \
        Eq(-Integral(-1/(u2*sqrt(1 - u2**2)), (u2, _a, x/f(x))) +
        log(C1*f(x)), 0)
    assert [constantsimp(i, [C1]) for i in [Eq(f(x), sqrt(-C1*x + x**2)), Eq(f(x), -sqrt(-C1*x + x**2))]] == \
        [Eq(f(x), sqrt(x*(C1 + x))), Eq(f(x), -sqrt(x*(C1 + x)))]

def test_ode_solutions():
    # only a few examples here, the rest will be tested in the actual dsolve tests
    assert constant_renumber(constantsimp(C1*exp(2*x)+exp(x)*(C2+C3), x, 3), 'C', 1, 3) == \
        constant_renumber(C1*exp(x)+C2*exp(2*x), 'C', 1, 2)
    assert constant_renumber(constantsimp(Eq(f(x),I*C1*sinh(x/3) + C2*cosh(x/3)), x, 2),
        'C', 1, 2) == constant_renumber(Eq(f(x), C1*sinh(x/3) + C2*cosh(x/3)), 'C', 1, 2)
    assert constant_renumber(constantsimp(Eq(f(x),acos((-C1)/cos(x))), x, 1), 'C', 1, 1) == \
        Eq(f(x),acos(C1/cos(x)))
    assert constant_renumber(constantsimp(Eq(log(f(x)/C1) + 2*exp(x/f(x)), 0), x, 1),
        'C', 1, 1) ==  Eq(log(C1*f(x)) + 2*exp(x/f(x)), 0)
    assert constant_renumber(constantsimp(Eq(log(x*sqrt(2)*sqrt(1/x)*sqrt(f(x))\
        /C1) + x**2/(2*f(x)**2), 0), x, 1), 'C', 1, 1) == \
        Eq(log(C1*x*sqrt(1/x)*sqrt(f(x))) + x**2/(2*f(x)**2), 0)
    assert constant_renumber(constantsimp(Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(x/C1) - \
        cos(f(x)/x)*exp(-f(x)/x)/2, 0), x, 1), 'C', 1, 1) == \
        Eq(-exp(-f(x)/x)*sin(f(x)/x)/2 + log(C1*x) - cos(f(x)/x)*exp(-f(x)/x)/2, 0)
    u2 = Symbol('u2')
    _a = Symbol('_a')
    assert constant_renumber(constantsimp(Eq(-Integral(-1/(sqrt(1 - u2**2)*u2), \
        (u2, _a, x/f(x))) + log(f(x)/C1), 0), x, 1), 'C', 1, 1) == \
        Eq(-Integral(-1/(u2*sqrt(1 - u2**2)), (u2, _a, x/f(x))) + \
        log(C1*f(x)), 0)
    assert [constant_renumber(constantsimp(i, x, 1), 'C', 1, 1) for i in
        [Eq(f(x), sqrt(-C1*x + x**2)), Eq(f(x), -sqrt(-C1*x +
        x**2))]] == [Eq(f(x), sqrt(C1*x + x**2)),
        Eq(f(x), -sqrt(C1*x + x**2))]

def test_issue_2850():
    assert manualintegrate(asin(x)*log(x), x) == -x*asin(x) - sqrt(-x**2 + 1) \
            + (x*asin(x) + sqrt(-x**2 + 1))*log(x) - Integral(sqrt(-x**2 + 1)/x, x)
    assert manualintegrate(acos(x)*log(x), x) == -x*acos(x) + sqrt(-x**2 + 1) + \
        (x*acos(x) - sqrt(-x**2 + 1))*log(x) + Integral(sqrt(-x**2 + 1)/x, x)
    assert manualintegrate(atan(x)*log(x), x) == -x*atan(x) + (x*atan(x) - \
            log(x**2 + 1)/2)*log(x) + log(x**2 + 1)/2 + Integral(log(x**2 + 1)/x, x)/2

def test_logcombine_1():
    x, y = symbols("x,y")
    a = Symbol("a")
    z, w = symbols("z,w", positive=True)
    b = Symbol("b", real=True)
    assert logcombine(log(x)+2*log(y)) == log(x) + 2*log(y)
    assert logcombine(log(x)+2*log(y), force=True) == log(x*y**2)
    assert logcombine(a*log(w)+log(z)) == a*log(w) + log(z)
    assert logcombine(b*log(z)+b*log(x)) == log(z**b) + b*log(x)
    assert logcombine(b*log(z)-log(w)) == log(z**b/w)
    assert logcombine(log(x)*log(z)) == log(x)*log(z)
    assert logcombine(log(w)*log(x)) == log(w)*log(x)
    assert logcombine(cos(-2*log(z)+b*log(w))) in [cos(log(w**b/z**2)),
                                                   cos(log(z**2/w**b))]
    assert logcombine(log(log(x)-log(y))-log(z), force=True) == \
        log(log((x/y)**(1/z)))
    assert logcombine((2+I)*log(x), force=True) == I*log(x)+log(x**2)
    assert logcombine((x**2+log(x)-log(y))/(x*y), force=True) == \
        log(x**(1/(x*y))*y**(-1/(x*y)))+x/y
    assert logcombine(log(x)*2*log(y)+log(z), force=True) == \
        log(z*y**log(x**2))
    assert logcombine((x*y+sqrt(x**4+y**4)+log(x)-log(y))/(pi*x**Rational(2, 3)*\
        sqrt(y)**3), force=True) == \
        log(x**(1/(pi*x**Rational(2, 3)*sqrt(y)**3))*y**(-1/(pi*\
        x**Rational(2, 3)*sqrt(y)**3))) + sqrt(x**4 + y**4)/(pi*\
        x**Rational(2, 3)*sqrt(y)**3) + x**Rational(1, 3)/(pi*sqrt(y))
    assert logcombine(Eq(log(x), -2*log(y)), force=True) == \
        Eq(log(x*y**2), Integer(0))
    assert logcombine(Eq(y, x*acos(-log(x/y))), force=True) == \
        Eq(y, x*acos(log(y/x)))
    assert logcombine(gamma(-log(x/y))*acos(-log(x/y)), force=True) == \
        acos(log(y/x))*gamma(log(y/x))
    assert logcombine((2+3*I)*log(x), force=True) == \
        log(x**2)+3*I*log(x)
    assert logcombine(Eq(y, -log(x)), force=True) == Eq(y, log(1/x))

def test_heurisch_trigonometric():
    assert heurisch(sin(x), x) == -cos(x)
    assert heurisch(pi*sin(x) + 1, x) == x - pi*cos(x)

    assert heurisch(cos(x), x) == sin(x)
    assert heurisch(tan(x), x) in [
        log(1 + tan(x)**2)/2,
        log(tan(x) + I) + I*x,
        log(tan(x) - I) - I*x,
    ]

    assert heurisch(sin(x)*sin(y), x) == -cos(x)*sin(y)
    assert heurisch(sin(x)*sin(y), y) == -cos(y)*sin(x)

    # gives sin(x) in answer when run via setup.py and cos(x) when run via py.test
    assert heurisch(sin(x)*cos(x), x) in [sin(x)**2 / 2, -cos(x)**2 / 2]
    assert heurisch(cos(x)/sin(x), x) == log(sin(x))

    assert heurisch(x*sin(7*x), x) == sin(7*x) / 49 - x*cos(7*x) / 7
    assert heurisch(1/pi/4 * x**2*cos(x), x) == 1/pi/4*(x**2*sin(x) -
                    2*sin(x) + 2*x*cos(x))

    assert heurisch(acos(x/4) * asin(x/4), x) == 2*x - (sqrt(16 - x**2))*asin(x/4) \
        + (sqrt(16 - x**2))*acos(x/4) + x*asin(x/4)*acos(x/4)

    assert heurisch(sin(x)/(cos(x)**2+1), x) == -atan(cos(x)) #fixes issue 13723
    assert heurisch(1/(cos(x)+2), x) == 2*sqrt(3)*atan(sqrt(3)*tan(x/2)/3)/3
    assert heurisch(2*sin(x)*cos(x)/(sin(x)**4 + 1), x) == atan(sqrt(2)*sin(x)
        - 1) - atan(sqrt(2)*sin(x) + 1)

    assert heurisch(1/cosh(x), x) == 2*atan(tanh(x/2))

    def eval(cls, arg):
        from sympy import acos
        arg = sympify(arg)

        if arg.is_Number:
            if arg is S.NaN:
                return S.NaN
            elif arg is S.Infinity:
                return S.Infinity
            elif arg is S.NegativeInfinity:
                return S.Infinity
            elif arg is S.Zero:
                return S.Pi*S.ImaginaryUnit / 2
            elif arg is S.One:
                return S.Zero
            elif arg is S.NegativeOne:
                return S.Pi*S.ImaginaryUnit

        if arg.is_number:
            cst_table = {
                S.ImaginaryUnit: log(S.ImaginaryUnit*(1 + sqrt(2))),
                -S.ImaginaryUnit: log(-S.ImaginaryUnit*(1 + sqrt(2))),
                S.Half: S.Pi/3,
                -S.Half: 2*S.Pi/3,
                sqrt(2)/2: S.Pi/4,
                -sqrt(2)/2: 3*S.Pi/4,
                1/sqrt(2): S.Pi/4,
                -1/sqrt(2): 3*S.Pi/4,
                sqrt(3)/2: S.Pi/6,
                -sqrt(3)/2: 5*S.Pi/6,
                (sqrt(3) - 1)/sqrt(2**3): 5*S.Pi/12,
                -(sqrt(3) - 1)/sqrt(2**3): 7*S.Pi/12,
                sqrt(2 + sqrt(2))/2: S.Pi/8,
                -sqrt(2 + sqrt(2))/2: 7*S.Pi/8,
                sqrt(2 - sqrt(2))/2: 3*S.Pi/8,
                -sqrt(2 - sqrt(2))/2: 5*S.Pi/8,
                (1 + sqrt(3))/(2*sqrt(2)): S.Pi/12,
                -(1 + sqrt(3))/(2*sqrt(2)): 11*S.Pi/12,
                (sqrt(5) + 1)/4: S.Pi/5,
                -(sqrt(5) + 1)/4: 4*S.Pi/5
            }

            if arg in cst_table:
                if arg.is_real:
                    return cst_table[arg]*S.ImaginaryUnit
                return cst_table[arg]

        if arg is S.ComplexInfinity:
            return S.Infinity

        i_coeff = arg.as_coefficient(S.ImaginaryUnit)

        if i_coeff is not None:
            if _coeff_isneg(i_coeff):
                return S.ImaginaryUnit * acos(i_coeff)
            return S.ImaginaryUnit * acos(-i_coeff)
        else:
            if _coeff_isneg(arg):
                return -cls(-arg)

def test_asech():
    x = Symbol('x')

    assert asech(-x) == asech(-x)

    # values at fixed points
    assert asech(1) == 0
    assert asech(-1) == pi*I
    assert asech(0) == oo
    assert asech(2) == I*pi/3
    assert asech(-2) == 2*I*pi / 3

    # at infinites
    assert asech(oo) == I*pi/2
    assert asech(-oo) == I*pi/2
    assert asech(zoo) == nan

    assert asech(I) == log(1 + sqrt(2)) - I*pi/2
    assert asech(-I) == log(1 + sqrt(2)) + I*pi/2
    assert asech(sqrt(2) - sqrt(6)) == 11*I*pi / 12
    assert asech(sqrt(2 - 2/sqrt(5))) == I*pi / 10
    assert asech(-sqrt(2 - 2/sqrt(5))) == 9*I*pi / 10
    assert asech(2 / sqrt(2 + sqrt(2))) == I*pi / 8
    assert asech(-2 / sqrt(2 + sqrt(2))) == 7*I*pi / 8
    assert asech(sqrt(5) - 1) == I*pi / 5
    assert asech(1 - sqrt(5)) == 4*I*pi / 5
    assert asech(-sqrt(2*(2 + sqrt(2)))) == 5*I*pi / 8

    # properties
    # asech(x) == acosh(1/x)
    assert asech(sqrt(2)) == acosh(1/sqrt(2))
    assert asech(2/sqrt(3)) == acosh(sqrt(3)/2)
    assert asech(2/sqrt(2 + sqrt(2))) == acosh(sqrt(2 + sqrt(2))/2)
    assert asech(S(2)) == acosh(1/S(2))

    # asech(x) == I*acos(1/x)
    # (Note: the exact formula is asech(x) == +/- I*acos(1/x))
    assert asech(-sqrt(2)) == I*acos(-1/sqrt(2))
    assert asech(-2/sqrt(3)) == I*acos(-sqrt(3)/2)
    assert asech(-S(2)) == I*acos(-S.Half)
    assert asech(-2/sqrt(2)) == I*acos(-sqrt(2)/2)

    # sech(asech(x)) / x == 1
    assert expand_mul(sech(asech(sqrt(6) - sqrt(2))) / (sqrt(6) - sqrt(2))) == 1
    assert expand_mul(sech(asech(sqrt(6) + sqrt(2))) / (sqrt(6) + sqrt(2))) == 1
    assert (sech(asech(sqrt(2 + 2/sqrt(5)))) / (sqrt(2 + 2/sqrt(5)))).simplify() == 1
    assert (sech(asech(-sqrt(2 + 2/sqrt(5)))) / (-sqrt(2 + 2/sqrt(5)))).simplify() == 1
    assert (sech(asech(sqrt(2*(2 + sqrt(2))))) / (sqrt(2*(2 + sqrt(2))))).simplify() == 1
    assert expand_mul(sech(asech((1 + sqrt(5)))) / ((1 + sqrt(5)))) == 1
    assert expand_mul(sech(asech((-1 - sqrt(5)))) / ((-1 - sqrt(5)))) == 1
    assert expand_mul(sech(asech((-sqrt(6) - sqrt(2)))) / ((-sqrt(6) - sqrt(2)))) == 1

    # numerical evaluation
    assert str(asech(5*I).n(6)) == '0.19869 - 1.5708*I'
    assert str(asech(-5*I).n(6)) == '0.19869 + 1.5708*I'

    def test_main_trig_functions_numeric(self):
        print "\n\n\n" + " Test if sin, cos and tan and inverses Work Numerically ".center(75, "#")
        from sympy import symbols, sin, cos, tan, asin, acos, atan
        x, y = symbols('x,y')
        test_expr = sin(x) + cos(x) + tan(x) + asin(y) + acos(y) + atan(y)
        target_expr = sin(x) + cos(x) + tan(x) + asin(y) + acos(y) + atan(y)

        print "Target expression: '%s'" % target_expr
        print "Test expression: '%s'" % test_expr
        equal = api.numeric_equality(test_expr, target_expr)

        self.assertTrue(equal, "Expected expressions to be found numerically equal!")
        print "   PASS   ".center(75, "#")

def test_acos_rewrite():
    assert acos(x).rewrite(log) == pi/2 + I*log(I*x + sqrt(1 - x**2))
    assert acos(x).rewrite(atan) == \
           atan(sqrt(1 - x**2)/x) + (pi/2)*(1 - x*sqrt(1/x**2))
    assert acos(0).rewrite(atan) == S.Pi/2
    assert acos(0.5).rewrite(atan) == acos(0.5).rewrite(log)
    assert acos(x).rewrite(asin) == S.Pi/2 - asin(x)
    assert acos(x).rewrite(acot) == -2*acot((sqrt(-x**2 + 1) + 1)/x) + pi/2
    assert acos(x).rewrite(asec) == asec(1/x)
    assert acos(x).rewrite(acsc) == -acsc(1/x) + pi/2

def get_rec_cell_parameters(a, b, c, alpha, beta, gamma):
    import sympy as sp
    """
        Returns cell vectors in direkt and reciprocal space and reciprocal
        lattice parameters of given cell parameters (vectors in direct crystal
        fixed carthesian system).
    """
    
    # Metric Tensor:
    G = sp.Matrix([[a**2,               a*b*sp.cos(gamma), a*c*sp.cos(beta)],
                   [a*b*sp.cos(gamma), b**2,              b*c*sp.cos(alpha)],
                   [a*c*sp.cos(beta),  b*c*sp.cos(alpha), c**2]])
    
    G_r = G.inv() # reciprocal Metric
    G_r.simplify()
    
    # volume of crystall system cell in carthesian system
    # V = sp.sqrt(G.det())
    
    
    # reciprocal cell lengths
    ar = sp.sqrt(G_r[0,0])
    br = sp.sqrt(G_r[1,1])
    cr = sp.sqrt(G_r[2,2])
    
    #alphar = sp.acos(G_r[1,2]/(br*cr)).simplify()
    alphar = sp.acos((-sp.cos(alpha) + sp.cos(beta)*sp.cos(gamma))/(sp.Abs(sp.sin(beta))*sp.Abs(sp.sin(gamma))))
    betar  = sp.acos(G_r[0,2]/(ar*cr))
    #gammar = sp.acos(G_r[0,1]/(ar*br)).simplify()
    gammar = sp.acos((sp.cos(alpha)*sp.cos(beta) - sp.cos(gamma))/(sp.Abs(sp.sin(alpha))*sp.Abs(sp.sin(beta))))
    
    # x parallel to a* and z parallel to a* x b*   (ITC Vol B Ch. 3.3.1.1.1)
    #B =   sp.Matrix([[ar, br*sp.cos(gammar),  cr*sp.cos(betar)],
    #                 [0,  br*sp.sin(gammar), -cr*sp.sin(betar)*sp.cos(alpha)],
    #                 [0,  0,                  1/c]])
    #B_0 = sp.Matrix([[1,  sp.cos(gammar),  sp.cos(betar)],
    #                 [0,  sp.sin(gammar), -sp.sin(betar)*sp.cos(alpha)],
    #                 [0,  0,               1]])
    #V_0 = sp.sqrt(1 - sp.cos(alphar)**2 - sp.cos(betar)**2 - sp.cos(gammar)**2 + 2*sp.cos(alphar)*sp.cos(betar)*sp.cos(gammar))
    
    # x parallel to a and z parallel to a x b   (ITC Vol B Ch. 3.3.1.1.1)
    V0 = sp.sin(alpha) * sp.sin(beta) * sp.sin(gammar)
    M = sp.Matrix([[a, b*sp.cos(gamma),  c*sp.cos(beta)],
                   [0, b*sp.sin(gamma),  c*(sp.cos(alpha) - sp.cos(beta)*sp.cos(gamma))/sp.sin(gamma)],
                   [0,               0,  c*V0/sp.sin(gamma)]])
    Mi = sp.Matrix([[1/a,-1/(a*sp.tan(gamma)),  (sp.cos(alpha)*sp.cos(gamma) - sp.cos(beta))/(a*V0*sp.sin(gamma))],
                    [0,   1/(b*sp.sin(gamma)),  (sp.cos(beta)*sp.cos(gamma) - sp.cos(alpha))/(b*V0*sp.sin(gamma))],
                    [0,   0,                     sp.sin(gamma)/(c*V0)]])
    
    return (ar, br, cr, alphar, betar, gammar, M, Mi, G, G_r)

def test_function__eval_nseries():
    n = Symbol('n')

    assert sin(x)._eval_nseries(x, 2, None) == x + O(x**2)
    assert sin(x + 1)._eval_nseries(x, 2, None) == x*cos(1) + sin(1) + O(x**2)
    assert sin(pi*(1 - x))._eval_nseries(x, 2, None) == pi*x + O(x**2)
    assert acos(1 - x**2)._eval_nseries(x, 2, None) == sqrt(2)*x + O(x**2)
    assert polygamma(n, x + 1)._eval_nseries(x, 2, None) == \
        polygamma(n, 1) + polygamma(n + 1, 1)*x + O(x**2)
    raises(PoleError, lambda: sin(1/x)._eval_nseries(x, 2, None))
    raises(PoleError, lambda: acos(1 - x)._eval_nseries(x, 2, None))
    raises(PoleError, lambda: acos(1 + x)._eval_nseries(x, 2, None))
    assert loggamma(1/x)._eval_nseries(x, 0, None) == \
        log(x)/2 - log(x)/x - 1/x + O(1, x)
    assert loggamma(log(1/x)).nseries(x, n=1, logx=y) == loggamma(-y)

def test_roots_cubic():
    assert roots_cubic(Poly(2*x**3, x)) == [0, 0, 0]
    assert roots_cubic(Poly(x**3 - 3*x**2 + 3*x - 1, x)) == [1, 1, 1]

    assert roots_cubic(Poly(x**3 + 1, x)) == \
        [-1, S.Half - I*sqrt(3)/2, S.Half + I*sqrt(3)/2]
    assert roots_cubic(Poly(2*x**3 - 3*x**2 - 3*x - 1, x))[0] == \
         S.Half + 3**Rational(1, 3)/2 + 3**Rational(2, 3)/2
    eq = -x**3 + 2*x**2 + 3*x - 2
    assert roots(eq, trig=True, multiple=True) == \
           roots_cubic(Poly(eq, x), trig=True) == [
        S(2)/3 + 2*sqrt(13)*cos(acos(8*sqrt(13)/169)/3)/3,
        -2*sqrt(13)*sin(-acos(8*sqrt(13)/169)/3 + pi/6)/3 + S(2)/3,
        -2*sqrt(13)*cos(-acos(8*sqrt(13)/169)/3 + pi/3)/3 + S(2)/3,
        ]

def log(T):
    (R,p) = breakTrans(T)
    if np.array_equiv(R, np.eye(3)):
        w = np.zeros((3,1))
        theta = 1
        u = p
        
    else:
        tr = R.trace()
        if tr == -1:
            theta = -math.pi
            w = np.array([[ sqrt((R[i][i]+1)/2) for i in range(3)]]).T
            W = bra(w)
        else:
            theta = acos((tr-1)/2)
            W = (1/(2*sin(theta)))*(R-R.T)
            w = ibra(W)
        
        W2 = np.dot(W,W)
        G1 = 1./theta*np.eye(3)
        G2 = 0.5*W
        G3 = (1./theta-cot(theta/2)/2)*W2
        G = G1+G2+G3
        # This corrects the bug!!!!!
        G = G.T
        u = np.dot(G,p)
    
    h = np.dot(w.T,u)
    return (w,u,theta,h)

def test_intrinsic_math1_codegen():
    # not included: log10
    from sympy import acos, asin, atan, ceiling, cos, cosh, floor, log, ln, \
        sin, sinh, sqrt, tan, tanh, N
    x = symbols('x')
    name_expr = [
        ("test_fabs", abs(x)),
        ("test_acos", acos(x)),
        ("test_asin", asin(x)),
        ("test_atan", atan(x)),
        ("test_cos", cos(x)),
        ("test_cosh", cosh(x)),
        ("test_log", log(x)),
        ("test_ln", ln(x)),
        ("test_sin", sin(x)),
        ("test_sinh", sinh(x)),
        ("test_sqrt", sqrt(x)),
        ("test_tan", tan(x)),
        ("test_tanh", tanh(x)),
    ]
    numerical_tests = []
    for name, expr in name_expr:
        for xval in 0.2, 0.5, 0.8:
            expected = N(expr.subs(x, xval))
            numerical_tests.append((name, (xval,), expected, 1e-14))
    for lang, commands in valid_lang_commands:
        if lang == "C":
            name_expr_C = [("test_floor", floor(x)), ("test_ceil", ceiling(x))]
        else:
            name_expr_C = []
        run_test("intrinsic_math1", name_expr + name_expr_C, numerical_tests, lang, commands)

def angle_between_vectors(a, b):
    """Return the minimum angle between two vectors. The angle returned for
    vectors a and -a is 0.
    """
    angle = (acos(dot(a, b)/(a.magnitude() * b.magnitude())) *
             180 / pi).evalf()
    return min(angle, 180 - angle)

    def _set_inv_trans_equations(curv_coord_name):
        """
        Store information about inverse transformation equations for
        pre-defined coordinate systems.

        Parameters
        ==========

        curv_coord_name : str
            Name of coordinate system

        """
        if curv_coord_name == 'cartesian':
            return lambda x, y, z: (x, y, z)

        if curv_coord_name == 'spherical':
            return lambda x, y, z: (
                sqrt(x**2 + y**2 + z**2),
                acos(z/sqrt(x**2 + y**2 + z**2)),
                atan2(y, x)
            )
        if curv_coord_name == 'cylindrical':
            return lambda x, y, z: (
                sqrt(x**2 + y**2),
                atan2(y, x),
                z
            )
        raise ValueError('Wrong set of parameters.'
                         'Type of coordinate system is defined')

def test_atan_rewrite():
    assert atan(x).rewrite(log) == I*log((1 - I*x)/(1 + I*x))/2
    assert atan(x).rewrite(asin) == (-asin(1/sqrt(x**2 + 1)) + pi/2)*sqrt(x**2)/x
    assert atan(x).rewrite(acos) == sqrt(x**2)*acos(1/sqrt(x**2 + 1))/x
    assert atan(x).rewrite(acot) == acot(1/x)
    assert atan(x).rewrite(asec) == sqrt(x**2)*asec(sqrt(x**2 + 1))/x
    assert atan(x).rewrite(acsc) == (-acsc(sqrt(x**2 + 1)) + pi/2)*sqrt(x**2)/x