Repository: sympy/sympy
Pull request: Make classes in diffgeom module immutable
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #19321
Fixes #19318

#### Brief description of what is fixed or changed

`Manifold.patches` attribute is deprecated.
`Patch.coord_systems` attribute is deprecated.
Class signature *names* of `CoordSystem` is deprecated. `CoordinateSymbol` class and `CoordSystem.symbols` attribute are introduced instead.
`CoordSystem.transforms` attribute,  `CoordSystem.connect_to` method, and `CoordSystem.coord_tuple_transform_to` method are deprecated. Class signature *relations*, `CoordSystem.transformation` and `CoordSystem.transform` methods are introduced instead.
`CoordSystem.jacobian_determinant` method is introduced.
 methods are introduced instead.

#### Other comments

`Manifold.patches`, `Patch.coord_systems`, and `CoordSystem.transforms` attributes are deprecated because these make classes mutable.  
`CoordSystem._dummies` is replaced with `CoordinateSymbol` class to use `Dummy` less. Also, this makes the coordinate symbols know which coordinate system they belong to.  

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->

- diffgeom
  - `Manifold.patches` attribute is deprecated.
  - `Patch.coord_systems` attribute is deprecated.
  - Class signature *names* of `CoordSystem` is deprecated. `CoordinateSymbol` class and `CoordSystem.symbols` attribute are introduced instead.
  - `CoordSystem.transforms`,  `CoordSystem.connect_to`, and `CoordSystem.coord_tuple_transform_to` are deprecated.  Class signature *relations*, `CoordSystem.transformation` and `CoordSystem.transform` are introduced instead.
  - `CoordSystem.jacobian_determinant` method is introduced.

<!-- END RELEASE NOTES -->

Problem statement:
Idea to implement CoordSystem's _fill_gaps_in_transformations
Currently, `CoordSystem` in `diffgeom` module has unimplemented method called `_fill_gaps_in_transformations`. This method is supposed to establish the transformation relation between two non-directly connected coordinate systems.

I am planning to implement this method. Using path finding algorithm, the shortest path to connect two systems (the coordinate systems between these) will be found. Then, the base scalars will be transformed successively. Once the relation is found, both coordinate systems will be added to each other's `transforms`.

Should there be anything more needed, or something that I should keep in mind, please leave a comment here.

Classes in diffgeom module could be made more strictly immutable
In `diffgeom` module, `Manifold`, `Patch` and `CoordSystem` have mutable features:

1. `Manifold.patches` and `Patch.coord_systems`. They store the patches and coordinate systems defined on themselves, respectively. However, they are not used really, and I doubt that these are necessary.

2. `CoordSystem.transforms`. This actually matters; it stores the coordinate trasform relations. If new coordinate system is defined and connected to existing one, old one's `transforms` will be mutated.

I think that instead of establishing the connection after construction, this `transforms` can be included to `CoordSystem`'s class signal. If new coordinate system B is constructed and connected to the existing one A, both relations "from A to B" and "from B to A" will be stored in `B.transforms`, not A's.

Hints from issue discussion:
I think that some classes in the `diffgeom` are using mutable objects like `transforms` or `coord_systems`, but they shouldn't be using them extensively than cacheing.
And I think that this makes the future unclear about whether `transforms` dict should be replaced by other stuff and be deprecated.
@sylee957 
I totally agree. Especially, I have tried many times to make `transforms` immutable, but none of them were successful. Perhaps I'll make another issue page, and resolve that matter first.
As @sylee957 mentioned in #19318, doing this will deprecate the old behavior. I think we need to, but I'd like to hear other opinions first.

Commit message: Fix References underline length

Patch to review:
diff --git a/doc/src/modules/diffgeom.rst b/doc/src/modules/diffgeom.rst
index 02fc0416af80..764d9f31412d 100644
--- a/doc/src/modules/diffgeom.rst
+++ b/doc/src/modules/diffgeom.rst
@@ -18,6 +18,9 @@ Base Class Reference
 .. autoclass:: CoordSystem
    :members:
 
+.. autoclass:: CoordinateSymbol
+   :members:
+
 .. autoclass:: Point
    :members:
 
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 6bdc452c62c2..0b8109e5854e 100644
--- a/sympy/core/tests/test_args.py
+++ b/sympy/core/tests/test_args.py
@@ -12,7 +12,7 @@
 
 from sympy.testing.pytest import XFAIL, SKIP
 
-x, y, z = symbols('x,y,z')
+a, b, c, x, y, z = symbols('a,b,c,x,y,z')
 
 
 def test_all_classes_are_tested():
@@ -4420,38 +4420,42 @@ def test_sympy__diffgeom__diffgeom__Patch():
 
 def test_sympy__diffgeom__diffgeom__CoordSystem():
     from sympy.diffgeom import Manifold, Patch, CoordSystem
-    assert _test_args(CoordSystem('name', Patch('name', Manifold('name', 3))))
+    assert _test_args(CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c]))
+
+
+def test_sympy__diffgeom__diffgeom__CoordinateSymbol():
+    from sympy.diffgeom import Manifold, Patch, CoordSystem, CoordinateSymbol
+    assert _test_args(CoordinateSymbol(CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c]), 0))
 
 
-@XFAIL
 def test_sympy__diffgeom__diffgeom__Point():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, Point
     assert _test_args(Point(
-        CoordSystem('name', Patch('name', Manifold('name', 3))), [x, y]))
+        CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c]), [x, y]))
 
 
 def test_sympy__diffgeom__diffgeom__BaseScalarField():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     assert _test_args(BaseScalarField(cs, 0))
 
 
 def test_sympy__diffgeom__diffgeom__BaseVectorField():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseVectorField
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     assert _test_args(BaseVectorField(cs, 0))
 
 
 def test_sympy__diffgeom__diffgeom__Differential():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField, Differential
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     assert _test_args(Differential(BaseScalarField(cs, 0)))
 
 
 def test_sympy__diffgeom__diffgeom__Commutator():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseVectorField, Commutator
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
-    cs1 = CoordSystem('name1', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
+    cs1 = CoordSystem('name1', Patch('name', Manifold('name', 3)), [a, b, c])
     v = BaseVectorField(cs, 0)
     v1 = BaseVectorField(cs1, 0)
     assert _test_args(Commutator(v, v1))
@@ -4459,14 +4463,14 @@ def test_sympy__diffgeom__diffgeom__Commutator():
 
 def test_sympy__diffgeom__diffgeom__TensorProduct():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField, Differential, TensorProduct
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     d = Differential(BaseScalarField(cs, 0))
     assert _test_args(TensorProduct(d, d))
 
 
 def test_sympy__diffgeom__diffgeom__WedgeProduct():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField, Differential, WedgeProduct
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     d = Differential(BaseScalarField(cs, 0))
     d1 = Differential(BaseScalarField(cs, 1))
     assert _test_args(WedgeProduct(d, d1))
@@ -4474,7 +4478,7 @@ def test_sympy__diffgeom__diffgeom__WedgeProduct():
 
 def test_sympy__diffgeom__diffgeom__LieDerivative():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField, Differential, BaseVectorField, LieDerivative
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     d = Differential(BaseScalarField(cs, 0))
     v = BaseVectorField(cs, 0)
     assert _test_args(LieDerivative(v, d))
@@ -4483,13 +4487,13 @@ def test_sympy__diffgeom__diffgeom__LieDerivative():
 @XFAIL
 def test_sympy__diffgeom__diffgeom__BaseCovarDerivativeOp():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseCovarDerivativeOp
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     assert _test_args(BaseCovarDerivativeOp(cs, 0, [[[0, ]*3, ]*3, ]*3))
 
 
 def test_sympy__diffgeom__diffgeom__CovarDerivativeOp():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseVectorField, CovarDerivativeOp
-    cs = CoordSystem('name', Patch('name', Manifold('name', 3)))
+    cs = CoordSystem('name', Patch('name', Manifold('name', 3)), [a, b, c])
     v = BaseVectorField(cs, 0)
     _test_args(CovarDerivativeOp(v, [[[0, ]*3, ]*3, ]*3))
 
diff --git a/sympy/diffgeom/__init__.py b/sympy/diffgeom/__init__.py
index 140f00539472..8846a9951060 100644
--- a/sympy/diffgeom/__init__.py
+++ b/sympy/diffgeom/__init__.py
@@ -1,18 +1,18 @@
 from .diffgeom import (
     BaseCovarDerivativeOp, BaseScalarField, BaseVectorField, Commutator,
-    contravariant_order, CoordSystem, CovarDerivativeOp, covariant_order,
-    Differential, intcurve_diffequ, intcurve_series, LieDerivative,
-    Manifold, metric_to_Christoffel_1st, metric_to_Christoffel_2nd,
-    metric_to_Ricci_components, metric_to_Riemann_components, Patch,
-    Point, TensorProduct, twoform_to_matrix, vectors_in_basis,
-    WedgeProduct,
+    contravariant_order, CoordSystem, CoordinateSymbol,
+    CovarDerivativeOp, covariant_order, Differential, intcurve_diffequ,
+    intcurve_series, LieDerivative, Manifold, metric_to_Christoffel_1st,
+    metric_to_Christoffel_2nd, metric_to_Ricci_components,
+    metric_to_Riemann_components, Patch, Point, TensorProduct, twoform_to_matrix,
+    vectors_in_basis, WedgeProduct,
 )
 
 __all__ = [
-    'BaseCovarDerivativeOp', 'BaseScalarField', 'BaseVectorField',
-    'Commutator', 'contravariant_order', 'CoordSystem', 'CovarDerivativeOp',
-    'covariant_order', 'Differential', 'intcurve_diffequ', 'intcurve_series',
-    'LieDerivative', 'Manifold', 'metric_to_Christoffel_1st',
+    'BaseCovarDerivativeOp', 'BaseScalarField', 'BaseVectorField', 'Commutator',
+    'contravariant_order', 'CoordSystem', 'CoordinateSymbol',
+    'CovarDerivativeOp', 'covariant_order', 'Differential', 'intcurve_diffequ',
+    'intcurve_series', 'LieDerivative', 'Manifold', 'metric_to_Christoffel_1st',
     'metric_to_Christoffel_2nd', 'metric_to_Ricci_components',
     'metric_to_Riemann_components', 'Patch', 'Point', 'TensorProduct',
     'twoform_to_matrix', 'vectors_in_basis', 'WedgeProduct',
diff --git a/sympy/diffgeom/diffgeom.py b/sympy/diffgeom/diffgeom.py
index 643b032399c4..d4e87f8fb3df 100644
--- a/sympy/diffgeom/diffgeom.py
+++ b/sympy/diffgeom/diffgeom.py
@@ -4,18 +4,20 @@
 
 from sympy.combinatorics import Permutation
 from sympy.core import (
-    Basic, Expr, Dummy, Function,  diff,
-    Pow, Mul, Add, Atom
+    Basic, Expr, Function, diff,
+    Pow, Mul, Add, Atom, Lambda, S, Tuple, Dict
 )
+from sympy.core.cache import cacheit
 from sympy.core.compatibility import reduce
-from sympy.core.numbers import Zero
+from sympy.core.symbol import Symbol, Dummy
 from sympy.core.symbol import Str
 from sympy.core.sympify import _sympify
 from sympy.functions import factorial
-from sympy.matrices import Matrix
+from sympy.matrices import ImmutableDenseMatrix as Matrix
 from sympy.simplify import simplify
 from sympy.solvers import solve
 
+from sympy.utilities.exceptions import SymPyDeprecationWarning
 
 # TODO you are a bit excessive in the use of Dummies
 # TODO dummy point, literal field
@@ -23,34 +25,58 @@
 # tests and find out why
 from sympy.tensor.array import ImmutableDenseNDimArray
 
-
 class Manifold(Atom):
     """A mathematical manifold.
 
     Explanation
     ===========
 
-    The only role that this object plays is to keep a list of all patches
-    defined on the manifold. It does not provide any means to study the
-    topological characteristics of the manifold that it represents.
+    A manifold is a topological space that locally resembles
+    Euclidean space near each point [1].
+
+    This class does not provide any means to study the topological
+    characteristics of the manifold that it represents, though.
 
     Parameters
     ==========
+
     name : str
         The name of the manifold.
 
     dim : int
         The dimension of the manifold.
+
+    Examples
+    ========
+
+    >>> from sympy.diffgeom import Manifold
+    >>> m = Manifold('M', 2)
+
+    >>> m.name
+    'M'
+    >>> m.dim
+    2
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Manifold
+
     """
 
-    def __new__(cls, name, dim):
+    def __new__(cls, name, dim, **kwargs):
         if not isinstance(name, Str):
             name = Str(name)
         dim = _sympify(dim)
         obj = super().__new__(cls, name, dim)
-        obj.patches = []
-        # The patches list is necessary if a Patch instance needs to enumerate
-        # other Patch instance on the same manifold.
+
+        obj.patches = _deprecated_list(
+            "Manifold.patches",
+            "external container for registry",
+            19321,
+            "1.7",
+            []
+        )
         return obj
 
     @property
@@ -67,17 +93,19 @@ class Patch(Atom):
     Explanation
     ===========
 
-    On a manifold one can have many patches that do not always include the
-    whole manifold. On these patches coordinate charts can be defined that
-    permit the parameterization of any point on the patch in terms of a tuple
-    of real numbers (the coordinates).
+    Coordinate patch, or patch in short, is a simply-connected open set around a point
+    in the manifold [1]. On a manifold one can have many patches that do not always
+    include the whole manifold. On these patches coordinate charts can be defined that
+    permit the parameterization of any point on the patch in terms of a tuple of
+    real numbers (the coordinates).
 
-    This object serves as a container/parent for all coordinate system charts
-    that can be defined on the patch it represents.
+    This class does not provide any means to study the topological
+    characteristics of the patch that it represents.
 
     Parameters
     ==========
-    name : string
+
+    name : str
         The name of the patch.
 
     manifold : Manifold
@@ -86,25 +114,34 @@ class Patch(Atom):
     Examples
     ========
 
-    Define a Manifold and a Patch on that Manifold:
-
     >>> from sympy.diffgeom import Manifold, Patch
-    >>> m = Manifold('M', 3)
+    >>> m = Manifold('M', 2)
     >>> p = Patch('P', m)
-    >>> p in m.patches
-    True
+
+    >>> p.name
+    'P'
+    >>> p.dim
+    2
+
+    References
+    ==========
+
+    .. [1] G. Sussman, J. Wisdom, W. Farr, Functional Differential Geometry (2013)
 
     """
-    # Contains a reference to the parent manifold in order to be able to access
-    # other patches.
-    def __new__(cls, name, manifold):
+    def __new__(cls, name, manifold, **kwargs):
         if not isinstance(name, Str):
             name = Str(name)
         obj = super().__new__(cls, name, manifold)
-        obj.manifold.patches.append(obj)
-        obj.coord_systems = []
-        # The list of coordinate systems is necessary for an instance of
-        # CoordSystem to enumerate other coord systems on the patch.
+
+        obj.manifold.patches.append(obj) # deprecated
+        obj.coord_systems = _deprecated_list(
+            "Patch.coord_systems",
+            "external container for registry",
+            19321,
+            "1.7",
+            []
+        )
         return obj
 
     @property
@@ -120,134 +157,155 @@ def dim(self):
         return self.manifold.dim
 
 class CoordSystem(Atom):
-    """A coordinate system defined on the patch
+    """A coordinate system defined on the patch.
 
     Explanation
     ===========
 
-    This class contains all coordinate transformation logic.
+    Coordinate system is a system that uses one or more coordinates to uniquely determine
+    the position of the points or other geometric elements on a manifold [1].
+
+    By passing Symbols to *symbols* parameter, user can define the name and assumptions
+    of coordinate symbols of the coordinate system.
+    By passing *relations* parameter, user can define the tranform relations of coordinate
+    systems. Inverse transformation and indirect transformation can be found automatically.
+    If this parameter is not passed, coordinate transformation cannot be done.
 
     Parameters
     ==========
 
-    name : string
+    name : str
         The name of the coordinate system.
 
     patch : Patch
         The patch where the coordinate system is defined.
 
-    names : list of strings, optional
-        Determines how base scalar fields will be printed.
+    symbols : list of Symbols
+        Defines the names and assumptions of coordinate symbols.
+
+    relations : dict, optional
+        - key : tuple of two strings, who are the names of systems where
+            the coordinates transform from and transform to.
+        - value : Lambda returning the transformed coordinates.
 
     Examples
     ========
 
-    Define a Manifold and a Patch, and then define two coord systems on that
-    patch:
-
-    >>> from sympy import symbols, sin, cos, pi
+    >>> from sympy import symbols, pi, Lambda, Matrix, sqrt, atan2, cos, sin
     >>> from sympy.diffgeom import Manifold, Patch, CoordSystem
-    >>> from sympy.simplify import simplify
-    >>> r, theta = symbols('r, theta')
     >>> m = Manifold('M', 2)
-    >>> patch = Patch('P', m)
-    >>> rect = CoordSystem('rect', patch)
-    >>> polar = CoordSystem('polar', patch)
-    >>> rect in patch.coord_systems
-    True
+    >>> p = Patch('P', m)
 
-    Connect the coordinate systems. An inverse transformation is automatically
-    found by ``solve`` when possible:
+    >>> x, y = symbols('x y', real=True)
+    >>> r, theta = symbols('r theta', nonnegative=True)
+    >>> relation_dict = {
+    ... ('Car2D', 'Pol'): Lambda((x, y), Matrix([sqrt(x**2 + y**2), atan2(y, x)])),
+    ... ('Pol', 'Car2D'): Lambda((r, theta), Matrix([r*cos(theta), r*sin(theta)]))
+    ... }
+    >>> Car2D = CoordSystem('Car2D', p, [x, y], relation_dict)
+    >>> Pol = CoordSystem('Pol', p, [r, theta], relation_dict)
+
+    >>> Car2D.name
+    'Car2D'
+    >>> Car2D.dim
+    2
+    >>> Car2D.symbols
+    [x, y]
+    >>> Car2D.transformation(Pol)
+    Lambda((x, y), Matrix([
+    [sqrt(x**2 + y**2)],
+    [      atan2(y, x)]]))
 
-    >>> polar.connect_to(rect, [r, theta], [r*cos(theta), r*sin(theta)])
-    >>> polar.coord_tuple_transform_to(rect, [0, 2])
+    >>> Car2D.transform(Pol)
     Matrix([
-    [0],
-    [0]])
-    >>> polar.coord_tuple_transform_to(rect, [2, pi/2])
-    Matrix([
-    [0],
-    [2]])
-    >>> rect.coord_tuple_transform_to(polar, [1, 1]).applyfunc(simplify)
+    [sqrt(x**2 + y**2)],
+    [      atan2(y, x)]])
+    >>> Car2D.transform(Pol, [1, 2])
     Matrix([
-    [sqrt(2)],
-    [   pi/4]])
+    [sqrt(5)],
+    [atan(2)]])
 
-    Calculate the jacobian of the polar to cartesian transformation:
-
-    >>> polar.jacobian(rect, [r, theta])
+    >>> Pol.jacobian_matrix(Car2D)
     Matrix([
     [cos(theta), -r*sin(theta)],
     [sin(theta),  r*cos(theta)]])
-
-    Define a point using coordinates in one of the coordinate systems:
-
-    >>> p = polar.point([1, 3*pi/4])
-    >>> rect.point_to_coords(p)
+    >>> Pol.jacobian_matrix(Car2D, [1, pi/2])
     Matrix([
-    [-sqrt(2)/2],
-    [ sqrt(2)/2]])
-
-    Define a basis scalar field (i.e. a coordinate function), that takes a
-    point and returns its coordinates. It is an instance of ``BaseScalarField``.
-
-    >>> rect.coord_function(0)(p)
-    -sqrt(2)/2
-    >>> rect.coord_function(1)(p)
-    sqrt(2)/2
+    [0, -1],
+    [1,  0]])
 
-    Define a basis vector field (i.e. a unit vector field along the coordinate
-    line). Vectors are also differential operators on scalar fields. It is an
-    instance of ``BaseVectorField``.
-
-    >>> v_x = rect.base_vector(0)
-    >>> x = rect.coord_function(0)
-    >>> v_x(x)
-    1
-    >>> v_x(v_x(x))
-    0
-
-    Define a basis oneform field:
-
-    >>> dx = rect.base_oneform(0)
-    >>> dx(v_x)
-    1
-
-    If you provide a list of names the fields will print nicely:
-    - without provided names:
-
-    >>> x, v_x, dx
-    (rect_0, e_rect_0, drect_0)
-
-    - with provided names
+    References
+    ==========
 
-    >>> rect = CoordSystem('rect', patch, ['x', 'y'])
-    >>> rect.coord_function(0), rect.base_vector(0), rect.base_oneform(0)
-    (x, e_x, dx)
+    .. [1] https://en.wikipedia.org/wiki/Coordinate_system
 
     """
-    #  Contains a reference to the parent patch in order to be able to access
-    # other coordinate system charts.
-    def __new__(cls, name, patch, names=None):
+    def __new__(cls, name, patch, symbols=None, relations={}, **kwargs):
         if not isinstance(name, Str):
             name = Str(name)
-        # names is not in args because it is related only to printing, not to
-        # identifying the CoordSystem instance.
-        obj = super().__new__(cls, name, patch)
-        if not names:
-            names = ['%s_%d' % (name, i) for i in range(patch.dim)]
-        obj._names = tuple(str(i) for i in names)
-        obj.patch.coord_systems.append(obj)
-        obj.transforms = {}
-        # All the coordinate transformation logic is in this dictionary in the
-        # form of:
-        #  key = other coordinate system
-        #  value = tuple of  # TODO make these Lambda instances
-        #          - list of `Dummy` coordinates in this coordinate system
-        #          - list of expressions as a function of the Dummies giving
-        #          the coordinates in another coordinate system
-        obj._dummies = [Dummy(str(n)) for n in names]
+
+        # canonicallize the symbols
+        if symbols is None or any(not isinstance(s, Symbol) for s in symbols) :
+            # support deprecated signature
+            SymPyDeprecationWarning(
+                feature="Class signature 'names' of CoordSystem",
+                useinstead="class signature 'symbols'",
+                issue=19321,
+                deprecated_since_version="1.7"
+            ).warn()
+
+            if symbols is None:
+                names = kwargs.get('names', None)
+            else:
+                names = symbols
+
+            if names is None:
+                symbols = Tuple(*[Symbol('%s_%d' % (name, i), real=True) for i in range(patch.dim)])
+            else:
+                symbols = Tuple(*[Symbol(n, real=True) for n in names])
+        else:
+            symbols = Tuple(*[Symbol(s.name, **s._assumptions.generator) for s in symbols])
+
+        # canonicallize the relations
+        rel_temp = {}
+        for k,v in relations.items():
+            s1, s2 = k
+            if not isinstance(s1, Str):
+                s1 = Str(s1)
+            if not isinstance(s2, Str):
+                s2 = Str(s2)
+            key = Tuple(s1, s2)
+            rel_temp[key] = v
+        relations = Dict(rel_temp)
+
+        # construct the object
+        obj = super().__new__(cls, name, patch, symbols, relations)
+
+        # Add deprecated attributes
+        obj.transforms = _deprecated_dict(
+            "Mutable CoordSystem.transforms",
+            "'relations' parameter in class signature",
+            19321,
+            "1.7",
+            {}
+        )
+        obj._names =_deprecated_list(
+                "CoordSystem._names",
+                "CoordSystem.symbols",
+                19321,
+                "1.7",
+                [str(n) for n in symbols]
+            )
+        obj.patch.coord_systems.append(obj) # deprecated
+        obj._dummies = _deprecated_list(
+            "CoordSystem._dummies",
+            "CoordSystem.symbols",
+            19321,
+            "1.7",
+            [Dummy(str(n)) for n in symbols]
+        )
         obj._dummy = Dummy()
+
         return obj
 
     @property
@@ -258,34 +316,132 @@ def name(self):
     def patch(self):
         return self.args[1]
 
+    @property
+    def manifold(self):
+        return self.patch.manifold
+
+    @property
+    def symbols(self):
+        return [
+            CoordinateSymbol(
+                self, i, **s._assumptions.generator
+            ) for i,s in enumerate(self.args[2])
+        ]
+
+    @property
+    def relations(self):
+        return self.args[3]
+
     @property
     def dim(self):
         return self.patch.dim
 
     ##########################################################################
-    # Coordinate transformations.
+    # Finding transformation relation
     ##########################################################################
 
+    def transformation(self, sys):
+        """
+        Return coordinate transform relation from *self* to *sys* as Lambda.
+        """
+        if self.relations != sys.relations:
+            raise TypeError(
+        "Two coordinate systems have different relations")
+
+        key = Tuple(self.name, sys.name)
+        if key in self.relations:
+            return self.relations[key]
+        elif key[::-1] in self.relations:
+            return self._inverse_transformation(sys, self)
+        else:
+            return self._indirect_transformation(self, sys)
+
+    @staticmethod
+    def _inverse_transformation(sys1, sys2):
+        # Find the transformation relation from sys2 to sys1
+        forward_transform = sys1.transform(sys2)
+        forward_syms, forward_results = forward_transform.args
+
+        inv_syms = [i.as_dummy() for i in forward_syms]
+        inv_results = solve(
+            [t[0] - t[1] for t in zip(inv_syms, forward_results)],
+            list(forward_syms), dict=True)[0]
+        inv_results = [inv_results[s] for s in forward_syms]
+
+        signature = tuple(inv_syms)
+        expr = Matrix(inv_results)
+        return Lambda(signature, expr)
+
+    @classmethod
+    @cacheit
+    def _indirect_transformation(cls, sys1, sys2):
+        # Find the transformation relation between two indirectly connected coordinate systems
+        path = cls._dijkstra(sys1, sys2)
+        Lambdas = []
+        for i in range(len(path) - 1):
+            s1, s2 = path[i], path[i + 1]
+            Lambdas.append(s1.transformation(s2))
+        syms = Lambdas[-1].signature
+        expr = syms
+        for l in reversed(Lambdas):
+            expr = l(*expr)
+        return Lambda(syms, expr)
+
+    @staticmethod
+    def _dijkstra(sys1, sys2):
+        # Use Dijkstra algorithm to find the shortest path between two indirectly-connected
+        # coordinate systems
+        relations = sys1.relations
+        graph = {}
+        for s1, s2 in relations.keys():
+            if s1 not in graph:
+                graph[s1] = {s2}
+            else:
+                graph[s1].add(s2)
+            if s2 not in graph:
+                graph[s2] = {s1}
+            else:
+                graph[s2].add(s1)
+
+        path_dict = {sys:[0, [], 0] for sys in graph} # minimum distance, path, times of visited
+
+        def visit(sys):
+            path_dict[sys][2] = 1
+            for newsys in graph[sys]:
+                distance = path_dict[sys][0] + 1
+                if path_dict[newsys][0] >= distance or not path_dict[newsys][1]:
+                    path_dict[newsys][0] = distance
+                    path_dict[newsys][1] = [i for i in path_dict[sys][1]]
+                    path_dict[newsys][1].append(sys)
+
+        visit(sys1)
+
+        while True:
+            min_distance = max(path_dict.values(), key=lambda x:x[0])[0]
+            newsys = None
+            for sys, lst in path_dict.items():
+                if 0 < lst[0] <= min_distance and not lst[2]:
+                    min_distance = lst[0]
+                    newsys = sys
+            if newsys is None:
+                break
+            visit(newsys)
+
+        result = path_dict[sys2][1]
+        result.append(sys2)
+
+        if result == [sys2]:
+            raise KeyError("Two coordinate systems are not connected.")
+        return result
+
     def connect_to(self, to_sys, from_coords, to_exprs, inverse=True, fill_in_gaps=False):
-        """Register the transformation used to switch to another coordinate system.
-
-        Parameters
-        ==========
-
-        to_sys
-            another instance of ``CoordSystem``
-        from_coords
-            list of symbols in terms of which ``to_exprs`` is given
-        to_exprs
-            list of the expressions of the new coordinate tuple
-        inverse
-            try to deduce and register the inverse transformation
-        fill_in_gaps
-            try to deduce other transformation that are made
-            possible by composing the present transformation with other already
-            registered transformation
+        SymPyDeprecationWarning(
+            feature="CoordSystem.connect_to",
+            useinstead="new instance generated with new 'transforms' parameter",
+            issue=19321,
+            deprecated_since_version="1.7"
+        ).warn()
 
-        """
         from_coords, to_exprs = dummyfy(from_coords, to_exprs)
         self.transforms[to_sys] = Matrix(from_coords), Matrix(to_exprs)
 
@@ -297,6 +453,7 @@ def connect_to(self, to_sys, from_coords, to_exprs, inverse=True, fill_in_gaps=F
 
     @staticmethod
     def _inv_transf(from_coords, to_exprs):
+        # Will be removed when connect_to is removed
         inv_from = [i.as_dummy() for i in from_coords]
         inv_to = solve(
             [t[0] - t[1] for t in zip(inv_from, to_exprs)],
@@ -306,13 +463,49 @@ def _inv_transf(from_coords, to_exprs):
 
     @staticmethod
     def _fill_gaps_in_transformations():
+        # Will be removed when connect_to is removed
         raise NotImplementedError
-        # TODO
+
+    ##########################################################################
+    # Coordinate transformations
+    ##########################################################################
+
+    def transform(self, sys, coordinates=None):
+        """
+        Return the result of coordinate transformation from *self* to *sys*.
+        If coordinates are not given, coordinate symbols of *self* are used.
+        """
+        if coordinates is None:
+            coordinates = Matrix(self.symbols)
+        else:
+            coordinates = Matrix(coordinates)
+        if self != sys:
+            transf = self.transformation(sys)
+            coordinates = transf(*coordinates)
+        return coordinates
+
+    def jacobian_matrix(self, sys, coordinates=None):
+        """
+        Return the jacobian matrix of a transformation.
+        """
+        result = self.transform(sys).jacobian(self.symbols)
+        if coordinates is not None:
+            result = result.subs(list(zip(self.symbols, coordinates)))
+        return result
+
+    def jacobian_determinant(self, sys, coordinates=None):
+        """Return the jacobian determinant of a transformation."""
+        return self.jacobian_matrix(sys, coordinates).det()
 
     def coord_tuple_transform_to(self, to_sys, coords):
-        """Transform ``coords`` to coord system ``to_sys``.
+        """Transform ``coords`` to coord system ``to_sys``."""
+        SymPyDeprecationWarning(
+            feature="CoordSystem.coord_tuple_transform_to",
+            useinstead="CoordSystem.transform",
+            issue=19321,
+            deprecated_since_version="1.7"
+        ).warn()
 
-        See the docstring of ``CoordSystem`` for examples."""
         coords = Matrix(coords)
         if self != to_sys:
             transf = self.transforms[to_sys]
@@ -321,77 +514,122 @@ def coord_tuple_transform_to(self, to_sys, coords):
 
     def jacobian(self, to_sys, coords):
         """Return the jacobian matrix of a transformation."""
+        SymPyDeprecationWarning(
+            feature="CoordSystem.jacobian",
+            useinstead="CoordSystem.jacobian_matrix",
+            issue=19321,
+            deprecated_since_version="1.7"
+        ).warn()
+
         with_dummies = self.coord_tuple_transform_to(
             to_sys, self._dummies).jacobian(self._dummies)
         return with_dummies.subs(list(zip(self._dummies, coords)))
 
     ##########################################################################
-    # Base fields.
+    # Points
     ##########################################################################
 
-    def coord_function(self, coord_index):
-        """Return a ``BaseScalarField`` that takes a point and returns one of the coords.
+    def point(self, coords):
+        """Create a ``Point`` with coordinates given in this coord system."""
+        return Point(self, coords)
+
+    def point_to_coords(self, point):
+        """Calculate the coordinates of a point in this coord system."""
+        return point.coords(self)
 
-        Takes a point and returns its coordinate in this coordinate system.
+    ##########################################################################
+    # Base fields.
+    ##########################################################################
 
-        See the docstring of ``CoordSystem`` for examples."""
+    def base_scalar(self, coord_index):
+        """Return ``BaseScalarField`` that takes a point and returns one of the coordinates."""
         return BaseScalarField(self, coord_index)
+    coord_function = base_scalar
 
-    def coord_functions(self):
+    def base_scalars(self):
         """Returns a list of all coordinate functions.
-
-        For more details see the ``coord_function`` method of this class."""
-        return [self.coord_function(i) for i in range(self.dim)]
+        For more details see the ``base_scalar`` method of this class."""
+        return [self.base_scalar(i) for i in range(self.dim)]
+    coord_functions = base_scalars
 
     def base_vector(self, coord_index):
         """Return a basis vector field.
-
         The basis vector field for this coordinate system. It is also an
-        operator on scalar fields.
-
-        See the docstring of ``CoordSystem`` for examples."""
+        operator on scalar fields."""
         return BaseVectorField(self, coord_index)
 
     def base_vectors(self):
         """Returns a list of all base vectors.
-
         For more details see the ``base_vector`` method of this class."""
         return [self.base_vector(i) for i in range(self.dim)]
 
     def base_oneform(self, coord_index):
         """Return a basis 1-form field.
-
         The basis one-form field for this coordinate system. It is also an
-        operator on vector fields.
-
-        See the docstring of ``CoordSystem`` for examples."""
+        operator on vector fields."""
         return Differential(self.coord_function(coord_index))
 
     def base_oneforms(self):
         """Returns a list of all base oneforms.
-
         For more details see the ``base_oneform`` method of this class."""
         return [self.base_oneform(i) for i in range(self.dim)]
 
-    ##########################################################################
-    # Points.
-    ##########################################################################
+class CoordinateSymbol(Symbol):
+    """A symbol which denotes an abstract value of i-th coordinate of
+    the coordinate system with given context.
 
-    def point(self, coords):
-        """Create a ``Point`` with coordinates given in this coord system.
+    Explanation
+    ===========
 
-        See the docstring of ``CoordSystem`` for examples."""
-        return Point(self, coords)
+    Each coordinates in coordinate system are represented by unique symbol,
+    such as x, y, z in Cartesian coordinate system.
 
-    def point_to_coords(self, point):
-        """Calculate the coordinates of a point in this coord system.
+    You may not construct this class directly. Instead, use `symbols` method
+    of CoordSystem.
 
-        See the docstring of ``CoordSystem`` for examples."""
-        return point.coords(self)
+    Parameters
+    ==========
 
-    ##########################################################################
-    # Printing.
-    ##########################################################################
+    coord_sys : CoordSystem
+
+    index : integer
+
+    Examples
+    ========
+
+    >>> from sympy import symbols
+    >>> from sympy.diffgeom import Manifold, Patch, CoordSystem
+    >>> m = Manifold('M', 2)
+    >>> p = Patch('P', m)
+    >>> _x, _y = symbols('x y', nonnegative=True)
+
+    >>> C = CoordSystem('C', p, [_x, _y])
+    >>> x, y = C.symbols
+
+    >>> x.name
+    'x'
+    >>> x.coord_sys == C
+    True
+    >>> x.index
+    0
+    >>> x.is_nonnegative
+    True
+
+    """
+    def __new__(cls, coord_sys, index, **assumptions):
+        name = coord_sys.args[2][index].name
+        obj = super().__new__(cls, name, **assumptions)
+        obj.coord_sys = coord_sys
+        obj.index = index
+        return obj
+
+    def __getnewargs__(self):
+        return (self.coord_sys, self.index)
+
+    def _hashable_content(self):
+        return (
+            self.coord_sys, self.index
+        ) + tuple(sorted(self.assumptions0.items()))
 
 class Point(Basic):
     """Point defined in a coordinate system.
@@ -399,7 +637,10 @@ class Point(Basic):
     Explanation
     ===========
 
-    To define a point you must supply coordinates and a coordinate system.
+    Mathematically, point is defined in the manifold and does not have any coordinates
+    by itself. Coordinate system is what imbues the coordinates to the point by coordinate
+    chart. However, due to the difficulty of realizing such logic, you must supply
+    a coordinate system and coordinates to define a Point here.
 
     The usage of this object after its definition is independent of the
     coordinate system that was used in order to define it, however due to
@@ -411,67 +652,74 @@ class Point(Basic):
 
     coord_sys : CoordSystem
 
-    coords: list of sympy expressions
+    coords : list
         The coordinates of the point.
 
     Examples
     ========
 
-    Define the boilerplate Manifold, Patch and coordinate systems:
+    >>> from sympy import pi
+    >>> from sympy.diffgeom import Point
+    >>> from sympy.diffgeom.rn import R2, R2_r, R2_p
+    >>> rho, theta = R2_p.symbols
 
-    >>> from sympy import symbols, sin, cos, pi
-    >>> from sympy.diffgeom import (
-    ...        Manifold, Patch, CoordSystem, Point)
-    >>> r, theta = symbols('r, theta')
-    >>> m = Manifold('M', 2)
-    >>> p = Patch('P', m)
-    >>> rect = CoordSystem('rect', p)
-    >>> polar = CoordSystem('polar', p)
-    >>> polar.connect_to(rect, [r, theta], [r*cos(theta), r*sin(theta)])
+    >>> p = Point(R2_p, [rho, 3*pi/4])
 
-    Define a point using coordinates from one of the coordinate systems:
+    >>> p.manifold == R2
+    True
 
-    >>> p = Point(polar, [r, 3*pi/4])
     >>> p.coords()
     Matrix([
-    [     r],
+    [   rho],
     [3*pi/4]])
-    >>> p.coords(rect)
+    >>> p.coords(R2_r)
     Matrix([
-    [-sqrt(2)*r/2],
-    [ sqrt(2)*r/2]])
+    [-sqrt(2)*rho/2],
+    [ sqrt(2)*rho/2]])
 
     """
-    def __new__(cls, coord_sys, coords):
+
+    def __new__(cls, coord_sys, coords, **kwargs):
         coords = Matrix(coords)
         obj = super().__new__(cls, coord_sys, coords)
         obj._coord_sys = coord_sys
         obj._coords = coords
         return obj
 
-    def coords(self, to_sys=None):
-        """Coordinates of the point in a given coordinate system.
+    @property
+    def patch(self):
+        return self._coord_sys.patch
 
-        If ``to_sys`` is ``None`` it returns the coordinates in the system in
-        which the point was defined."""
-        if to_sys:
-            return self._coord_sys.coord_tuple_transform_to(to_sys, self._coords)
-        else:
+    @property
+    def manifold(self):
+        return self._coord_sys.manifold
+
+    @property
+    def dim(self):
+        return self.manifold.dim
+
+    def coords(self, sys=None):
+        """
+        Coordinates of the point in given coordinate system. If coordinate system
+        is not passed, it returns the coordinates in the coordinate system in which
+        the poin was defined.
+        """
+        if sys is None:
             return self._coords
+        else:
+            return self._coord_sys.transform(sys, self._coords)
 
     @property
     def free_symbols(self):
         return self._coords.free_symbols
 
-
 class BaseScalarField(Expr):
-    """Base Scalar Field over a Manifold for a given Coordinate System.
+    """Base scalar field over a manifold for a given coordinate system.
 
     Explanation
     ===========
 
     A scalar field takes a point as an argument and returns a scalar.
-
     A base scalar field of a coordinate system takes a point and returns one of
     the coordinates of that point in the coordinate system in question.
 
@@ -482,7 +730,6 @@ class BaseScalarField(Expr):
     coordinate system in which it was defined, however due to limitations in
     the simplification routines you may arrive at more complicated
     expression if you use unappropriate coordinate systems.
-
     You can build complicated scalar fields by just building up SymPy
     expressions containing ``BaseScalarField`` instances.
 
@@ -496,30 +743,23 @@ class BaseScalarField(Expr):
     Examples
     ========
 
-    Define boilerplate Manifold, Patch and coordinate systems:
-
-    >>> from sympy import symbols, sin, cos, pi, Function
-    >>> from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
-    >>> r0, theta0 = symbols('r0, theta0')
-    >>> m = Manifold('M', 2)
-    >>> p = Patch('P', m)
-    >>> rect = CoordSystem('rect', p)
-    >>> polar = CoordSystem('polar', p)
-    >>> polar.connect_to(rect, [r0, theta0], [r0*cos(theta0), r0*sin(theta0)])
-
-    Point to be used as an argument for the filed:
-
-    >>> point = polar.point([r0, 0])
-
-    Examples of fields:
+    >>> from sympy import Function, pi
+    >>> from sympy.diffgeom import BaseScalarField
+    >>> from sympy.diffgeom.rn import R2_r, R2_p
+    >>> rho, _ = R2_p.symbols
+    >>> point = R2_p.point([rho, 0])
+    >>> fx, fy = R2_r.base_scalars()
+    >>> ftheta = BaseScalarField(R2_r, 1)
+
+    >>> fx(point)
+    rho
+    >>> fy(point)
+    0
 
-    >>> fx = BaseScalarField(rect, 0)
-    >>> fy = BaseScalarField(rect, 1)
     >>> (fx**2+fy**2).rcall(point)
-    r0**2
+    rho**2
 
     >>> g = Function('g')
-    >>> ftheta = BaseScalarField(polar, 1)
     >>> fg = g(ftheta-pi)
     >>> fg.rcall(point)
     g(-pi)
@@ -528,16 +768,35 @@ class BaseScalarField(Expr):
 
     is_commutative = True
 
-    def __new__(cls, coord_sys, index):
+    def __new__(cls, coord_sys, index, **kwargs):
         index = _sympify(index)
         obj = super().__new__(cls, coord_sys, index)
         obj._coord_sys = coord_sys
         obj._index = index
         return obj
 
+    @property
+    def coord_sys(self):
+        return self.args[0]
+
+    @property
+    def index(self):
+        return self.args[1]
+
+    @property
+    def patch(self):
+        return self.coord_sys.patch
+
+    @property
+    def manifold(self):
+        return self.coord_sys.manifold
+
+    @property
+    def dim(self):
+        return self.manifold.dim
+
     def __call__(self, *args):
         """Evaluating the field at a point or doing nothing.
-
         If the argument is a ``Point`` instance, the field is evaluated at that
         point. The field is returned itself if the argument is any other
         object. It is so in order to have working recursive calling mechanics
@@ -559,14 +818,13 @@ def doit(self):
 
 
 class BaseVectorField(Expr):
-    r"""Vector Field over a Manifold.
+    r"""Base vector field over a manifold for a given coordinate system.
 
     Explanation
     ===========
 
     A vector field is an operator taking a scalar field and returning a
     directional derivative (which is also a scalar field).
-
     A base vector field is the same type of operator, however the derivation is
     specifically done with respect to a chosen coordinate.
 
@@ -580,7 +838,6 @@ class BaseVectorField(Expr):
 
     Parameters
     ==========
-
     coord_sys : CoordSystem
 
     index : integer
@@ -588,29 +845,19 @@ class BaseVectorField(Expr):
     Examples
     ========
 
-    Use the predefined R2 manifold, setup some boilerplate.
-
-    >>> from sympy import symbols, Function
-    >>> from sympy.diffgeom.rn import R2, R2_p, R2_r
+    >>> from sympy import Function
+    >>> from sympy.diffgeom.rn import R2_p, R2_r
     >>> from sympy.diffgeom import BaseVectorField
     >>> from sympy import pprint
-    >>> x0, y0, r0, theta0 = symbols('x0, y0, r0, theta0')
 
-    Points to be used as arguments for the field:
-
-    >>> point_p = R2_p.point([r0, theta0])
-    >>> point_r = R2_r.point([x0, y0])
-
-    Scalar field to operate on:
+    >>> x, y = R2_r.symbols
+    >>> rho, theta = R2_p.symbols
+    >>> fx, fy = R2_r.base_scalars()
+    >>> point_p = R2_p.point([rho, theta])
+    >>> point_r = R2_r.point([x, y])
 
     >>> g = Function('g')
-    >>> s_field = g(R2.x, R2.y)
-    >>> s_field.rcall(point_r)
-    g(x0, y0)
-    >>> s_field.rcall(point_p)
-    g(r0*cos(theta0), r0*sin(theta0))
-
-    Vector field:
+    >>> s_field = g(fx, fy)
 
     >>> v = BaseVectorField(R2_r, 1)
     >>> pprint(v(s_field))
@@ -618,31 +865,49 @@ class BaseVectorField(Expr):
     |---(g(x, xi))||
     \dxi          /|xi=y
     >>> pprint(v(s_field).rcall(point_r).doit())
-     d
-    ---(g(x0, y0))
-    dy0
+    d
+    --(g(x, y))
+    dy
     >>> pprint(v(s_field).rcall(point_p))
     / d                        \|
-    |---(g(r0*cos(theta0), xi))||
-    \dxi                       /|xi=r0*sin(theta0)
+    |---(g(rho*cos(theta), xi))||
+    \dxi                       /|xi=rho*sin(theta)
 
     """
 
     is_commutative = False
 
-    def __new__(cls, coord_sys, index):
+    def __new__(cls, coord_sys, index, **kwargs):
         index = _sympify(index)
         obj = super().__new__(cls, coord_sys, index)
         obj._coord_sys = coord_sys
         obj._index = index
         return obj
 
+    @property
+    def coord_sys(self):
+        return self.args[0]
+
+    @property
+    def index(self):
+        return self.args[1]
+
+    @property
+    def patch(self):
+        return self.coord_sys.patch
+
+    @property
+    def manifold(self):
+        return self.coord_sys.manifold
+
+    @property
+    def dim(self):
+        return self.manifold.dim
+
     def __call__(self, scalar_field):
         """Apply on a scalar field.
-
         The action of a vector field on a scalar field is a directional
         differentiation.
-
         If the argument is not a scalar field an error is raised.
         """
         if covariant_order(scalar_field) or contravariant_order(scalar_field):
@@ -662,11 +927,11 @@ def __call__(self, scalar_field):
         d_result = d_result.diff(d_var)
 
         # Second step: e_x(x) -> 1 and e_x(r) -> cos(atan2(x, y))
-        coords = self._coord_sys._dummies
+        coords = self._coord_sys.symbols
         d_funcs_deriv = [f.diff(d_var) for f in d_funcs]
         d_funcs_deriv_sub = []
         for b in base_scalars:
-            jac = self._coord_sys.jacobian(b._coord_sys, coords)
+            jac = self._coord_sys.jacobian_matrix(b._coord_sys, coords)
             d_funcs_deriv_sub.append(jac[b._index, self._index])
         d_result = d_result.subs(list(zip(d_funcs_deriv, d_funcs_deriv_sub)))
 
@@ -686,6 +951,9 @@ def _find_coords(expr):
 class Commutator(Expr):
     r"""Commutator of two vector fields.
 
+    Explanation
+    ===========
+
     The commutator of two vector fields `v_1` and `v_2` is defined as the
     vector field `[v_1, v_2]` that evaluated on each scalar field `f` is equal
     to `v_1(v_2(f)) - v_2(v_1(f))`.
@@ -693,15 +961,15 @@ class Commutator(Expr):
     Examples
     ========
 
-    Use the predefined R2 manifold, setup some boilerplate.
 
-    >>> from sympy.diffgeom.rn import R2
+    >>> from sympy.diffgeom.rn import R2_p, R2_r
     >>> from sympy.diffgeom import Commutator
     >>> from sympy.simplify import simplify
 
-    Vector fields:
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> e_r = R2_p.base_vector(0)
 
-    >>> e_x, e_y, e_r = R2.e_x, R2.e_y, R2.e_r
     >>> c_xy = Commutator(e_x, e_y)
     >>> c_xr = Commutator(e_x, e_r)
     >>> c_xy
@@ -710,9 +978,8 @@ class Commutator(Expr):
     Unfortunately, the current code is not able to compute everything:
 
     >>> c_xr
-    Commutator(e_x, e_r)
-
-    >>> simplify(c_xr(R2.y**2))
+    Commutator(e_x, e_rho)
+    >>> simplify(c_xr(fy**2))
     -2*cos(theta)*y**2/(x**2 + y**2)
 
     """
@@ -722,13 +989,13 @@ def __new__(cls, v1, v2):
             raise ValueError(
                 'Only commutators of vector fields are supported.')
         if v1 == v2:
-            return Zero()
+            return S.Zero
         coord_sys = set().union(*[_find_coords(v) for v in (v1, v2)])
         if len(coord_sys) == 1:
             # Only one coordinate systems is used, hence it is easy enough to
             # actually evaluate the commutator.
             if all(isinstance(v, BaseVectorField) for v in (v1, v2)):
-                return Zero()
+                return S.Zero
             bases_1, bases_2 = [list(v.atoms(BaseVectorField))
                                 for v in (v1, v2)]
             coeffs_1 = [v1.expand().coeff(b) for b in bases_1]
@@ -739,53 +1006,51 @@ def __new__(cls, v1, v2):
                     res += c1*b1(c2)*b2 - c2*b2(c1)*b1
             return res
         else:
-            return super().__new__(cls, v1, v2)
+            obj = super().__new__(cls, v1, v2)
+            obj._v1 = v1 # deprecated assignment
+            obj._v2 = v2 # deprecated assignment
+            return obj
 
-    def __init__(self, v1, v2):
-        super().__init__()
-        self._args = (v1, v2)
-        self._v1 = v1
-        self._v2 = v2
+    @property
+    def v1(self):
+        return self.args[0]
+
+    @property
+    def v2(self):
+        return self.args[1]
 
     def __call__(self, scalar_field):
         """Apply on a scalar field.
-
         If the argument is not a scalar field an error is raised.
         """
-        return self._v1(self._v2(scalar_field)) - self._v2(self._v1(scalar_field))
+        return self.v1(self.v2(scalar_field)) - self.v2(self.v1(scalar_field))
 
 
 class Differential(Expr):
     r"""Return the differential (exterior derivative) of a form field.
 
+    Explanation
+    ===========
+
     The differential of a form (i.e. the exterior derivative) has a complicated
     definition in the general case.
-
     The differential `df` of the 0-form `f` is defined for any vector field `v`
     as `df(v) = v(f)`.
 
     Examples
     ========
 
-    Use the predefined R2 manifold, setup some boilerplate.
-
     >>> from sympy import Function
-    >>> from sympy.diffgeom.rn import R2
+    >>> from sympy.diffgeom.rn import R2_r
     >>> from sympy.diffgeom import Differential
     >>> from sympy import pprint
 
-    Scalar field (0-forms):
-
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
     >>> g = Function('g')
-    >>> s_field = g(R2.x, R2.y)
-
-    Vector fields:
-
-    >>> e_x, e_y, = R2.e_x, R2.e_y
-
-    Differentials:
-
+    >>> s_field = g(fx, fy)
     >>> dg = Differential(s_field)
+
     >>> dg
     d(g(x, y))
     >>> pprint(dg(e_x))
@@ -801,7 +1066,6 @@ class Differential(Expr):
 
     >>> Differential(dg)
     0
-
     """
 
     is_commutative = False
@@ -811,18 +1075,22 @@ def __new__(cls, form_field):
             raise ValueError(
                 'A vector field was supplied as an argument to Differential.')
         if isinstance(form_field, Differential):
-            return Zero()
+            return S.Zero
         else:
-            return super().__new__(cls, form_field)
+            obj = super().__new__(cls, form_field)
+            obj._form_field = form_field # deprecated assignment
+            return obj
 
-    def __init__(self, form_field):
-        super().__init__()
-        self._form_field = form_field
-        self._args = (self._form_field, )
+    @property
+    def form_field(self):
+        return self.args[0]
 
     def __call__(self, *vector_fields):
         """Apply on a list of vector_fields.
 
+        Explanation
+        ===========
+
         If the number of vector fields supplied is not equal to 1 + the order of
         the form field inside the differential the result is undefined.
 
@@ -867,10 +1135,12 @@ def __call__(self, *vector_fields):
                         ret += (-1)**(i + j)*t
             return ret
 
-
 class TensorProduct(Expr):
     """Tensor product of forms.
 
+    Explanation
+    ===========
+
     The tensor product permits the creation of multilinear functionals (i.e.
     higher order tensors) out of lower order fields (e.g. 1-forms and vector
     fields). However, the higher tensors thus created lack the interesting
@@ -880,27 +1150,28 @@ class TensorProduct(Expr):
     Examples
     ========
 
-    Use the predefined R2 manifold, setup some boilerplate.
-
-    >>> from sympy.diffgeom.rn import R2
+    >>> from sympy.diffgeom.rn import R2_r
     >>> from sympy.diffgeom import TensorProduct
 
-    >>> TensorProduct(R2.dx, R2.dy)(R2.e_x, R2.e_y)
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> dx, dy = R2_r.base_oneforms()
+
+    >>> TensorProduct(dx, dy)(e_x, e_y)
     1
-    >>> TensorProduct(R2.dx, R2.dy)(R2.e_y, R2.e_x)
+    >>> TensorProduct(dx, dy)(e_y, e_x)
     0
-    >>> TensorProduct(R2.dx, R2.x*R2.dy)(R2.x*R2.e_x, R2.e_y)
+    >>> TensorProduct(dx, fx*dy)(fx*e_x, e_y)
     x**2
-    >>> TensorProduct(R2.e_x, R2.e_y)(R2.x**2, R2.y**2)
+    >>> TensorProduct(e_x, e_y)(fx**2, fy**2)
     4*x*y
-    >>> TensorProduct(R2.e_y, R2.dx)(R2.y)
+    >>> TensorProduct(e_y, dx)(fy)
     dx
 
-
     You can nest tensor products.
 
-    >>> tp1 = TensorProduct(R2.dx, R2.dy)
-    >>> TensorProduct(tp1, R2.dx)(R2.e_x, R2.e_y, R2.e_x)
+    >>> tp1 = TensorProduct(dx, dy)
+    >>> TensorProduct(tp1, dx)(e_x, e_y, e_x)
     1
 
     You can make partial contraction for instance when 'raising an index'.
@@ -908,13 +1179,13 @@ class TensorProduct(Expr):
     respective position in the tensor product is left as it is.
 
     >>> TP = TensorProduct
-    >>> metric = TP(R2.dx, R2.dx) + 3*TP(R2.dy, R2.dy)
-    >>> metric.rcall(R2.e_y, None)
+    >>> metric = TP(dx, dx) + 3*TP(dy, dy)
+    >>> metric.rcall(e_y, None)
     3*dy
 
     Or automatically pad the args with ``None`` without specifying them.
 
-    >>> metric.rcall(R2.e_y)
+    >>> metric.rcall(e_y)
     3*dy
 
     """
@@ -928,10 +1199,6 @@ def __new__(cls, *args):
         else:
             return scalar
 
-    def __init__(self, *args):
-        super().__init__()
-        self._args = args
-
     def __call__(self, *fields):
         """Apply on a list of fields.
 
@@ -942,6 +1209,7 @@ def __call__(self, *fields):
         the forms inside the tensor product. The sublists are provided as
         arguments to these forms and the resulting expressions are given to the
         constructor of ``TensorProduct``.
+
         """
         tot_order = covariant_order(self) + contravariant_order(self)
         tot_args = len(fields)
@@ -957,30 +1225,35 @@ def __call__(self, *fields):
 class WedgeProduct(TensorProduct):
     """Wedge product of forms.
 
+    Explanation
+    ===========
+
     In the context of integration only completely antisymmetric forms make
     sense. The wedge product permits the creation of such forms.
 
     Examples
     ========
 
-    Use the predefined R2 manifold, setup some boilerplate.
-
-    >>> from sympy.diffgeom.rn import R2
+    >>> from sympy.diffgeom.rn import R2_r
     >>> from sympy.diffgeom import WedgeProduct
 
-    >>> WedgeProduct(R2.dx, R2.dy)(R2.e_x, R2.e_y)
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> dx, dy = R2_r.base_oneforms()
+
+    >>> WedgeProduct(dx, dy)(e_x, e_y)
     1
-    >>> WedgeProduct(R2.dx, R2.dy)(R2.e_y, R2.e_x)
+    >>> WedgeProduct(dx, dy)(e_y, e_x)
     -1
-    >>> WedgeProduct(R2.dx, R2.x*R2.dy)(R2.x*R2.e_x, R2.e_y)
+    >>> WedgeProduct(dx, fx*dy)(fx*e_x, e_y)
     x**2
-    >>> WedgeProduct(R2.e_x,R2.e_y)(R2.y,None)
+    >>> WedgeProduct(e_x, e_y)(fy, None)
     -e_x
 
     You can nest wedge products.
 
-    >>> wp1 = WedgeProduct(R2.dx, R2.dy)
-    >>> WedgeProduct(wp1, R2.dx)(R2.e_x, R2.e_y, R2.e_x)
+    >>> wp1 = WedgeProduct(dx, dy)
+    >>> WedgeProduct(wp1, dx)(e_x, e_y, e_x)
     0
 
     """
@@ -988,7 +1261,6 @@ class WedgeProduct(TensorProduct):
     # TODO you do not need all these permutations (neither the prefactor)
     def __call__(self, *fields):
         """Apply on a list of vector_fields.
-
         The expression is rewritten internally in terms of tensor products and evaluated."""
         orders = (covariant_order(e) + contravariant_order(e) for e in self.args)
         mul = 1/Mul(*(factorial(o) for o in orders))
@@ -1002,6 +1274,9 @@ def __call__(self, *fields):
 class LieDerivative(Expr):
     """Lie derivative with respect to a vector field.
 
+    Explanation
+    ===========
+
     The transport operator that defines the Lie derivative is the pushforward of
     the field to be derived along the integral curve of the field with respect
     to which one derives.
@@ -1009,27 +1284,35 @@ class LieDerivative(Expr):
     Examples
     ========
 
+    >>> from sympy.diffgeom.rn import R2_r, R2_p
     >>> from sympy.diffgeom import (LieDerivative, TensorProduct)
-    >>> from sympy.diffgeom.rn import R2
-    >>> LieDerivative(R2.e_x, R2.y)
+
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> e_rho, e_theta = R2_p.base_vectors()
+    >>> dx, dy = R2_r.base_oneforms()
+
+    >>> LieDerivative(e_x, fy)
     0
-    >>> LieDerivative(R2.e_x, R2.x)
+    >>> LieDerivative(e_x, fx)
     1
-    >>> LieDerivative(R2.e_x, R2.e_x)
+    >>> LieDerivative(e_x, e_x)
     0
 
     The Lie derivative of a tensor field by another tensor field is equal to
     their commutator:
 
-    >>> LieDerivative(R2.e_x, R2.e_r)
-    Commutator(e_x, e_r)
-    >>> LieDerivative(R2.e_x + R2.e_y, R2.x)
+    >>> LieDerivative(e_x, e_rho)
+    Commutator(e_x, e_rho)
+    >>> LieDerivative(e_x + e_y, fx)
     1
-    >>> tp = TensorProduct(R2.dx, R2.dy)
-    >>> LieDerivative(R2.e_x, tp)
+
+    >>> tp = TensorProduct(dx, dy)
+    >>> LieDerivative(e_x, tp)
     LieDerivative(e_x, TensorProduct(dx, dy))
-    >>> LieDerivative(R2.e_x, tp)
+    >>> LieDerivative(e_x, tp)
     LieDerivative(e_x, TensorProduct(dx, dy))
+
     """
     def __new__(cls, v_field, expr):
         expr_form_ord = covariant_order(expr)
@@ -1038,21 +1321,27 @@ def __new__(cls, v_field, expr):
                              ' vector fields. The supplied argument was not a '
                              'vector field.')
         if expr_form_ord > 0:
-            return super().__new__(cls, v_field, expr)
+            obj = super().__new__(cls, v_field, expr)
+            # deprecated assignments
+            obj._v_field = v_field
+            obj._expr = expr
+            return obj
         if expr.atoms(BaseVectorField):
             return Commutator(v_field, expr)
         else:
             return v_field.rcall(expr)
 
-    def __init__(self, v_field, expr):
-        super().__init__()
-        self._v_field = v_field
-        self._expr = expr
-        self._args = (self._v_field, self._expr)
+    @property
+    def v_field(self):
+        return self.args[0]
+
+    @property
+    def expr(self):
+        return self.args[1]
 
     def __call__(self, *args):
-        v = self._v_field
-        expr = self._expr
+        v = self.v_field
+        expr = self.expr
         lead_term = v(expr(*args))
         rest = Add(*[Mul(*args[:i] + (Commutator(v, args[i]),) + args[i + 1:])
                      for i in range(len(args))])
@@ -1065,32 +1354,52 @@ class BaseCovarDerivativeOp(Expr):
     Examples
     ========
 
-    >>> from sympy.diffgeom.rn import R2, R2_r
+    >>> from sympy.diffgeom.rn import R2_r
     >>> from sympy.diffgeom import BaseCovarDerivativeOp
     >>> from sympy.diffgeom import metric_to_Christoffel_2nd, TensorProduct
+
     >>> TP = TensorProduct
-    >>> ch = metric_to_Christoffel_2nd(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> dx, dy = R2_r.base_oneforms()
+
+    >>> ch = metric_to_Christoffel_2nd(TP(dx, dx) + TP(dy, dy))
     >>> ch
     [[[0, 0], [0, 0]], [[0, 0], [0, 0]]]
     >>> cvd = BaseCovarDerivativeOp(R2_r, 0, ch)
-    >>> cvd(R2.x)
+    >>> cvd(fx)
     1
-    >>> cvd(R2.x*R2.e_x)
+    >>> cvd(fx*e_x)
     e_x
     """
-    def __init__(self, coord_sys, index, christoffel):
-        super().__init__()
-        self._coord_sys = coord_sys
-        self._index = index
-        self._christoffel = christoffel
-        self._args = self._coord_sys, self._index, self._christoffel
+
+    def __new__(cls, coord_sys, index, christoffel):
+        index = _sympify(index)
+        christoffel = ImmutableDenseNDimArray(christoffel)
+        obj = super().__new__(cls, coord_sys, index, christoffel)
+        # deprecated assignments
+        obj._coord_sys = coord_sys
+        obj._index = index
+        obj._christoffel = christoffel
+        return obj
+
+    @property
+    def coord_sys(self):
+        return self.args[0]
+
+    @property
+    def index(self):
+        return self.args[1]
+
+    @property
+    def christoffel(self):
+        return self.args[2]
 
     def __call__(self, field):
         """Apply on a scalar field.
 
         The action of a vector field on a scalar field is a directional
         differentiation.
-
         If the argument is not a scalar field the behaviour is undefined.
         """
         if covariant_order(field) != 0:
@@ -1136,31 +1445,44 @@ class CovarDerivativeOp(Expr):
     Examples
     ========
 
-    >>> from sympy.diffgeom.rn import R2
+    >>> from sympy.diffgeom.rn import R2_r
     >>> from sympy.diffgeom import CovarDerivativeOp
     >>> from sympy.diffgeom import metric_to_Christoffel_2nd, TensorProduct
     >>> TP = TensorProduct
-    >>> ch = metric_to_Christoffel_2nd(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
+    >>> fx, fy = R2_r.base_scalars()
+    >>> e_x, e_y = R2_r.base_vectors()
+    >>> dx, dy = R2_r.base_oneforms()
+    >>> ch = metric_to_Christoffel_2nd(TP(dx, dx) + TP(dy, dy))
+
     >>> ch
     [[[0, 0], [0, 0]], [[0, 0], [0, 0]]]
-    >>> cvd = CovarDerivativeOp(R2.x*R2.e_x, ch)
-    >>> cvd(R2.x)
+    >>> cvd = CovarDerivativeOp(fx*e_x, ch)
+    >>> cvd(fx)
     x
-    >>> cvd(R2.x*R2.e_x)
+    >>> cvd(fx*e_x)
     x*e_x
 
     """
-    def __init__(self, wrt, christoffel):
-        super().__init__()
+
+    def __new__(cls, wrt, christoffel):
         if len({v._coord_sys for v in wrt.atoms(BaseVectorField)}) > 1:
             raise NotImplementedError()
         if contravariant_order(wrt) != 1 or covariant_order(wrt):
             raise ValueError('Covariant derivatives are defined only with '
                              'respect to vector fields. The supplied argument '
                              'was not a vector field.')
-        self._wrt = wrt
-        self._christoffel = christoffel
-        self._args = self._wrt, self._christoffel
+        obj = super().__new__(cls, wrt, christoffel)
+        # deprecated assigments
+        obj._wrt = wrt
+        obj._christoffel = christoffel
+        return obj
+
+    @property
+    def wrt(self):
+        return self.args[0]
+
+    def christoffel(self):
+        return self.args[1]
 
     def __call__(self, field):
         vectors = list(self._wrt.atoms(BaseVectorField))
@@ -1175,6 +1497,9 @@ def __call__(self, field):
 def intcurve_series(vector_field, param, start_point, n=6, coord_sys=None, coeffs=False):
     r"""Return the series expansion for an integral curve of the field.
 
+    Explanation
+    ===========
+
     Integral curve is a function `\gamma` taking a parameter in `R` to a point
     in the manifold. It verifies the equation:
 
@@ -1184,7 +1509,6 @@ def intcurve_series(vector_field, param, start_point, n=6, coord_sys=None, coeff
     value `t` for the parameter and any scalar field `f`.
 
     This equation can also be decomposed of a basis of coordinate functions
-
     `V(f_i)\big(\gamma(t)\big) = \frac{d}{dt}f_i\big(\gamma(t)\big) \quad \forall i`
 
     This function returns a series expansion of `\gamma(t)` in terms of the
@@ -1193,22 +1517,20 @@ def intcurve_series(vector_field, param, start_point, n=6, coord_sys=None, coeff
     represent movement between points on the manifold (i.e. there is no such
     thing as a difference of points for a general manifold).
 
-    See Also
-    ========
-
-    intcurve_diffequ
-
     Parameters
     ==========
-
     vector_field
         the vector field for which an integral curve will be given
+
     param
         the argument of the function `\gamma` from R to the curve
+
     start_point
         the point which corresponds to `\gamma(0)`
+
     n
         the order to which to expand
+
     coord_sys
         the coordinate system in which to expand
         coeffs (default False) - if True return a list of elements of the expansion
@@ -1267,6 +1589,11 @@ def intcurve_series(vector_field, param, start_point, n=6, coord_sys=None, coeff
     [t**2*(-x**2/(x**2 + y**2)**(3/2) + 1/sqrt(x**2 + y**2))/2],
     [                                t**2*x*y/(x**2 + y**2)**2]])
 
+    See Also
+    ========
+
+    intcurve_diffequ
+
     """
     if contravariant_order(vector_field) != 1 or covariant_order(vector_field):
         raise ValueError('The supplied field was not a vector field.')
@@ -1291,6 +1618,9 @@ def taylor_terms_per_coord(coord_function):
 def intcurve_diffequ(vector_field, param, start_point, coord_sys=None):
     r"""Return the differential equation for an integral curve of the field.
 
+    Explanation
+    ===========
+
     Integral curve is a function `\gamma` taking a parameter in `R` to a point
     in the manifold. It verifies the equation:
 
@@ -1305,25 +1635,24 @@ def intcurve_diffequ(vector_field, param, start_point, coord_sys=None):
     represent movement between points on the manifold (i.e. there is no such
     thing as a difference of points for a general manifold).
 
-    See Also
-    ========
-
-    intcurve_series
-
     Parameters
     ==========
 
     vector_field
         the vector field for which an integral curve will be given
+
     param
         the argument of the function `\gamma` from R to the curve
+
     start_point
         the point which corresponds to `\gamma(0)`
+
     coord_sys
         the coordinate system in which to give the equations
 
     Returns
     =======
+
     a tuple of (equations, initial conditions)
 
     Examples
@@ -1356,6 +1685,11 @@ def intcurve_diffequ(vector_field, param, start_point, coord_sys=None):
     >>> init_cond
     [f_0(0) - 1, f_1(0) - pi/2]
 
+    See Also
+    ========
+
+    intcurve_series
+
     """
     if contravariant_order(vector_field) != 1 or covariant_order(vector_field):
         raise ValueError('The supplied field was not a vector field.')
@@ -1394,6 +1728,7 @@ def contravariant_order(expr, _strict=False):
     >>> from sympy.diffgeom import contravariant_order
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.abc import a
+
     >>> contravariant_order(a)
     0
     >>> contravariant_order(a*R2.x + 2)
@@ -1439,6 +1774,7 @@ def covariant_order(expr, _strict=False):
     >>> from sympy.diffgeom import covariant_order
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.abc import a
+
     >>> covariant_order(a)
     0
     >>> covariant_order(a*R2.x + 2)
@@ -1480,7 +1816,6 @@ def covariant_order(expr, _strict=False):
 ###############################################################################
 def vectors_in_basis(expr, to_sys):
     """Transform all base vectors in base vectors of a specified coord basis.
-
     While the new base vectors are in the new coordinate system basis, any
     coefficients are kept in the old system.
 
@@ -1489,16 +1824,18 @@ def vectors_in_basis(expr, to_sys):
 
     >>> from sympy.diffgeom import vectors_in_basis
     >>> from sympy.diffgeom.rn import R2_r, R2_p
+
     >>> vectors_in_basis(R2_r.e_x, R2_p)
-    -y*e_theta/(x**2 + y**2) + x*e_r/sqrt(x**2 + y**2)
+    -y*e_theta/(x**2 + y**2) + x*e_rho/sqrt(x**2 + y**2)
     >>> vectors_in_basis(R2_p.e_r, R2_r)
     sin(theta)*e_y + cos(theta)*e_x
+
     """
     vectors = list(expr.atoms(BaseVectorField))
     new_vectors = []
     for v in vectors:
         cs = v._coord_sys
-        jac = cs.jacobian(to_sys, cs.coord_functions())
+        jac = cs.jacobian_matrix(to_sys, cs.coord_functions())
         new = (jac.T*Matrix(to_sys.base_vectors()))[v._index]
         new_vectors.append(new)
     return expr.subs(list(zip(vectors, new_vectors)))
@@ -1520,6 +1857,7 @@ def twoform_to_matrix(expr):
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import twoform_to_matrix, TensorProduct
     >>> TP = TensorProduct
+
     >>> twoform_to_matrix(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
     Matrix([
     [1, 0],
@@ -1551,7 +1889,6 @@ def twoform_to_matrix(expr):
 
 def metric_to_Christoffel_1st(expr):
     """Return the nested list of Christoffel symbols for the given metric.
-
     This returns the Christoffel symbol of first kind that represents the
     Levi-Civita connection for the given metric.
 
@@ -1561,6 +1898,7 @@ def metric_to_Christoffel_1st(expr):
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import metric_to_Christoffel_1st, TensorProduct
     >>> TP = TensorProduct
+
     >>> metric_to_Christoffel_1st(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
     [[[0, 0], [0, 0]], [[0, 0], [0, 0]]]
     >>> metric_to_Christoffel_1st(R2.x*TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
@@ -1584,7 +1922,6 @@ def metric_to_Christoffel_1st(expr):
 
 def metric_to_Christoffel_2nd(expr):
     """Return the nested list of Christoffel symbols for the given metric.
-
     This returns the Christoffel symbol of second kind that represents the
     Levi-Civita connection for the given metric.
 
@@ -1594,6 +1931,7 @@ def metric_to_Christoffel_2nd(expr):
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import metric_to_Christoffel_2nd, TensorProduct
     >>> TP = TensorProduct
+
     >>> metric_to_Christoffel_2nd(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
     [[[0, 0], [0, 0]], [[0, 0], [0, 0]]]
     >>> metric_to_Christoffel_2nd(R2.x*TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
@@ -1611,7 +1949,7 @@ def metric_to_Christoffel_2nd(expr):
     for e in matrix:
         s_fields.update(e.atoms(BaseScalarField))
     s_fields = list(s_fields)
-    dums = coord_sys._dummies
+    dums = coord_sys.symbols
     matrix = matrix.subs(list(zip(s_fields, dums))).inv().subs(list(zip(dums, s_fields)))
     # XXX end of workaround
     christoffel = [[[Add(*[matrix[i, l]*ch_1st[l, j, k] for l in indices])
@@ -1635,18 +1973,18 @@ def metric_to_Riemann_components(expr):
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import metric_to_Riemann_components, TensorProduct
     >>> TP = TensorProduct
+
     >>> metric_to_Riemann_components(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
     [[[[0, 0], [0, 0]], [[0, 0], [0, 0]]], [[[0, 0], [0, 0]], [[0, 0], [0, 0]]]]
-
     >>> non_trivial_metric = exp(2*R2.r)*TP(R2.dr, R2.dr) + \
         R2.r**2*TP(R2.dtheta, R2.dtheta)
     >>> non_trivial_metric
-    exp(2*r)*TensorProduct(dr, dr) + r**2*TensorProduct(dtheta, dtheta)
+    exp(2*rho)*TensorProduct(drho, drho) + rho**2*TensorProduct(dtheta, dtheta)
     >>> riemann = metric_to_Riemann_components(non_trivial_metric)
     >>> riemann[0, :, :, :]
-    [[[0, 0], [0, 0]], [[0, exp(-2*r)*r], [-exp(-2*r)*r, 0]]]
+    [[[0, 0], [0, 0]], [[0, exp(-2*rho)*rho], [-exp(-2*rho)*rho, 0]]]
     >>> riemann[1, :, :, :]
-    [[[0, -1/r], [1/r, 0]], [[0, 0], [0, 0]]]
+    [[[0, -1/rho], [1/rho, 0]], [[0, 0], [0, 0]]]
 
     """
     ch_2nd = metric_to_Christoffel_2nd(expr)
@@ -1676,6 +2014,7 @@ def metric_to_Riemann_components(expr):
 
 
 def metric_to_Ricci_components(expr):
+
     """Return the components of the Ricci tensor expressed in a given basis.
 
     Given a metric it calculates the components of the Ricci tensor in the
@@ -1689,15 +2028,15 @@ def metric_to_Ricci_components(expr):
     >>> from sympy.diffgeom.rn import R2
     >>> from sympy.diffgeom import metric_to_Ricci_components, TensorProduct
     >>> TP = TensorProduct
+
     >>> metric_to_Ricci_components(TP(R2.dx, R2.dx) + TP(R2.dy, R2.dy))
     [[0, 0], [0, 0]]
-
     >>> non_trivial_metric = exp(2*R2.r)*TP(R2.dr, R2.dr) + \
                              R2.r**2*TP(R2.dtheta, R2.dtheta)
     >>> non_trivial_metric
-    exp(2*r)*TensorProduct(dr, dr) + r**2*TensorProduct(dtheta, dtheta)
+    exp(2*rho)*TensorProduct(drho, drho) + rho**2*TensorProduct(dtheta, dtheta)
     >>> metric_to_Ricci_components(non_trivial_metric)
-    [[1/r, 0], [0, exp(-2*r)*r]]
+    [[1/rho, 0], [0, exp(-2*rho)*rho]]
 
     """
     riemann = metric_to_Riemann_components(expr)
@@ -1707,3 +2046,43 @@ def metric_to_Ricci_components(expr):
               for j in indices]
              for i in indices]
     return ImmutableDenseNDimArray(ricci)
+
+###############################################################################
+# Classes for deprecation
+###############################################################################
+
+class _deprecated_container(object):
+    # This class gives deprecation warning.
+    # When deprecated features are completely deleted, this should be removed as well.
+    # See https://github.com/sympy/sympy/pull/19368
+    def __init__(self, feature, useinstead, issue, version, data):
+        super().__init__(data)
+        self.feature = feature
+        self.useinstead = useinstead
+        self.issue = issue
+        self.version = version
+
+    def warn(self):
+        SymPyDeprecationWarning(
+                    feature=self.feature,
+                    useinstead=self.useinstead,
+                    issue=self.issue,
+                    deprecated_since_version=self.version).warn()
+
+    def __iter__(self):
+        self.warn()
+        return super().__iter__()
+
+    def __getitem__(self, key):
+        self.warn()
+        return super().__getitem__(key)
+
+    def __contains__(self, key):
+        self.warn()
+        return super().__contains__(key)
+
+class _deprecated_list(_deprecated_container, list):
+    pass
+
+class _deprecated_dict(_deprecated_container, dict):
+    pass
diff --git a/sympy/diffgeom/rn.py b/sympy/diffgeom/rn.py
index 1226a63f0377..55c964ccbd2b 100644
--- a/sympy/diffgeom/rn.py
+++ b/sympy/diffgeom/rn.py
@@ -9,28 +9,44 @@
 """
 
 from typing import Any
+import warnings
 
+from sympy import sqrt, atan2, acos, sin, cos, Lambda, Matrix, symbols, Dummy
 from .diffgeom import Manifold, Patch, CoordSystem
-from sympy import sqrt, atan2, acos, sin, cos, Dummy
+
+__all__ = [
+    'R2', 'R2_origin', 'relations_2d', 'R2_r', 'R2_p',
+    'R3', 'R3_origin', 'relations_3d', 'R3_r', 'R3_c', 'R3_s'
+]
 
 ###############################################################################
 # R2
 ###############################################################################
 R2 = Manifold('R^2', 2)  # type: Any
-# Patch and coordinate systems.
+
 R2_origin = Patch('origin', R2)  # type: Any
-R2_r = CoordSystem('rectangular', R2_origin, ['x', 'y'])  # type: Any
-R2_p = CoordSystem('polar', R2_origin, ['r', 'theta'])  # type: Any
-
-# Connecting the coordinate charts.
-x, y, r, theta = [Dummy(s) for s in ['x', 'y', 'r', 'theta']]
-R2_r.connect_to(R2_p, [x, y],
-                      [sqrt(x**2 + y**2), atan2(y, x)],
-                inverse=False, fill_in_gaps=False)
-R2_p.connect_to(R2_r, [r, theta],
-                      [r*cos(theta), r*sin(theta)],
-                inverse=False, fill_in_gaps=False)
-del x, y, r, theta
+
+x, y = symbols('x y', real=True)
+r, theta = symbols('rho theta', nonnegative=True)
+
+relations_2d = {
+    ('rectangular', 'polar'): Lambda((x, y), Matrix([sqrt(x**2 + y**2), atan2(y, x)])),
+    ('polar', 'rectangular'): Lambda((r, theta), Matrix([r*cos(theta), r*sin(theta)])),
+}
+
+R2_r = CoordSystem('rectangular', R2_origin, [x, y], relations_2d)  # type: Any
+R2_p = CoordSystem('polar', R2_origin, [r, theta], relations_2d)  # type: Any
+
+# support deprecated feature
+with warnings.catch_warnings():
+    warnings.simplefilter("ignore")
+    x, y, r, theta = symbols('x y r theta', cls=Dummy)
+    R2_r.connect_to(R2_p, [x, y],
+                        [sqrt(x**2 + y**2), atan2(y, x)],
+                    inverse=False, fill_in_gaps=False)
+    R2_p.connect_to(R2_r, [r, theta],
+                        [r*cos(theta), r*sin(theta)],
+                    inverse=False, fill_in_gaps=False)
 
 # Defining the basis coordinate functions and adding shortcuts for them to the
 # manifold and the patch.
@@ -51,39 +67,69 @@
 # R3
 ###############################################################################
 R3 = Manifold('R^3', 3)  # type: Any
-# Patch and coordinate systems.
+
 R3_origin = Patch('origin', R3)  # type: Any
-R3_r = CoordSystem('rectangular', R3_origin, ['x', 'y', 'z'])  # type: Any
-R3_c = CoordSystem('cylindrical', R3_origin, ['rho', 'psi', 'z'])  # type: Any
-R3_s = CoordSystem('spherical', R3_origin, ['r', 'theta', 'phi'])  # type: Any
-
-# Connecting the coordinate charts.
-x, y, z, rho, psi, r, theta, phi = [Dummy(s) for s in ['x', 'y', 'z',
-                                          'rho', 'psi', 'r', 'theta', 'phi']]
-## rectangular <-> cylindrical
-R3_r.connect_to(R3_c, [x, y, z],
-                      [sqrt(x**2 + y**2), atan2(y, x), z],
-                inverse=False, fill_in_gaps=False)
-R3_c.connect_to(R3_r, [rho, psi, z],
-                      [rho*cos(psi), rho*sin(psi), z],
-                inverse=False, fill_in_gaps=False)
-## rectangular <-> spherical
-R3_r.connect_to(R3_s, [x, y, z],
-                      [sqrt(x**2 + y**2 + z**2), acos(z/
-                            sqrt(x**2 + y**2 + z**2)), atan2(y, x)],
-                inverse=False, fill_in_gaps=False)
-R3_s.connect_to(R3_r, [r, theta, phi],
-                      [r*sin(theta)*cos(phi), r*sin(
-                          theta)*sin(phi), r*cos(theta)],
-                inverse=False, fill_in_gaps=False)
-## cylindrical <-> spherical
-R3_c.connect_to(R3_s, [rho, psi, z],
-                      [sqrt(rho**2 + z**2), acos(z/sqrt(rho**2 + z**2)), psi],
-                inverse=False, fill_in_gaps=False)
-R3_s.connect_to(R3_c, [r, theta, phi],
-                      [r*sin(theta), phi, r*cos(theta)],
-                inverse=False, fill_in_gaps=False)
-del x, y, z, rho, psi, r, theta, phi
+
+x, y, z = symbols('x y z', real=True)
+rho, psi, r, theta, phi = symbols('rho psi r theta phi', nonnegative=True)
+
+relations_3d = {
+    ('rectangular', 'cylindrical'):
+        Lambda((x, y, z), Matrix([sqrt(x**2 + y**2), atan2(y, x), z])),
+    ('cylindrical', 'rectangular'):
+        Lambda((rho, psi, z), Matrix([rho*cos(psi), rho*sin(psi), z])),
+    ('rectangular', 'spherical'):
+        Lambda(
+            (x, y, z),
+            Matrix([
+                sqrt(x**2 + y**2 + z**2),
+                acos(z/sqrt(x**2 + y**2 + z**2)),
+                atan2(y, x)])
+        ),
+    ('spherical', 'rectangular'):
+        Lambda(
+            (r, theta, phi),
+            Matrix([r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta)])
+        ),
+    ('cylindrical', 'spherical'):
+        Lambda(
+            (rho, psi, z),
+            Matrix([sqrt(rho**2 + z**2), acos(z/sqrt(rho**2 + z**2)), psi])
+        ),
+    ('spherical', 'cylindrical'):
+        Lambda((r, theta, phi), Matrix([r*sin(theta), phi, r*cos(theta)])),
+}
+
+R3_r = CoordSystem('rectangular', R3_origin, [x, y, z], relations_3d)  # type: Any
+R3_c = CoordSystem('cylindrical', R3_origin, [rho, psi, z], relations_3d)  # type: Any
+R3_s = CoordSystem('spherical', R3_origin, [r, theta, phi], relations_3d)  # type: Any
+
+# support deprecated feature
+with warnings.catch_warnings():
+    warnings.simplefilter("ignore")
+    x, y, z, rho, psi, r, theta, phi = symbols('x y z rho psi r theta phi', cls=Dummy)
+    R3_r.connect_to(R3_c, [x, y, z],
+                        [sqrt(x**2 + y**2), atan2(y, x), z],
+                    inverse=False, fill_in_gaps=False)
+    R3_c.connect_to(R3_r, [rho, psi, z],
+                        [rho*cos(psi), rho*sin(psi), z],
+                    inverse=False, fill_in_gaps=False)
+    ## rectangular <-> spherical
+    R3_r.connect_to(R3_s, [x, y, z],
+                        [sqrt(x**2 + y**2 + z**2), acos(z/
+                                sqrt(x**2 + y**2 + z**2)), atan2(y, x)],
+                    inverse=False, fill_in_gaps=False)
+    R3_s.connect_to(R3_r, [r, theta, phi],
+                        [r*sin(theta)*cos(phi), r*sin(
+                            theta)*sin(phi), r*cos(theta)],
+                    inverse=False, fill_in_gaps=False)
+    ## cylindrical <-> spherical
+    R3_c.connect_to(R3_s, [rho, psi, z],
+                        [sqrt(rho**2 + z**2), acos(z/sqrt(rho**2 + z**2)), psi],
+                    inverse=False, fill_in_gaps=False)
+    R3_s.connect_to(R3_c, [r, theta, phi],
+                        [r*sin(theta), phi, r*cos(theta)],
+                    inverse=False, fill_in_gaps=False)
 
 # Defining the basis coordinate functions.
 R3_r.x, R3_r.y, R3_r.z = R3_r.coord_functions()
diff --git a/sympy/diffgeom/tests/test_class_structure.py b/sympy/diffgeom/tests/test_class_structure.py
index 519eed5bfa42..68ab461ed61e 100644
--- a/sympy/diffgeom/tests/test_class_structure.py
+++ b/sympy/diffgeom/tests/test_class_structure.py
@@ -1,33 +1,23 @@
 from sympy.diffgeom import Manifold, Patch, CoordSystem, Point
 from sympy import symbols, Function
+from sympy.testing.pytest import warns_deprecated_sympy
 
 m = Manifold('m', 2)
 p = Patch('p', m)
-cs = CoordSystem('cs', p, ['a', 'b'])
-cs_noname = CoordSystem('cs', p)
+a, b = symbols('a b')
+cs = CoordSystem('cs', p, [a, b])
 x, y = symbols('x y')
 f = Function('f')
 s1, s2 = cs.coord_functions()
 v1, v2 = cs.base_vectors()
 f1, f2 = cs.base_oneforms()
 
-
 def test_point():
     point = Point(cs, [x, y])
-    assert point == point.func(*point.args)
     assert point != Point(cs, [2, y])
     #TODO assert point.subs(x, 2) == Point(cs, [2, y])
     #TODO assert point.free_symbols == set([x, y])
 
-def test_rebuild():
-    assert m == m.func(*m.args)
-    assert p == p.func(*p.args)
-    assert cs == cs.func(*cs.args)
-    assert cs_noname == cs_noname.func(*cs_noname.args)
-    assert s1 == s1.func(*s1.args)
-    assert v1 == v1.func(*v1.args)
-    assert f1 == f1.func(*f1.args)
-
 def test_subs():
     assert s1.subs(s1, s2) == s2
     assert v1.subs(v1, v2) == v2
@@ -35,3 +25,10 @@ def test_subs():
     assert (x*f(s1) + y).subs(s1, s2) == x*f(s2) + y
     assert (f(s1)*v1).subs(v1, v2) == f(s1)*v2
     assert (y*f(s1)*f1).subs(f1, f2) == y*f(s1)*f2
+
+def test_deprecated():
+    with warns_deprecated_sympy():
+        cs_wname = CoordSystem('cs', p, ['a', 'b'])
+        cs_noname = CoordSystem('cs', p)
+        assert cs_wname == cs_wname.func(*cs_wname.args)
+        assert cs_noname == cs_wname.func(*cs_noname.args)
diff --git a/sympy/diffgeom/tests/test_diffgeom.py b/sympy/diffgeom/tests/test_diffgeom.py
index 7c63cc67e5d1..0c9cc397cb12 100644
--- a/sympy/diffgeom/tests/test_diffgeom.py
+++ b/sympy/diffgeom/tests/test_diffgeom.py
@@ -9,6 +9,7 @@
 from sympy.functions import sqrt, atan2, sin
 from sympy.matrices import Matrix
 from sympy.testing.pytest import raises, nocache_fail
+from sympy.testing.pytest import warns_deprecated_sympy
 
 TP = TensorProduct
 
@@ -26,23 +27,38 @@ def test_R2():
     # 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)
 
+    #TODO assert m == R2_r.transform(R2_p, R2_p.transform(R2_r, [a, b])).applyfunc(simplify)
+    assert m == R2_p.transform(R2_r, R2_r.transform(R2_p, m)).applyfunc(simplify)
+
+    # deprecated method
+    with warns_deprecated_sympy():
+        assert m == R2_p.coord_tuple_transform_to(
+            R2_r, R2_r.coord_tuple_transform_to(R2_p, m)).applyfunc(simplify)
 
 def test_R3():
     a, b, c = symbols('a b c', positive=True)
     m = Matrix([[a], [b], [c]])
-    assert m == R3_c.coord_tuple_transform_to(
-        R3_r, R3_r.coord_tuple_transform_to(R3_c, m)).applyfunc(simplify)
-    #TODO assert m == R3_r.coord_tuple_transform_to(R3_c, R3_c.coord_tuple_transform_to(R3_r, m)).applyfunc(simplify)
-    assert m == R3_s.coord_tuple_transform_to(
-        R3_r, R3_r.coord_tuple_transform_to(R3_s, m)).applyfunc(simplify)
-    #TODO assert m == R3_r.coord_tuple_transform_to(R3_s, R3_s.coord_tuple_transform_to(R3_r, m)).applyfunc(simplify)
-    assert m == R3_s.coord_tuple_transform_to(
-        R3_c, R3_c.coord_tuple_transform_to(R3_s, m)).applyfunc(simplify)
-    #TODO assert m == R3_c.coord_tuple_transform_to(R3_s, R3_s.coord_tuple_transform_to(R3_c, m)).applyfunc(simplify)
+
+    assert m == R3_c.transform(R3_r, R3_r.transform(R3_c, m)).applyfunc(simplify)
+    #TODO assert m == R3_r.transform(R3_c, R3_c.transform(R3_r, m)).applyfunc(simplify)
+    assert m == R3_s.transform(
+        R3_r, R3_r.transform(R3_s, m)).applyfunc(simplify)
+    #TODO assert m == R3_r.transform(R3_s, R3_s.transform(R3_r, m)).applyfunc(simplify)
+    assert m == R3_s.transform(
+        R3_c, R3_c.transform(R3_s, m)).applyfunc(simplify)
+    #TODO assert m == R3_c.transform(R3_s, R3_s.transform(R3_c, m)).applyfunc(simplify)
+
+    with warns_deprecated_sympy():
+        assert m == R3_c.coord_tuple_transform_to(
+            R3_r, R3_r.coord_tuple_transform_to(R3_c, m)).applyfunc(simplify)
+        #TODO assert m == R3_r.coord_tuple_transform_to(R3_c, R3_c.coord_tuple_transform_to(R3_r, m)).applyfunc(simplify)
+        assert m == R3_s.coord_tuple_transform_to(
+            R3_r, R3_r.coord_tuple_transform_to(R3_s, m)).applyfunc(simplify)
+        #TODO assert m == R3_r.coord_tuple_transform_to(R3_s, R3_s.coord_tuple_transform_to(R3_r, m)).applyfunc(simplify)
+        assert m == R3_s.coord_tuple_transform_to(
+            R3_c, R3_c.coord_tuple_transform_to(R3_s, m)).applyfunc(simplify)
+        #TODO assert m == R3_c.coord_tuple_transform_to(R3_s, R3_s.coord_tuple_transform_to(R3_c, m)).applyfunc(simplify)
 
 
 def test_point():
diff --git a/sympy/diffgeom/tests/test_function_diffgeom_book.py b/sympy/diffgeom/tests/test_function_diffgeom_book.py
index a7317a7d1191..20297cd6e1eb 100644
--- a/sympy/diffgeom/tests/test_function_diffgeom_book.py
+++ b/sympy/diffgeom/tests/test_function_diffgeom_book.py
@@ -4,6 +4,7 @@
 from sympy.simplify import trigsimp, simplify
 from sympy.functions import sqrt, atan2, sin, cos
 from sympy.matrices import Matrix
+from sympy.testing.pytest import warns_deprecated_sympy
 
 # Most of the functionality is covered in the
 # test_functional_diffgeom_ch* tests which are based on the
@@ -25,8 +26,12 @@ def test_functional_diffgeom_ch2():
     assert (R2_r.point_to_coords(R2_p.point([r0, theta0])) ==
            Matrix([r0*cos(theta0), r0*sin(theta0)]))
 
-    assert R2_p.jacobian(R2_r, [r0, theta0]) == Matrix(
-        [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])
+    assert R2_p.jacobian_matrix(R2_r, [r0, theta0]) == Matrix(
+            [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])
+
+    with warns_deprecated_sympy():
+        assert R2_p.jacobian(R2_r, [r0, theta0]) == Matrix(
+            [[cos(theta0), -r0*sin(theta0)], [sin(theta0), r0*cos(theta0)]])
 
     field = f(R2.x, R2.y)
     p1_in_rect = R2_r.point([x0, y0])
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index c695b344fc5a..6d28f6005909 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -2506,31 +2506,52 @@ def _print_MatrixHomomorphism(self, h):
             self._print(h.domain), self._print(h.codomain))
 
     def _print_Manifold(self, manifold):
-        return r'\text{%s}' % manifold.name
+        string = manifold.name.name
+        if '{' in string:
+            name, supers, subs = string, [], []
+        else:
+            name, supers, subs = split_super_sub(string)
+
+            name = translate(name)
+            supers = [translate(sup) for sup in supers]
+            subs = [translate(sub) for sub in subs]
+
+        name = r'\text{%s}' % name
+        if supers:
+            name += "^{%s}" % " ".join(supers)
+        if subs:
+            name += "_{%s}" % " ".join(subs)
+
+        return name
 
     def _print_Patch(self, patch):
-        return r'\text{%s}_{\text{%s}}' % (patch.name, patch.manifold.name)
+        return r'\text{%s}_{%s}' % (self._print(patch.name), self._print(patch.manifold))
+
+    def _print_CoordSystem(self, coordsys):
+        return r'\text{%s}^{\text{%s}}_{%s}' % (
+            self._print(coordsys.name), self._print(coordsys.patch.name), self._print(coordsys.manifold)
+        )
 
-    def _print_CoordSystem(self, coords):
-        return r'\text{%s}^{\text{%s}}_{\text{%s}}' % (
-            coords.name, coords.patch.name, coords.patch.manifold.name
+    def _print_Point(self, point):
+        return '{%s}_{%s}' % (
+            self._print(point.coordinates()), self._print(point._coord_sys)
         )
 
     def _print_CovarDerivativeOp(self, cvd):
         return r'\mathbb{\nabla}_{%s}' % self._print(cvd._wrt)
 
     def _print_BaseScalarField(self, field):
-        string = field._coord_sys._names[field._index]
+        string = field._coord_sys.symbols[field._index].name
         return r'\mathbf{{{}}}'.format(self._print(Symbol(string)))
 
     def _print_BaseVectorField(self, field):
-        string = field._coord_sys._names[field._index]
+        string = field._coord_sys.symbols[field._index].name
         return r'\partial_{{{}}}'.format(self._print(Symbol(string)))
 
     def _print_Differential(self, diff):
         field = diff._form_field
         if hasattr(field, '_coord_sys'):
-            string = field._coord_sys._names[field._index]
+            string = field._coord_sys.symbols[field._index].name
             return r'\operatorname{{d}}{}'.format(self._print(Symbol(string)))
         else:
             string = self._print(field)
diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index 5fe19d241f16..c4c3630fd38c 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -2608,17 +2608,17 @@ def _print_CoordSystem(self, coords):
         return self._print(coords.name)
 
     def _print_BaseScalarField(self, field):
-        string = field._coord_sys._names[field._index]
+        string = field._coord_sys.symbols[field._index].name
         return self._print(pretty_symbol(string))
 
     def _print_BaseVectorField(self, field):
-        s = U('PARTIAL DIFFERENTIAL') + '_' + field._coord_sys._names[field._index]
+        s = U('PARTIAL DIFFERENTIAL') + '_' + field._coord_sys.symbols[field._index].name
         return self._print(pretty_symbol(s))
 
     def _print_Differential(self, diff):
         field = diff._form_field
         if hasattr(field, '_coord_sys'):
-            string = field._coord_sys._names[field._index]
+            string = field._coord_sys.symbols[field._index].name
             return self._print(u'\N{DOUBLE-STRUCK ITALIC SMALL D} ' + pretty_symbol(string))
         else:
             pform = self._print(field)
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index a7cdbd79e3d2..742fb34b81b1 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -7051,11 +7051,12 @@ def test_issue_18272():
 
 def test_diffgeom():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
+    x,y = symbols('x y', real=True)
     m = Manifold('M', 2)
     assert pretty(m) == 'M'
     p = Patch('P', m)
     assert pretty(p) == "P"
-    rect = CoordSystem('rect', p)
+    rect = CoordSystem('rect', p, [x, y])
     assert pretty(rect) == "rect"
     b = BaseScalarField(rect, 0)
-    assert pretty(b) == "rect_0"
+    assert pretty(b) == "x"
diff --git a/sympy/printing/repr.py b/sympy/printing/repr.py
index d156c6f3121d..1d33861800e9 100644
--- a/sympy/printing/repr.py
+++ b/sympy/printing/repr.py
@@ -248,6 +248,24 @@ def _print_Symbol(self, expr):
             return "%s(%s, %s)" % (expr.__class__.__name__,
                                    self._print(expr.name), ', '.join(attr))
 
+    def _print_CoordinateSymbol(self, expr):
+        d = expr._assumptions.generator
+
+        if d == {}:
+            return "%s(%s, %s)" % (
+                expr.__class__.__name__,
+                self._print(expr.coordinate_system),
+                self._print(expr.index)
+            )
+        else:
+            attr = ['%s=%s' % (k, v) for k, v in d.items()]
+            return "%s(%s, %s, %s)" % (
+                expr.__class__.__name__,
+                self._print(expr.coordinate_system),
+                self._print(expr.index),
+                ', '.join(attr)
+            )
+
     def _print_Predicate(self, expr):
         return "%s(%s)" % (expr.__class__.__name__, self._print(expr.name))
 
@@ -328,19 +346,6 @@ def _print_ExtensionElement(self, f):
         ext = self._print(f.ext)
         return "ExtElem(%s, %s)" % (rep, ext)
 
-    def _print_CoordSystem(self, coords):
-        class_name = coords.func.__name__
-        name = self._print(coords.name)
-        patch = self._print(coords.patch)
-        names = self._print(coords._names)
-        return "%s(%s, %s, %s)" % (class_name, name, patch, names)
-
-    def _print_BaseScalarField(self, bsf):
-        class_name = bsf.func.__name__
-        coords = self._print(bsf._coord_sys)
-        idx = self._print(bsf._index)
-        return "%s(%s, %s)" % (class_name, coords, idx)
-
 def srepr(expr, **settings):
     """return expr in repr form"""
     return ReprPrinter(settings).doprint(expr)
diff --git a/sympy/printing/str.py b/sympy/printing/str.py
index 2639a9a5d570..16e5d8cf905c 100644
--- a/sympy/printing/str.py
+++ b/sympy/printing/str.py
@@ -876,24 +876,24 @@ def _print_Category(self, category):
         return 'Category("%s")' % category.name
 
     def _print_Manifold(self, manifold):
-        return manifold.name
+        return self._print(manifold.name)
 
     def _print_Patch(self, patch):
-        return patch.name
+        return self._print(patch.name)
 
     def _print_CoordSystem(self, coords):
-        return coords.name
+        return self._print(coords.name)
 
     def _print_BaseScalarField(self, field):
-        return field._coord_sys._names[field._index]
+        return field._coord_sys.symbols[field._index].name
 
     def _print_BaseVectorField(self, field):
-        return 'e_%s' % field._coord_sys._names[field._index]
+        return 'e_%s' % field._coord_sys.symbols[field._index].name
 
     def _print_Differential(self, diff):
         field = diff._form_field
         if hasattr(field, '_coord_sys'):
-            return 'd%s' % field._coord_sys._names[field._index]
+            return 'd%s' % field._coord_sys.symbols[field._index].name
         else:
             return 'd(%s)' % self._print(field)
 
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index 9c2d37e97546..675ac0ab76d5 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -2535,14 +2535,15 @@ def test_text_re_im():
 def test_latex_diffgeom():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField, Differential
     from sympy.diffgeom.rn import R2
+    x,y = symbols('x y', real=True)
     m = Manifold('M', 2)
     assert latex(m) == r'\text{M}'
     p = Patch('P', m)
     assert latex(p) == r'\text{P}_{\text{M}}'
-    rect = CoordSystem('rect', p)
+    rect = CoordSystem('rect', p, [x, y])
     assert latex(rect) == r'\text{rect}^{\text{P}}_{\text{M}}'
     b = BaseScalarField(rect, 0)
-    assert latex(b) ==  r'\mathbf{rect_{0}}'
+    assert latex(b) ==  r'\mathbf{x}'
 
     g = Function('g')
     s_field = g(R2.x, R2.y)
diff --git a/sympy/printing/tests/test_repr.py b/sympy/printing/tests/test_repr.py
index ed892c1801ca..481547cd2646 100644
--- a/sympy/printing/tests/test_repr.py
+++ b/sympy/printing/tests/test_repr.py
@@ -312,16 +312,6 @@ def test_Permutation():
     import_stmt = "from sympy.combinatorics import Permutation"
     sT(Permutation(1, 2), "Permutation(1, 2)", import_stmt)
 
-
-def test_diffgeom():
-    from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
-    m = Manifold('M', 2)
-    p = Patch('P', m)
-    rect = CoordSystem('rect', p)
-    assert srepr(rect) == "CoordSystem(Str('rect'), Patch(Str('P'), Manifold(Str('M'), Integer(2))), ('rect_0', 'rect_1'))"
-    b = BaseScalarField(rect, 0)
-    assert srepr(b) == "BaseScalarField(CoordSystem(Str('rect'), Patch(Str('P'), Manifold(Str('M'), Integer(2))), ('rect_0', 'rect_1')), Integer(0))"
-
 def test_dict():
     from sympy import srepr
     from sympy.abc import x, y, z
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index d05052fe0c8c..bc185de760d1 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -936,11 +936,12 @@ def test_issue_14567():
 
 def test_diffgeom():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
+    x,y = symbols('x y', real=True)
     m = Manifold('M', 2)
     assert str(m) == "M"
     p = Patch('P', m)
     assert str(p) == "P"
-    rect = CoordSystem('rect', p)
+    rect = CoordSystem('rect', p, [x, y])
     assert str(rect) == "rect"
     b = BaseScalarField(rect, 0)
-    assert str(b) == "rect_0"
+    assert str(b) == "x"


Code review comments:
diff_hunk: @@ -321,85 +514,133 @@ def coord_tuple_transform_to(self, to_sys, coords):
 
     def jacobian(self, to_sys, coords):
         """Return the jacobian matrix of a transformation."""
+        SymPyDeprecationWarning(
path: sympy/diffgeom/diffgeom.py
text: @user1:
Although I acknowledge that the term `jacobian` was not used consistently across sympy as well as the general scientific community, I'm not sure this had been necessary to deprecate and introduce new methods

@author:
I'm convinced. Instead, I made `jacobian_matrix` as an alias of `jacobian`.

diff_hunk: @@ -120,134 +157,155 @@ def dim(self):
         return self.manifold.dim
 
 class CoordSystem(Atom):
-    """A coordinate system defined on the patch
+    """A coordinate system defined on the patch.
 
     Explanation
     ===========
 
-    This class contains all coordinate transformation logic.
+    Coordinate system is a system that uses one or more coordinates to uniquely determine
+    the position of the points or other geometric elements on a manifold [1].
+
+    By passing Symbols to *symbols* parameter, user can define the name and assumptions
+    of coordinate symbols of the coordinate system.
+    By passing *relations* parameter, user can define the tranform relations of coordinate
+    systems. Inverse transformation and indirect transformation can be found automatically.
+    If this parameter is not passed, coordinate transformation cannot be done.
 
     Parameters
     ==========
 
-    name : string
+    name : str
         The name of the coordinate system.
 
     patch : Patch
         The patch where the coordinate system is defined.
 
-    names : list of strings, optional
-        Determines how base scalar fields will be printed.
+    symbols : list of Symbols
+        Defines the names and assumptions of coordinate symbols.
+
+    relations : dict, optional
+        - key : tuple of two strings, who are the names of systems where
+            the coordinates transform from and transform to.
+        - value : Lambda returning the transformed coordinates.
 
     Examples
     ========
 
-    Define a Manifold and a Patch, and then define two coord systems on that
-    patch:
-
-    >>> from sympy import symbols, sin, cos, pi
+    >>> from sympy import symbols, pi, Lambda, Matrix, sqrt, atan2, cos, sin
     >>> from sympy.diffgeom import Manifold, Patch, CoordSystem
-    >>> from sympy.simplify import simplify
-    >>> r, theta = symbols('r, theta')
     >>> m = Manifold('M', 2)
-    >>> patch = Patch('P', m)
-    >>> rect = CoordSystem('rect', patch)
-    >>> polar = CoordSystem('polar', patch)
-    >>> rect in patch.coord_systems
-    True
+    >>> p = Patch('P', m)
 
-    Connect the coordinate systems. An inverse transformation is automatically
-    found by ``solve`` when possible:
+    >>> x, y = symbols('x y', real=True)
+    >>> r, theta = symbols('r theta', nonnegative=True)
+    >>> relation_dict = {
+    ... ('Car2D', 'Pol'): Lambda((x, y), Matrix([sqrt(x**2 + y**2), atan2(y, x)])),
+    ... ('Pol', 'Car2D'): Lambda((r, theta), Matrix([r*cos(theta), r*sin(theta)]))
+    ... }
+    >>> Car2D = CoordSystem('Car2D', p, [x, y], relation_dict)
+    >>> Pol = CoordSystem('Pol', p, [r, theta], relation_dict)
+
+    >>> Car2D.name
+    'Car2D'
+    >>> Car2D.dim
+    2
+    >>> Car2D.symbols
+    [x, y]
+    >>> Car2D.transformation(Pol)
+    Lambda((x, y), Matrix([
+    [sqrt(x**2 + y**2)],
+    [      atan2(y, x)]]))
 
-    >>> polar.connect_to(rect, [r, theta], [r*cos(theta), r*sin(theta)])
-    >>> polar.coord_tuple_transform_to(rect, [0, 2])
+    >>> Car2D.transform(Pol)
     Matrix([
-    [0],
-    [0]])
-    >>> polar.coord_tuple_transform_to(rect, [2, pi/2])
-    Matrix([
-    [0],
-    [2]])
-    >>> rect.coord_tuple_transform_to(polar, [1, 1]).applyfunc(simplify)
+    [sqrt(x**2 + y**2)],
+    [      atan2(y, x)]])
+    >>> Car2D.transform(Pol, [1, 2])
     Matrix([
-    [sqrt(2)],
-    [   pi/4]])
+    [sqrt(5)],
+    [atan(2)]])
 
-    Calculate the jacobian of the polar to cartesian transformation:
-
-    >>> polar.jacobian(rect, [r, theta])
+    >>> Pol.jacobian_matrix(Car2D)
     Matrix([
     [cos(theta), -r*sin(theta)],
     [sin(theta),  r*cos(theta)]])
-
-    Define a point using coordinates in one of the coordinate systems:
-
-    >>> p = polar.point([1, 3*pi/4])
-    >>> rect.point_to_coords(p)
+    >>> Pol.jacobian_matrix(Car2D, [1, pi/2])
     Matrix([
-    [-sqrt(2)/2],
-    [ sqrt(2)/2]])
-
-    Define a basis scalar field (i.e. a coordinate function), that takes a
-    point and returns its coordinates. It is an instance of ``BaseScalarField``.
-
-    >>> rect.coord_function(0)(p)
-    -sqrt(2)/2
-    >>> rect.coord_function(1)(p)
-    sqrt(2)/2
+    [0, -1],
+    [1,  0]])
 
-    Define a basis vector field (i.e. a unit vector field along the coordinate
-    line). Vectors are also differential operators on scalar fields. It is an
-    instance of ``BaseVectorField``.
-
-    >>> v_x = rect.base_vector(0)
-    >>> x = rect.coord_function(0)
-    >>> v_x(x)
-    1
-    >>> v_x(v_x(x))
-    0
-
-    Define a basis oneform field:
-
-    >>> dx = rect.base_oneform(0)
-    >>> dx(v_x)
-    1
-
-    If you provide a list of names the fields will print nicely:
-    - without provided names:
-
-    >>> x, v_x, dx
-    (rect_0, e_rect_0, drect_0)
-
-    - with provided names
+    References
+    ==========
 
-    >>> rect = CoordSystem('rect', patch, ['x', 'y'])
-    >>> rect.coord_function(0), rect.base_vector(0), rect.base_oneform(0)
-    (x, e_x, dx)
+    .. [1] https://en.wikipedia.org/wiki/Coordinate_system
 
     """
-    #  Contains a reference to the parent patch in order to be able to access
-    # other coordinate system charts.
-    def __new__(cls, name, patch, names=None):
+    def __new__(cls, name, patch, symbols=None, relations={}, **kwargs):
         if not isinstance(name, Str):
             name = Str(name)
-        # names is not in args because it is related only to printing, not to
-        # identifying the CoordSystem instance.
-        obj = super().__new__(cls, name, patch)
-        if not names:
-            names = ['%s_%d' % (name, i) for i in range(patch.dim)]
-        obj._names = tuple(str(i) for i in names)
-        obj.patch.coord_systems.append(obj)
-        obj.transforms = {}
-        # All the coordinate transformation logic is in this dictionary in the
-        # form of:
-        #  key = other coordinate system
-        #  value = tuple of  # TODO make these Lambda instances
-        #          - list of `Dummy` coordinates in this coordinate system
-        #          - list of expressions as a function of the Dummies giving
-        #          the coordinates in another coordinate system
-        obj._dummies = [Dummy(str(n)) for n in names]
+
+        # canonicallize the symbols
+        if symbols is None or any(not isinstance(s, Symbol) for s in symbols) :
+            # support deprecated signature
+            SymPyDeprecationWarning(
+                feature="Class signature 'names' of CoordSystem",
+                useinstead="class signature 'symbols'",
+                issue=19321,
+                deprecated_since_version="1.7"
+            ).warn()
+
+            if symbols is None:
+                names = kwargs.get('names', None)
+            else:
+                names = symbols
+
+            if names is None:
+                symbols = Tuple(*[Symbol('%s_%d' % (name, i), real=True) for i in range(patch.dim)])
+            else:
+                symbols = Tuple(*[Symbol(n, real=True) for n in names])
+        else:
+            symbols = Tuple(*[Symbol(s.name, **s._assumptions.generator) for s in symbols])
+
+        # canonicallize the relations
+        rel_temp = {}
+        for k,v in relations.items():
+            s1, s2 = k
+            if not isinstance(s1, Str):
+                s1 = Str(s1)
+            if not isinstance(s2, Str):
+                s2 = Str(s2)
+            key = Tuple(s1, s2)
+            rel_temp[key] = v
+        relations = Dict(rel_temp)
+
+        # construct the object
+        obj = super().__new__(cls, name, patch, symbols, relations)
+
+        # Add deprecated attributes
+        obj.transforms = _deprecated_dict(
+            "Mutable CoordSystem.transforms",
+            "'relations' parameter in class signature",
+            19321,
+            "1.7",
+            {}
+        )
+        obj._names =_deprecated_list(
+                "CoordSystem._names",
+                "CoordSystem.symbols",
+                19321,
+                "1.7",
+                [str(n) for n in symbols]
+            )
+        obj.patch.coord_systems.append(obj) # deprecated
+        obj._dummies = _deprecated_list(
path: sympy/diffgeom/diffgeom.py
text: @user1:
These are private methods, so it can be overwhelming to raise deprecation warnings for `_names` and `_dummies`

@author:
I guess you're right. I un-deprecated them.

diff_hunk: @@ -2506,31 +2506,52 @@ def _print_MatrixHomomorphism(self, h):
             self._print(h.domain), self._print(h.codomain))
 
     def _print_Manifold(self, manifold):
-        return r'\text{%s}' % manifold.name
+        string = manifold.name.name
+        if '{' in string:
+            name, supers, subs = string, [], []
+        else:
+            name, supers, subs = split_super_sub(string)
+
+            name = translate(name)
+            supers = [translate(sup) for sup in supers]
+            subs = [translate(sub) for sub in subs]
+
+        name = r'\text{%s}' % name
+        if supers:
+            name += "^{%s}" % " ".join(supers)
+        if subs:
+            name += "_{%s}" % " ".join(subs)
+
+        return name
 
     def _print_Patch(self, patch):
-        return r'\text{%s}_{\text{%s}}' % (patch.name, patch.manifold.name)
+        return r'\text{%s}_{%s}' % (self._print(patch.name), self._print(patch.manifold))
+
+    def _print_CoordSystem(self, coordsys):
+        return r'\text{%s}^{\text{%s}}_{%s}' % (
+            self._print(coordsys.name), self._print(coordsys.patch.name), self._print(coordsys.manifold)
+        )
 
-    def _print_CoordSystem(self, coords):
-        return r'\text{%s}^{\text{%s}}_{\text{%s}}' % (
-            coords.name, coords.patch.name, coords.patch.manifold.name
+    def _print_Point(self, point):
+        return '{%s}_{%s}' % (
+            self._print(point.coordinates()), self._print(point._coord_sys)
path: sympy/printing/latex.py
text: @user1:
I think that `_print_Point` should not be defined because it conflicts with printing of `Point` of `sympy.geometry` or `sympy.vector`. Maybe it was a bad idea to have duplicate names for classes at the first place, but it works better if you just leave it undefined and have some reasonable tuple notation of args.

@user1:
I have thought as a same idea to rename the Point in diffgeom. But it may be noted in the release notes.

@author:
Deleted this. I think we can give different name to class and introduce aliasing, but that would better be handled in another PR.