def test_functional_diffgeom_ch4():
    x0, y0, theta0 = symbols('x0, y0, theta0', real=True)
    x, y, r, theta = symbols('x, y, r, theta', real=True)
    r0 = symbols('r0', positive=True)
    f = Function('f')
    b1 = Function('b1')
    b2 = Function('b2')
    p_r = R2_r.point([x0, y0])
    p_p = R2_p.point([r0, theta0])

    f_field = b1(R2.x,R2.y)*R2.dx + b2(R2.x,R2.y)*R2.dy
    assert f_field(R2.e_x)(p_r) == b1(x0, y0)
    assert f_field(R2.e_y)(p_r) == b2(x0, y0)

    s_field_r = f(R2.x,R2.y)
    df = Differential(s_field_r)
    assert df(R2.e_x)(p_r).doit() == Derivative(f(x0, y0), x0)
    assert df(R2.e_y)(p_r).doit() == Derivative(f(x0, y0), y0)

    s_field_p = f(R2.r,R2.theta)
    df = Differential(s_field_p)
    assert trigsimp(df(R2.e_x)(p_p).doit()) == cos(theta0)*Derivative(f(r0, theta0), r0) - sin(theta0)*Derivative(f(r0, theta0), theta0)/r0
    assert trigsimp(df(R2.e_y)(p_p).doit()) == sin(theta0)*Derivative(f(r0, theta0), r0) + cos(theta0)*Derivative(f(r0, theta0), theta0)/r0

    assert R2.dx(R2.e_x)(p_r) == 1
    assert R2.dx(R2.e_y)(p_r) == 0

    circ = -R2.y*R2.e_x + R2.x*R2.e_y
    assert R2.dx(circ)(p_r).doit() == -y0
    assert R2.dy(circ)(p_r) == x0
    assert R2.dr(circ)(p_r) == 0
    assert simplify(R2.dtheta(circ)(p_r)) == 1

    assert (circ - R2.e_theta)(s_field_r)(p_r) == 0

def test_vector_simplify():
    A, s, k, m = symbols("A, s, k, m")

    test1 = (1 / a + 1 / b) * i
    assert (test1 & i) != (a + b) / (a * b)
    test1 = simplify(test1)
    assert (test1 & i) == (a + b) / (a * b)
    assert test1.simplify() == simplify(test1)

    test2 = (A ** 2 * s ** 4 / (4 * pi * k * m ** 3)) * i
    test2 = simplify(test2)
    assert (test2 & i) == (A ** 2 * s ** 4 / (4 * pi * k * m ** 3))

    test3 = ((4 + 4 * a - 2 * (2 + 2 * a)) / (2 + 2 * a)) * i
    test3 = simplify(test3)
    assert (test3 & i) == 0

    test4 = ((-4 * a * b ** 2 - 2 * b ** 3 - 2 * a ** 2 * b) / (a + b) ** 2) * i
    test4 = simplify(test4)
    assert (test4 & i) == -2 * b

    v = (sin(a) + cos(a)) ** 2 * i - j
    assert trigsimp(v) == (2 * sin(a + pi / 4) ** 2) * i + (-1) * j
    assert trigsimp(v) == v.trigsimp()

    assert simplify(Vector.zero) == Vector.zero

def test_R2():
    x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0', real=True)
    point_r = R2_r.point([x0, y0])
    point_p = R2_p.point([r0, theta0])

    # r**2 = x**2 + y**2
    assert (R2.r**2 - R2.x**2 - R2.y**2)(point_r) == 0
    assert trigsimp( (R2.r**2 - R2.x**2 - R2.y**2)(point_p) ) == 0
    assert trigsimp(R2.e_r(R2.x**2+R2.y**2)(point_p).doit()) == 2*r0

    # polar->rect->polar == Id
    a, b = symbols('a b', positive=True)
    m = Matrix([[a], [b]])
    #TODO assert m == R2_r.coord_tuple_transform_to(R2_p, R2_p.coord_tuple_transform_to(R2_r, [a, b])).applyfunc(simplify)
    assert m == R2_p.coord_tuple_transform_to(R2_r, R2_r.coord_tuple_transform_to(R2_p, m)).applyfunc(simplify)

    def express_coordinates(self, coordinate_system):
        """
        Returns the Cartesian/rectangular coordinates of this point
        wrt the origin of the given CoordSysCartesian instance.

        Parameters
        ==========

        coordinate_system : CoordSysCartesian
            The coordinate system to express the coordinates of this
            Point in.

        Examples
        ========

        >>> from sympy.vector import Point, CoordSysCartesian
        >>> N = CoordSysCartesian('N')
        >>> p1 = N.origin.locate_new('p1', 10 * N.i)
        >>> p2 = p1.locate_new('p2', 5 * N.j)
        >>> p2.express_coordinates(N)
        (10, 5, 0)

        """

        #Determine the position vector
        pos_vect = self.position_wrt(coordinate_system.origin)
        #Express it in the given coordinate system
        pos_vect = trigsimp(express(pos_vect, coordinate_system,
                                    variables = True))
        coords = []
        for vect in coordinate_system.base_vectors():
            coords.append(pos_vect.dot(vect))
        return tuple(coords)

 def _eval_rewrite_as_cos(self, n, m, theta, phi):
     # This method can be expensive due to extensive use of simplification!
     from sympy.simplify import simplify, trigsimp
     # TODO: Make sure n \in N
     # TODO: Assert |m| <= n ortherwise we should return 0
     term = simplify(self.expand(func=True))
     # We can do this because of the range of theta
     term = term.xreplace({Abs(sin(theta)):sin(theta)})
     return simplify(trigsimp(term))

    def __contains__(self, o):
        if isinstance(o, Point):
            x = Dummy('x', real=True)
            y = Dummy('y', real=True)

            res = self.equation(x, y).subs({x: o.x, y: o.y})
            return trigsimp(simplify(res)) is S.Zero
        elif isinstance(o, Ellipse):
            return self == o
        return False

    def __contains__(self, o):
        if isinstance(o, Point):
            x = C.Dummy("x", real=True)
            y = C.Dummy("y", real=True)

            res = self.equation(x, y).subs({x: o[0], y: o[1]})
            return trigsimp(simplify(res)) is S.Zero
        elif isinstance(o, Ellipse):
            return self == o
        return False

 def __contains__(self, o):
     if isinstance(o, Point):
         x = Basic.Symbol('x', real=True)
         y = Basic.Symbol('y', real=True)
         res = self.equation('x', 'y').subs_dict({x: o[0], y: o[1]})
         res = trigsimp(simplify(res))
         return res == 0
     elif isinstance(o, Ellipse):
         return (self == o)
     return False

    def __contains__(self, o):
        if isinstance(o, Point):
            x = C.Symbol("x", real=True, dummy=True)
            y = C.Symbol("y", real=True, dummy=True)

            res = self.equation(x, y).subs({x: o[0], y: o[1]})
            res = trigsimp(simplify(res))
            return res == 0
        elif isinstance(o, Ellipse):
            return self == o
        return False

    def __contains__(self, o):
        if isinstance(o, Point):
            x = C.Dummy('x', real=True)
            y = C.Dummy('y', real=True)

            res = self.equation(x, y).subs({x: o[0], y: o[1]})
            res = trigsimp(simplify(res))
            return res == 0
        elif isinstance(o, Ellipse):
            return (self == o)
        return False

 def trigsimp(self, deep=False, recursive=False):
     """See the trigsimp function in sympy.simplify"""
     from sympy.simplify import trigsimp
     return trigsimp(self, deep, recursive)

 def _eval_trigsimp(self, **opts):
     from sympy.simplify import trigsimp
     return self.func(trigsimp(self.lhs, **opts), trigsimp(self.rhs, **opts))