Repository: sympy/sympy

Issue:
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.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v158). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* diffgeom
  - `Manifold.patches` is deleted. ([#19368](https://github.com/sympy/sympy/pull/19368) by [@JSS95](https://github.com/JSS95) and [@mcpl-sympy](https://github.com/mcpl-sympy))

  - `Patch.coord_systems` is deleted. ([#19368](https://github.com/sympy/sympy/pull/19368) by [@JSS95](https://github.com/JSS95) and [@mcpl-sympy](https://github.com/mcpl-sympy))

  - `CoordSystem.transforms` and `CoordSystem.connect_to` are deprecated. Pass coordinate transform relation into the construction signal instead. ([#19368](https://github.com/sympy/sympy/pull/19368) by [@JSS95](https://github.com/JSS95) and [@mcpl-sympy](https://github.com/mcpl-sympy))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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

    #### Brief description of what is fixed or changed
    `Manifold.patches` is deleted.  
    `Patch.coord_systems` is deleted.  
    `CoordSystem.transforms` and `CoordSystem.connect_to` are deprecated. Pass coordinate transform relation into the construction signal instead.

    #### Other comments
    Before this PR, `diffgeom` module's `Manifold` and `Patch` had mutable attributes to register the other instances. This PR deletes that attributes; user will need to use external containers to register these.  
    Also, `CoordSystem.transforms` is deprecated. This `dict` used to store the transform relation with other coordinate systems. Since this relation is crucial for the module, it is better to make it immutable so that the information is included in `srepr` printing and `_hashable_content`. Now, user need to pass the transform relation to class construction signature, rather than using `connect_to` method after instance construction.

    #### 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` is deleted.  
        - `Patch.coord_systems` is deleted.  
        - `CoordSystem.transforms` and `CoordSystem.connect_to` are deprecated. Pass coordinate transform relation into the construction signal instead.

    <!-- END RELEASE NOTES -->

</details><p>

mcpl-sympy: @sylee957 Please review this.
codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/19368?src=pr&el=h1) Report
> Merging [#19368](https://codecov.io/gh/sympy/sympy/pull/19368?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/360755780a8074896eed4b5a6e4a9840be0ce5f0&el=desc) will **decrease** coverage by `0.002%`.
> The diff coverage is `64.444%`.

```diff
@@              Coverage Diff              @@
##            master    #19368       +/-   ##
=============================================
- Coverage   75.618%   75.616%   -0.003%     
=============================================
  Files          651       652        +1     
  Lines       169418    169580      +162     
  Branches     39973     40015       +42     
=============================================
+ Hits        128112    128230      +118     
- Misses       35703     35737       +34     
- Partials      5603      5613       +10     
```

sylee957: Regarding the deprecation issues, I think that `Manifold.patches` and `Patch.coord_systems` can raise some deprecation warnings by using some custom objects that wraps on `list` and `dict` that raises some warnings on item access or iterations.
But if the original mutable designs have other critical issues like memory leak, this can be ignored.
mcpl-sympy: @sylee957  
Sounds alright. I'll try that one.
mcpl-sympy: @sylee957

I added deprecation wrapper for list and dict. Please review this
sylee957: I also think that the future about #19318 should be addressed here because `fill_in_gaps` will not be valid after `connect_to` is deprecated.
I'd have an idea as implementing that as a external function that takes multiple `dict` instance and can figure out a composition, but I'm not sure if you may say that this is not sufficient for your application.
mcpl-sympy: @sylee957 
Yes, I was planning for that, either.  
My idea was like this: If we have explicit relation between A-B and B-C, then A-C are implicitly related as well. It's just that deriving the relation will be time-consuming because `solve` will be involved.  
So, instead of developing `fill_in_gaps`, I was going to make a method which solves and returns the relation between A and C then caches it. (Perhaps by using `cachit` decorator).
mcpl-sympy: I originally wrote this comment [here](https://github.com/mcpl-sympy/sympy/pull/2), but I'll quote it in this PR so that other people can see and share their opinions.

> I'm sorry, but I've been thinking of this awhile and I have to confess; I'm still not convinced that transform relation should be kept external.
> 
> For [curvilinear coordinates](https://en.wikipedia.org/wiki/Curvilinear_coordinates), relation with its 'parent' Cartesian coordinate system is the most important attribute; it literally defines the very nature of the system. Not only its Jacobian matrix, but its differential operator, metric tensor, everything depends on this relation.
> 
> I just can't convince myself that a single coordinate system is indetermined and can be Cartesian, curvilinear, affine, or anything else depending on the external dictionary. For me, this is tantamount to making the `args` of `Add`, `Mul`, or any other class undetermined and letting it dependent to external container - making it effectively mutable.
> 
> If we want to make a directed graph between the classes, we might as well make an algorithm to 'crawl' every connected system by querying whether one system 'knows' the other. If this is unacceptable because of performance reason, then it would be better to just make `CoordSystem` the subclass of `object` and let it be freely mutable.
sylee957: > For curvilinear coordinates, relation with its 'parent' Cartesian coordinate system is the most important attribute; it literally defines the very nature of the system.

I don't get what is such 'nature'. We can think as a whole collection of `(CoordSys1, CoordSys2, ..., Transforms)` to be definition of the system.

> Not only its Jacobian matrix, but its differential operator, metric tensor, everything depends on this relation.

As far as I see how `diffgeom` offers for now, all other stuff are derived from the jacobian matrix. I'm not sure that if are some more advanced concepts that cannot be derived from it, But again I believe that anything can be wrapped up in an another object that defines the concept, if the stuff have to be mutually defined.

> For me, this is tantamount to making the args of Add, Mul, or any other class undetermined and letting it dependent to external container

Addition and multiplication is not well defined from the context as elements of sets, but from the context of rings or fields. 
And it is possible to let such algebra as an optional context.

> making it effectively mutable.

I don't think that it should be sympy's scope to have some incremental registry because I have seen that it had been becoming problems.
But I left is as user's choice to use the mutable context or not. If you just pass on mutable dictionary to `transforms` keyword argument, it would suffice.

> If we want to make a directed graph between the classes, we might as well make an algorithm to 'crawl' every connected system by querying whether one system 'knows' the other. 

Why not crawl the graph defined in transformation dict?

> then it would be better to just make CoordSystem the subclass of object and let it be freely mutable.

I believe that my approach supersets any designs you suggest. Basically, if you pass some part of dictionary, it should be equivalent to wrapping the dictionary with object, and if you pass some mutable dictionary, it should be equivalent to defining mutable object.
mcpl-sympy: @sylee957  
I appreciate your comment. Thanks to it, I have a new perspective for this matter now.

Maintaining your design, I've come up with an idea of introducing `CoordSystems`class (its name may change), which takes coordinate systems and transformation dictionary as arguments.

I'll work on this and show you when it's done. Thanks again!
mcpl-sympy: I've merged your PR, but please don't merge this yet.
mcpl-sympy: @sylee957 
OK, it's done. Some changes have been made after your work, though.

1. I introduced `CoordSystemStructure` class. This is a subclass of `Dict`, but its `__getitem__` method is overrided so that inverse transform can be derived. After this PR, issue #19318 will be dealt with by adding additional method to this class.

2. Rather than including `CoordSystemStructure` in method signature, I added it in class signatures of `Point`, `BaseScalarField`, etc. With this, these classes now represent "geometric object defined with certain context". This structure exists in `args`, so one can introduce another coordinate system - by defining new `CoordSystemStructure` and replacing it into geometric entities.

Please review this!
mcpl-sympy: @sylee957, I added a feature to resolve #19318. I will be much obliged if you kindly review this.
sylee957: I notice a lot of slowdowns in `test_function_diffgeom_book` when I set the environment variable `SYMPY_USE_CACHE=no`, 
mcpl-sympy: @sylee957 
I tried some past commits, and here are the times took in `test_function_diffgeom_book` with cache off.

|HEAD|time (s)|
|---|---|
|master|9.7|
|87e49fd|9.3|
|1783523|12.6|
|c24965d (current)| 13.6|

Any suggestion, please?
sylee957: master, cache - 1.17 seconds
master, nocache - 11.09 seconds
head, cache - 1.74 seconds
head, nocache - 15.79 seconds

Okay, it could already been 10 times slower in the master. But I want to figure out which whether any changes here are introducing some slowdowns.
mcpl-sympy: @sylee957 
The first slowdown is encountered in 1783523, so I would like to start the investigation there.
mcpl-sympy: In master, each test took

|function|time (s)|
|---|---|
|test_functional_diffgeom_ch2|0.26|
|test_functional_diffgeom_ch3|1.53|
|test_functional_diffgeom_ch4|4.73|
|test_functional_diffgeom_ch6|0.26|


In 1783523, each test took

|function|time (s)|
|---|---|
|test_functional_diffgeom_ch2|0.23|
|test_functional_diffgeom_ch3|2.08|
|test_functional_diffgeom_ch4|**6.37**|
|test_functional_diffgeom_ch6|0.33|

In HEAD, each test took

|function|time (s)|
|---|---|
|test_functional_diffgeom_ch2|0.24|
|test_functional_diffgeom_ch3|**2.17**|
|test_functional_diffgeom_ch4|**6.57**|
|test_functional_diffgeom_ch6|**1.41**|
mcpl-sympy: `df(R2.e_x)`, `df(R2.e_y)`, `R2.dr(circ)`, `R2.dtheta(circ)` were the most time-consuming in `test_functional_diffgeom_ch4`.
sylee957: I only recognize that some behavior of substitution changed in `CoordSystem.jacobian`. And it looks like `subs` is doing more work than before.
But I'm not sure this alone had made things slower.
sylee957: ```diff
$ git diff
diff --git a/sympy/diffgeom/diffgeom.py b/sympy/diffgeom/diffgeom.py
index aca1d1322c..51c82b828b 100644
--- a/sympy/diffgeom/diffgeom.py
+++ b/sympy/diffgeom/diffgeom.py
@@ -750,7 +750,7 @@ def __call__(self, scalar_field):
         d_result = scalar_field.subs(list(zip(base_scalars, d_funcs)))
         d_result = d_result.diff(d_var)
         # Second step: e_x(x) -> 1 and e_x(r) -> cos(atan2(x, y))
-        coords = [Dummy() for _ in range(self._coord_sys.dim)]
+        coords = self._coord_sys._dummies
         d_funcs_deriv = [f.diff(d_var) for f in d_funcs]
         d_funcs_deriv_sub = []
         for b in base_scalars:
```

I guess this just made things slower because it doesn't equal the local dummies used by `jacobian`, and subs is expensive operation.
So I'd dump the attempt to replace `_dummies` cache for now.
mcpl-sympy: @sylee957 Indeed it was! Reverting that decresed the time spent to 8.44 second.
sylee957: I had been reading *"Functional Differential Geometry", Sussman and Wisdom.*, which contains the philosophy behind the designs of `diffgeom` module,
But I approve that the book doesn't say anything about whether coordinate transformations should be mutable or not, and can only be a matter of implementation details.

I think that the doc directive for `CoordSystemStructure` should be added in `diffgeom.rst`, because it should be one of the core feature of the `diffgeom` module after the API changes.
Since `Introduction` section completely empty there, a basic tutorial can be added to how to define coordinate systems for beginner users, but this can be optional for now.
mcpl-sympy: @sylee957 
Sure thing. I will add right away.
mcpl-sympy: @sylee957
I think #19451 might as well be implemented and documented in this PR.
sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v160). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* 7be7b3db236ec5d2ab2891952448b1553e08659b:
  - `release/compare_tar_against_git.py`
  - `sympy/deprecated/tests/test_deprecated_imports.py`
  - `sympy/matrices/expressions/companion.py`
  - `sympy/matrices/expressions/tests/test_companion.py`

If these files were added/deleted on purpose, you can ignore this message.


Patch:
diff --git a/doc/src/modules/diffgeom.rst b/doc/src/modules/diffgeom.rst
index 02fc0416af..764d9f3141 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/power.py b/sympy/core/power.py
index 0a3b8b564c..fe24eb8585 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1488,7 +1488,7 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         #    g has order O(x**d) where d is strictly positive.
         # 2) Then b**e = (f**e)*((1 + g)**e).
         #    (1 + g)**e is computed using binomial series.
-        from sympy import im, I, ceiling, polygamma, limit, logcombine, EulerGamma, exp, nan, zoo, log, factorial, ff, PoleError, O, powdenest, Wild
+        from sympy import im, I, ceiling, polygamma, logcombine, EulerGamma, exp, nan, zoo, log, factorial, ff, PoleError, O, powdenest, Wild
         from itertools import product
         self = powdenest(self, force=True).trigsimp()
         b, e = self.as_base_exp()
@@ -1553,14 +1553,7 @@ def mul(d1, d2):
                     res[ex] = res.get(ex, S.Zero) + d1[e1]*d2[e2]
             return res
 
-        try:
-            _, d = g.leadterm(x)
-        except (ValueError, NotImplementedError):
-            if limit(g/x**maxpow, x, 0) == 0:
-                # g has higher order zero
-                return f**e + e*f**e*g  # first term of binomial series
-            else:
-                raise NotImplementedError()
+        _, d = g.leadterm(x)
         if not d.is_positive:
             g = (b - f).simplify()/f
             _, d = g.leadterm(x)
diff --git a/sympy/diffgeom/__init__.py b/sympy/diffgeom/__init__.py
index 140f005394..8846a99510 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 643b032399..1d2ab0b27a 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,154 @@ 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. If not passed, these symbols are
+    generated automatically and are assumed to be real valued.
+    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, optional
+        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(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(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:
+            names = kwargs.get('names', None)
+            if names is None:
+                symbols = Tuple(
+                    *[Symbol('%s_%s' % (name.name, i), real=True) for i in range(patch.dim)]
+                )
+            else:
+                SymPyDeprecationWarning(
+                    feature="Class signature 'names' of CoordSystem",
+                    useinstead="class signature 'symbols'",
+                    issue=19321,
+                    deprecated_since_version="1.7"
+                ).warn()
+                symbols = Tuple(
+                    *[Symbol(n, real=True) for n in names]
+                )
+        else:
+            syms = []
+            for s in symbols:
+                if isinstance(s, Symbol):
+                    syms.append(Symbol(s.name, **s._assumptions.generator))
+                elif isinstance(s, str):
+                    SymPyDeprecationWarning(
+                        feature="Passing str as coordinate symbol's name",
+                        useinstead="Symbol which contains the name and assumption for coordinate symbol",
+                        issue=19321,
+                        deprecated_since_version="1.7"
+                    ).warn()
+                    syms.append(Symbol(s, real=True))
+            symbols = Tuple(*syms)
+
+        # 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 = [str(n) for n in symbols]
+        obj.patch.coord_systems.append(obj) # deprecated
+        obj._dummies = [Dummy(str(n)) for n in symbols] # deprecated
         obj._dummy = Dummy()
+
         return obj
 
     @property
@@ -258,34 +315,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 +452,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,92 +462,162 @@ 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 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]
             coords = transf[1].subs(list(zip(transf[0], coords)))
         return coords
 
-    def jacobian(self, to_sys, coords):
-        """Return the jacobian matrix of a transformation."""
-        with_dummies = self.coord_tuple_transform_to(
-            to_sys, self._dummies).jacobian(self._dummies)
-        return with_dummies.subs(list(zip(self._dummies, coords)))
+    def jacobian(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
+    jacobian_matrix = jacobian
+
+    def jacobian_determinant(self, sys, coordinates=None):
+        """Return the jacobian determinant of a transformation."""
+        return self.jacobian(sys, coordinates).det()
+
 
     ##########################################################################
-    # 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 +625,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 +640,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 +718,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 +731,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 +756,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 +806,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 +826,6 @@ class BaseVectorField(Expr):
 
     Parameters
     ==========
-
     coord_sys : CoordSystem
 
     index : integer
@@ -588,29 +833,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 +853,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,7 +915,7 @@ 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:
@@ -686,6 +939,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 +949,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 +966,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 +977,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 +994,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
+
+    @property
+    def v1(self):
+        return self.args[0]
 
-    def __init__(self, v1, v2):
-        super().__init__()
-        self._args = (v1, v2)
-        self._v1 = v1
-        self._v2 = v2
+    @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 +1054,6 @@ class Differential(Expr):
 
     >>> Differential(dg)
     0
-
     """
 
     is_commutative = False
@@ -811,18 +1063,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 +1123,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 +1138,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 +1167,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 +1187,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 +1197,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 +1213,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 +1249,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 +1262,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 +1272,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 +1309,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 +1342,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 +1433,45 @@ 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]
+
+    @property
+    def christoffel(self):
+        return self.args[1]
 
     def __call__(self, field):
         vectors = list(self._wrt.atoms(BaseVectorField))
@@ -1175,6 +1486,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 +1498,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 +1506,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 +1578,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 +1607,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 +1624,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 +1674,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 +1717,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 +1763,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 +1805,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,10 +1813,12 @@ 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 = []
@@ -1520,6 +1846,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 +1878,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 +1887,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 +1911,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 +1920,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 +1938,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 +1962,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 +2003,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 +2017,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 +2035,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 1226a63f03..55c964ccbd 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/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 98a8387a12..21b26f61f4 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -1004,7 +1004,6 @@ class LambertW(Function):
 
     .. [1] https://en.wikipedia.org/wiki/Lambert_W_function
     """
-    _singularities = (-Pow(S.Exp1, -1, evaluate=False), S.ComplexInfinity)
 
     @classmethod
     def eval(cls, x, k=None):
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index c695b344fc..f17c042d1b 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -2506,31 +2506,47 @@ 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, coords):
-        return r'\text{%s}^{\text{%s}}_{\text{%s}}' % (
-            coords.name, coords.patch.name, coords.patch.manifold.name
+    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_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 5fe19d241f..c4c3630fd3 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/repr.py b/sympy/printing/repr.py
index d156c6f312..1d33861800 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 2639a9a5d5..16e5d8cf90 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/series/limits.py b/sympy/series/limits.py
index 3b7dd4bcf7..a30161c309 100644
--- a/sympy/series/limits.py
+++ b/sympy/series/limits.py
@@ -241,7 +241,7 @@ def remove_abs(expr):
                 newe = e.subs(z, z + z0)
             try:
                 coeff, ex = newe.leadterm(z, cdir)
-            except (ValueError, NotImplementedError):
+            except ValueError:
                 pass
             else:
                 if ex > 0:
diff --git a/sympy/tensor/tensor.py b/sympy/tensor/tensor.py
index 8b04f6151f..4bb9b680f2 100644
--- a/sympy/tensor/tensor.py
+++ b/sympy/tensor/tensor.py
@@ -919,7 +919,7 @@ class TensorIndexType(Basic):
     dummy_name : name of the head of dummy indices
     dim : dimension, it can be a symbol or an integer or ``None``
     eps_dim : dimension of the epsilon tensor
-    metric_symmetry : integer that denotes metric symmetry or ``None`` for no metirc
+    metric_symmetry : integer that denotes metric symmetry or `None` for no metirc
     metric_name : string with the name of the metric tensor
 
     Attributes
@@ -933,15 +933,15 @@ class TensorIndexType(Basic):
     Notes
     =====
 
-    The possible values of the ``metric_symmetry`` parameter are:
+    The possible values of the `metric_symmetry` parameter are:
 
         ``1``   :   metric tensor is fully symmetric
         ``0``   :   metric tensor possesses no index symmetry
         ``-1``  :   metric tensor is fully antisymmetric
-        ``None``:   there is no metric tensor (metric equals to ``None``)
+        ``None``:   there is no metric tensor (metric equals to `None`)
 
     The metric is assumed to be symmetric by default. It can also be set
-    to a custom tensor by the ``.set_metric()`` method.
+    to a custom tensor by the `.set_metric()` method.
 
     If there is a metric the metric is used to raise and lower indices.
 


Test patch:
diff --git a/sympy/core/tests/test_args.py b/sympy/core/tests/test_args.py
index 6bdc452c62..998af7ea47 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():
@@ -4421,37 +4421,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 +4464,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 +4479,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 +4488,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/core/tests/test_power.py b/sympy/core/tests/test_power.py
index 6aa6ef2e57..f07107f5fa 100644
--- a/sympy/core/tests/test_power.py
+++ b/sympy/core/tests/test_power.py
@@ -286,7 +286,6 @@ def test_nseries():
     assert cbrt(I*x - 1)._eval_nseries(x, 4, None, -1) == (-1)**(S(1)/3)*exp(-2*I*pi/3) - \
     (-1)**(S(5)/6)*x*exp(-2*I*pi/3)/3 + (-1)**(S(1)/3)*x**2*exp(-2*I*pi/3)/9 + \
     5*(-1)**(S(5)/6)*x**3*exp(-2*I*pi/3)/81 + O(x**4)
-    assert (1 / (exp(-1/x) + 1/x))._eval_nseries(x, 2, None) == -x**2*exp(-1/x) + x
 
 
 def test_issue_6100_12942_4473():
diff --git a/sympy/diffgeom/tests/test_class_structure.py b/sympy/diffgeom/tests/test_class_structure.py
index 519eed5bfa..920180be66 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,8 @@ 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'])
+        assert cs_wname == cs_wname.func(*cs_wname.args)
diff --git a/sympy/diffgeom/tests/test_diffgeom.py b/sympy/diffgeom/tests/test_diffgeom.py
index 7c63cc67e5..0c9cc397cb 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/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index a7cdbd79e3..742fb34b81 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/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index 9c2d37e975..675ac0ab76 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 ed892c1801..481547cd26 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 d05052fe0c..bc185de760 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"
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 1169c1102f..733844e7be 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -814,12 +814,6 @@ def test_issue_18508():
     assert limit(sin(x)/sqrt(1-cos(x)), x, 0, dir='-') == -sqrt(2)
 
 
-def test_issue_18969():
-    a, b = symbols('a b', positive=True)
-    assert limit(LambertW(a), a, b) == LambertW(b)
-    assert limit(exp(LambertW(a)), a, b) == exp(LambertW(b))
-
-
 def test_issue_18992():
     assert limit(n/(factorial(n)**(1/n)), n, oo) == exp(1)