def arbitrary_point(self, parameter_name='t'):
        """A parametrised point on the Line.

        Parameters
        ----------
        parameter_name : str, optional
            The name of the parameter which will be used for the parametric
            point. The default value is 't'.

        Returns
        -------
        point : Point

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy import Point, Line
        >>> p1, p2 = Point(1, 0), Point(5, 3)
        >>> l1 = Line(p1, p2)
        >>> l1.arbitrary_point()
        Point(1 + 4*t, 3*t)

        """
        t = C.Symbol(parameter_name, real=True)
        x = simplify(self.p1[0] + t*(self.p2[0] - self.p1[0]))
        y = simplify(self.p1[1] + t*(self.p2[1] - self.p1[1]))
        return Point(x, y)

    def centroid(self):
        """The centroid of the polygon.

        Returns
        -------
        centroid : Point

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy import Point, Polygon
        >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)])
        >>> poly = Polygon(p1, p2, p3, p4)
        >>> poly.centroid
        Point(31/18, 11/18)

        """
        A = 1 / (6*self.area)
        cx,cy = 0,0
        for ind in xrange(-1, len(self.vertices)-1):
            pi = self.vertices[ind]
            pii = self.vertices[ind+1]
            v = pi[0]*pii[1]-pii[0]*pi[1]
            cx += v*(pi[0] + pii[0])
            cy += v*(pi[1] + pii[1])
        return Point(simplify(A*cx), simplify(A*cy))

    def centroid(self):
        """The centroid of the polygon.

        Returns
        -------
        centroid : Point

        See Also
        --------
        Point

        Examples
        --------
        >>> from sympy import Point, Polygon
        >>> p1, p2, p3, p4 = map(Point, [(0, 0), (1, 0), (5, 1), (0, 1)])
        >>> poly = Polygon(p1, p2, p3, p4)
        >>> poly.centroid
        Point(31/18, 11/18)

        """
        A = 1/(6*self.area)
        cx, cy = 0, 0
        for i in xrange(len(self)):
            pt1 = self[i - 1]
            pt2 = self[i]
            v = pt1[0]*pt2[1] - pt2[0]*pt1[1]
            cx += v*(pt1[0] + pt2[0])
            cy += v*(pt1[1] + pt2[1])
        return Point(simplify(A*cx), simplify(A*cy))

    def incenter(self):
        """The center of the incircle.

        The incircle is the circle which lies inside the triangle and touches
        all three sides.

        Returns
        -------
        incenter : Point

        See Also
        --------
        incircle
        Point

        Examples
        --------
        >>> from sympy.geometry import Point, Triangle
        >>> p1, p2, p3 = Point(0, 0), Point(1, 0), Point(0, 1)
        >>> t = Triangle(p1, p2, p3)
        >>> t.incenter
        Point(-sqrt(2)/2 + 1, -sqrt(2)/2 + 1)

        """
        s = self.sides
        v = self.vertices
        A,B,C = v[0],v[1],v[2]
        a,b,c = s[1].length,s[2].length,s[0].length
        x = simplify((a*A[0] + b*B[0] + c*C[0]) / (a+b+c))
        y = simplify((a*A[1] + b*B[1] + c*C[1]) / (a+b+c))
        return Point(x, y)

 def __new__(cls, radius=1, center=[0,0,0], direction=[0,0,1], closed=False, **kwargs):
     """
     >>> from sympy import *
     >>> from symplus.strplus import init_mprinting
     >>> init_mprinting()
     >>> InfiniteCylinder()
     InfiniteCylinder(1, [0 0 0]', [0 0 1]', False)
     >>> InfiniteCylinder(2, [0,0,0], [0,1,1])
     InfiniteCylinder(2, [0 0 0]', [0 -sqrt(2)/2 -sqrt(2)/2]', False)
     >>> InfiniteCylinder().contains((1,1,1))
     False
     >>> InfiniteCylinder(2, [0,0,0], [0,1,1]).contains((1,1,1))
     True
     """
     normalization = kwargs.pop("normalization", True)
     radius = sympify(abs(radius))
     direction = Mat(direction)
     if normalization:
         if norm(direction) == 0:
             raise ValueError
         direction = simplify(normalize(direction))
     direction = max(direction, -direction, key=hash)
     center = Mat(center)
     if normalization:
         center = simplify(center - project(center, direction))
     closed = sympify(bool(closed))
     return Basic.__new__(cls, radius, center, direction, closed)

def test_reciprocal_frame_test():
    with GA_Printer():
        metric = '1 # #,' + \
                 '# 1 #,' + \
                 '# # 1,'

        (e1, e2, e3) = MV.setup('e1 e2 e3', metric)

        E = e1 ^ e2 ^ e3
        Esq = (E*E).scalar()
        assert str(E) == 'e1^e2^e3'
        assert str(Esq) == '(e1.e2)**2 - 2*(e1.e2)*(e1.e3)*(e2.e3) + (e1.e3)**2 + (e2.e3)**2 - 1'
        Esq_inv = 1/Esq

        E1 = (e2 ^ e3)*E
        E2 = (-1)*(e1 ^ e3)*E
        E3 = (e1 ^ e2)*E

        assert str(E1) == '((e2.e3)**2 - 1)*e1 + ((e1.e2) - (e1.e3)*(e2.e3))*e2 + (-(e1.e2)*(e2.e3) + (e1.e3))*e3'
        assert str(E2) == '((e1.e2) - (e1.e3)*(e2.e3))*e1 + ((e1.e3)**2 - 1)*e2 + (-(e1.e2)*(e1.e3) + (e2.e3))*e3'
        assert str(E3) == '(-(e1.e2)*(e2.e3) + (e1.e3))*e1 + (-(e1.e2)*(e1.e3) + (e2.e3))*e2 + ((e1.e2)**2 - 1)*e3'

        w = (E1 | e2)
        w = w.expand()
        assert str(w) == '0'

        w = (E1 | e3)
        w = w.expand()
        assert str(w) == '0'

        w = (E2 | e1)
        w = w.expand()
        assert str(w) == '0'

        w = (E2 | e3)
        w = w.expand()
        assert str(w) == '0'

        w = (E3 | e1)
        w = w.expand()
        assert str(w) == '0'

        w = (E3 | e2)
        w = w.expand()
        assert str(w) == '0'

        w = (E1 | e1)
        w = (w.expand()).scalar()
        Esq = expand(Esq)
        assert str(simplify(w/Esq)) == '1'

        w = (E2 | e2)
        w = (w.expand()).scalar()
        assert str(simplify(w/Esq)) == '1'

        w = (E3 | e3)
        w = (w.expand()).scalar()
        assert str(simplify(w/Esq)) == '1'

    return

    def is_collinear(*points):
        """Is a sequence of points collinear?

        Test whether or not a set of points are collinear. Returns True if
        the set of points are collinear, or False otherwise.

        Parameters
        ----------
        points : sequence of Point

        Returns
        -------
        is_collinear : boolean

        Notes
        --------------------------
        Slope is preserved everywhere on a line, so the slope between
        any two points on the line should be the same. Take the first
        two points, p1 and p2, and create a translated point v1
        with p1 as the origin. Now for every other point we create
        a translated point, vi with p1 also as the origin. Note that
        these translations preserve slope since everything is
        consistently translated to a new origin of p1. Since slope
        is preserved then we have the following equality:
                v1_slope = vi_slope
          =>    v1.y/v1.x = vi.y/vi.x (due to translation)
          =>    v1.y*vi.x = vi.y*v1.x
          =>    v1.y*vi.x - vi.y*v1.x = 0           (*)
        Hence, if we have a vi such that the equality in (*) is False
        then the points are not collinear. We do this test for every
        point in the list, and if all pass then they are collinear.

        Examples
        --------
        >>> from sympy import Point
        >>> from sympy.abc import x
        >>> p1, p2 = Point(0, 0), Point(1, 1)
        >>> p3, p4, p5 = Point(2, 2), Point(x, x), Point(1, 2)
        >>> Point.is_collinear(p1, p2, p3, p4)
        True
        >>> Point.is_collinear(p1, p2, p3, p5)
        False

        """
        points = GeometryEntity.extract_entities(points)
        if len(points) == 0: return False
        if len(points) <= 2: return True # two points always form a line

        # XXX Cross product is used now, but that only extends to three
        #     dimensions. If the concept needs to extend to greater
        #     dimensions then another method would have to be used
        p1 = points[0]
        p2 = points[1]
        v1 = p2 - p1
        for p3 in points[2:]:
            v2 = p3 - p1
            test = simplify(v1[0]*v2[1] - v1[1]*v2[0])
            if simplify(test) != 0:
                return False
        return True

def custom_intersection(c1,c2):
    #check if they are the same circle
    if c1.center == c2.center:
        if c2.radius == c1.radius:
            return c2
        return []

    dx, dy = (c2.center - c1.center).args
    d = (dy**2 + dx**2)**(.5)
    a = (c1.radius**2 - c2.radius**2 + d**2) / (2*d)

    x2 = c1.center.x + (dx * a/d)
    y2 = c1.center.y + (dy * a/d)

    h = (c1.radius**2 - a**2)**(.5)
    rx = -dy * (h/d)
    ry = dx * (h/d)

    xi_1 = simplify(x2 + rx)
    xi_2 = simplify(x2 - rx)
    yi_1 = simplify(y2 + ry)
    yi_2 = simplify(y2 - ry)

    ret = [Point(xi_1, yi_1)]
    if xi_1 != xi_2 or yi_1 != yi_2:
        ret.append(Point(xi_2, yi_2))
    return ret

    def tangent_line(self, p):
        """
        If p is on the ellipse, returns the tangent line through point p.
        Otherwise, returns the tangent line(s) from p to the ellipse, or
        None if no tangent line is possible (e.g., p inside ellipse).

        Example:

        In [1]: e = Ellipse(Point(0,0), 3, 2)

        In [2]: t = e.tangent_line(e.random_point())

        In [3]: p = Plot()

        In [4]: p[0] = e

        In [5]: p[1] = t

        The above will plot an ellipse together with a tangent line.
        """
        if p in self:
            rise = (self.vradius ** 2)*(self.center[0] - p[0])
            run = (self.hradius ** 2)*(p[1] - self.center[1])
            p2 = Point(simplify(p[0] + run),
                       simplify(p[1] + rise))
            return Line(p, p2)
        else:
            # TODO If p is not on the ellipse, attempt to create the
            #      tangent(s) from point p to the ellipse..?
            raise NotImplementedError("Cannot find tangent lines when p is not on the ellipse")

def test_reciprocal_frame():
    """
    Test of formula for general reciprocal frame of three vectors.
    Let three independent vectors be e1, e2, and e3. The reciprocal
    vectors E1, E2, and E3 obey the relations:

    e_i.E_j = delta_ij*(e1^e2^e3)**2
    """
    g = '1 # #,'+ \
        '# 1 #,'+ \
        '# # 1'

    g3dn = Ga('e1 e2 e3',g=g)

    (e1,e2,e3) = g3dn.mv()

    E = e1^e2^e3
    Esq = (E*E).scalar()
    Esq_inv = 1 / Esq

    E1 = (e2^e3)*E
    E2 = (-1)*(e1^e3)*E
    E3 = (e1^e2)*E

    w = (E1|e2)
    w = w.expand()
    assert w.scalar() == 0

    w = (E1|e3)
    w = w.expand()
    assert w.scalar() == 0

    w = (E2|e1)
    w = w.expand()
    assert w.scalar() == 0

    w = (E2|e3)
    w = w.expand()
    assert w.scalar() == 0

    w = (E3|e1)
    w = w.expand()
    assert w.scalar() == 0

    w = (E3|e2)
    w = w.expand()
    assert w.scalar() == 0

    w = (E1|e1)
    w = (w.expand()).scalar()
    Esq = expand(Esq)
    assert simplify(w/Esq) == 1

    w = (E2|e2)
    w = (w.expand()).scalar()
    assert simplify(w/Esq) == 1

    w = (E3|e3)
    w = (w.expand()).scalar()
    assert simplify(w/Esq) == 1

def test_reciprocal_frame_test():
    with GA_Printer():
        metric = "1 # #," + "# 1 #," + "# # 1,"

        (e1, e2, e3) = MV.setup("e1 e2 e3", metric)

        E = e1 ^ e2 ^ e3
        Esq = (E * E).scalar()
        assert str(E) == "e1^e2^e3"
        assert str(Esq) == "(e1.e2)**2 - 2*(e1.e2)*(e1.e3)*(e2.e3) + (e1.e3)**2 + (e2.e3)**2 - 1"
        Esq_inv = 1 / Esq

        E1 = (e2 ^ e3) * E
        E2 = (-1) * (e1 ^ e3) * E
        E3 = (e1 ^ e2) * E

        assert str(E1) == "((e2.e3)**2 - 1)*e1 + ((e1.e2) - (e1.e3)*(e2.e3))*e2 + (-(e1.e2)*(e2.e3) + (e1.e3))*e3"
        assert str(E2) == "((e1.e2) - (e1.e3)*(e2.e3))*e1 + ((e1.e3)**2 - 1)*e2 + (-(e1.e2)*(e1.e3) + (e2.e3))*e3"
        assert str(E3) == "(-(e1.e2)*(e2.e3) + (e1.e3))*e1 + (-(e1.e2)*(e1.e3) + (e2.e3))*e2 + ((e1.e2)**2 - 1)*e3"

        w = E1 | e2
        w = w.expand()
        assert str(w) == "0"

        w = E1 | e3
        w = w.expand()
        assert str(w) == "0"

        w = E2 | e1
        w = w.expand()
        assert str(w) == "0"

        w = E2 | e3
        w = w.expand()
        assert str(w) == "0"

        w = E3 | e1
        w = w.expand()
        assert str(w) == "0"

        w = E3 | e2
        w = w.expand()
        assert str(w) == "0"

        w = E1 | e1
        w = (w.expand()).scalar()
        Esq = expand(Esq)
        assert str(simplify(w / Esq)) == "1"

        w = E2 | e2
        w = (w.expand()).scalar()
        assert str(simplify(w / Esq)) == "1"

        w = E3 | e3
        w = (w.expand()).scalar()
        assert str(simplify(w / Esq)) == "1"

    return

def _separate(eq, dep, others):
    """Separate expression into two parts based on dependencies of variables."""

    # FIRST PASS
    # Extract derivatives depending our separable variable...
    terms = set()
    for term in eq.args:
        if term.is_Mul:
            for i in term.args:
                if i.is_Derivative and not i.has_any_symbols(*others):
                    terms.add(term)
                    continue
        elif term.is_Derivative and not term.has_any_symbols(*others):
            terms.add(term)
    # Find the factor that we need to divide by
    div = set()
    for term in terms:
        ext, sep = term.expand().as_independent(dep)
        # Failed?
        if sep.has_any_symbols(*others):
            return None
        div.add(ext)
    # FIXME: Find lcm() of all the divisors and divide with it, instead of
    # current hack :(
    # http://code.google.com/p/sympy/issues/detail?id=1498
    if len(div) > 0:
        final = 0
        for term in eq.args:
            eqn = 0
            for i in div:
                eqn += term / i
            final += simplify(eqn)
        eq = final

    # SECOND PASS - separate the derivatives
    div = set()
    lhs = rhs = 0
    for term in eq.args:
        # Check, whether we have already term with independent variable...
        if not term.has_any_symbols(*others):
            lhs += term
            continue
        # ...otherwise, try to separate
        temp, sep = term.expand().as_independent(dep)
        # Failed?
        if sep.has_any_symbols(*others):
            return None
        # Extract the divisors
        div.add(sep)
        rhs -= term.expand()
    # Do the division
    fulldiv = reduce(operator.add, div)
    lhs = simplify(lhs/fulldiv).expand()
    rhs = simplify(rhs/fulldiv).expand()
    # ...and check whether we were successful :)
    if lhs.has_any_symbols(*others) or rhs.has_any_symbols(dep):
        return None
    return [lhs, rhs]

 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 incenter(self):
     """The incenter of the triangle."""
     s = self.sides
     v = self.vertices
     A,B,C = v[0],v[1],v[2]
     a,b,c = s[1].length,s[2].length,s[0].length
     x = simplify( (a*A[0] + b*B[0] + c*C[0]) / (a+b+c) )
     y = simplify( (a*A[1] + b*B[1] + c*C[1]) / (a+b+c) )
     return Point(x, y)

def test_simplify():
    x, y = R2_r.coord_functions()
    dx, dy = R2_r.base_oneforms()
    ex, ey = R2_r.base_vectors()
    assert simplify(x) == x
    assert simplify(x*y) == x*y
    assert simplify(dx*dy) == dx*dy
    assert simplify(ex*ey) == ex*ey
    assert ((1-x)*dx)/(1-x)**2 == dx/(1-x)

	def cleanup(self):
		self.cleanState           = symbols(self.dofNames)
		self.cleanVelocity        = symbols(['d'+name for name in self.dofNames])
		symdict                   = dict(zip(self.state+self.velocity, self.cleanState+self.cleanVelocity))
		cleanMatrix               = Matrix([simplify(element).subs(symdict) for element in  self.mass])
		self.cleanMass            = list([cleanMatrix.reshape(self.degreesOfFreedom, self.degreesOfFreedom)])
		self.cleanForce           = [simplify(element).subs(symdict) for element in self.force]
		self.cleanConstraint      = simplify(self.constraint_equation.subs(symdict))
		self.cleanConstraintForce = [simplify(element).subs(symdict) for element in self.virtualForce]
		self.cleanNullForce       = self.nullForce

 def _image(self, func):
     if isinstance(func, EuclideanTransformation):
         direction = simplify(qrotate(func.rquat, func.parity*self.direction))
         offset = simplify(self.offset + dot(func.tvec, direction))
         closed = self.closed
         return Halfspace(
             offset=offset,
             direction=direction,
             closed=closed,
             normalization=False)

 def __add__(self, other):
     """Add two points, or add a factor to this point's coordinates."""
     if isinstance(other, Point):
         if len(other) == len(self):
             return Point([simplify(a + b) for a, b in zip(self, other)])
         else:
             raise TypeError("Points must have the same number of dimensions")
     else:
         other = sympify(other)
         return Point([simplify(a + other) for a in self])

    def intersection(self, o):
        if isinstance(o, Circle):
            dx,dy = o.center - self.center
            d = sqrt( simplify(dy**2 + dx**2) )
            a = simplify((self.radius**2 - o.radius**2 + d**2) / (2*d))

            x2 = self.center[0] + (dx * a/d)
            y2 = self.center[1] + (dy * a/d)

            h = sqrt( simplify(self.radius**2 - a**2) )
            rx = -dy * (h/d)
            ry =  dx * (h/d)

            xi_1 = simplify(x2 + rx)
            xi_2 = simplify(x2 - rx)
            yi_1 = simplify(y2 + ry)
            yi_2 = simplify(y2 - ry)

            ret = [Point(xi_1, yi_1)]
            if xi_1 != xi_2 or yi_1 != yi_2:
                ret.append(Point(xi_2, yi_2))
            return ret
        elif isinstance(o, Ellipse):
            a, b, r = o.hradius, o.vradius, self.radius
            x = a*sqrt(simplify((r**2 - b**2)/(a**2 - b**2)))
            y = b*sqrt(simplify((a**2 - r**2)/(a**2 - b**2)))
            return list(set([Point(x,y), Point(x,-y), Point(-x,y), Point(-x,-y)]))

        return Ellipse.intersection(self, o)

    def arbitrary_point(self, parameter='t'):
        """A parameterized point on the Segment.

        Parameters
        ==========

        parameter : str, optional
            The name of the parameter which will be used for the parametric
            point. The default value is 't'.

        Returns
        =======

        point : Point


        Parameters
        ==========

        parameter : str, optional
            The name of the parameter which will be used for the parametric
            point. The default value is 't'.

        Returns
        =======

        point : Point

        Raises
        ======

        ValueError
            When ``parameter`` already appears in the Segment's definition.

        See Also
        ========

        sympy.geometry.point.Point

        Examples
        ========

        >>> from sympy import Point, Segment
        >>> p1, p2 = Point(1, 0), Point(5, 3)
        >>> s1 = Segment(p1, p2)
        >>> s1.arbitrary_point()
        Point(4*t + 1, 3*t)

        """
        t = _symbol(parameter)
        if t.name in (f.name for f in self.free_symbols):
            raise ValueError('Symbol %s already appears in object and cannot be used as a parameter.' % t.name)
        x = simplify(self.p1.x + t*(self.p2.x - self.p1.x))
        y = simplify(self.p1.y + t*(self.p2.y - self.p1.y))
        return Point(x, y)