Repository: einsteinpy/einsteinpy

Issue:
Change the internal coordinate handling functions in metric modules
🐞 **Problem**
`.py` files in `metric` still rely on `utils` for coordinate handling.

🎯 **Goal**
  - Change them to newly created `coordinates` module functions

📋  **Steps to solve the problem**

 * Comment below about what you've started working on.
 * Add, commit, push your changes
 * Submit a pull request and add this in comments - `Addresses #<put issue number here>`
 * Ask for a review in comments section of pull request
 * Celebrate your contribution to this project 🎉


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=h1) Report
> Merging [#203](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=desc) into [master](https://codecov.io/gh/einsteinpy/einsteinpy/commit/20a2979b4601d755b0f7417b8f369bb61e39e309?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/graphs/tree.svg?width=650&token=H3oDVYA11f&height=150&src=pr)](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #203   +/-   ##
=======================================
  Coverage   92.78%   92.78%           
=======================================
  Files          22       22           
  Lines        1040     1040           
=======================================
  Hits          965      965           
  Misses         75       75
```


| [Impacted Files](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [src/einsteinpy/coordinates/velocity.py](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/diff?src=pr&el=tree#diff-c3JjL2VpbnN0ZWlucHkvY29vcmRpbmF0ZXMvdmVsb2NpdHkucHk=) | `100% <ø> (ø)` | :arrow_up: |
| [src/einsteinpy/coordinates/core.py](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/diff?src=pr&el=tree#diff-c3JjL2VpbnN0ZWlucHkvY29vcmRpbmF0ZXMvY29yZS5weQ==) | `94.59% <ø> (ø)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=footer). Last update [20a2979...54f4888](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

codecov[bot]: # [Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=h1) Report
> Merging [#203](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=desc) into [master](https://codecov.io/gh/einsteinpy/einsteinpy/commit/20a2979b4601d755b0f7417b8f369bb61e39e309?src=pr&el=desc) will **not change** coverage.
> The diff coverage is `n/a`.

[![Impacted file tree graph](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/graphs/tree.svg?width=650&token=H3oDVYA11f&height=150&src=pr)](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #203   +/-   ##
=======================================
  Coverage   92.78%   92.78%           
=======================================
  Files          22       22           
  Lines        1040     1040           
=======================================
  Hits          965      965           
  Misses         75       75
```


| [Impacted Files](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [src/einsteinpy/coordinates/velocity.py](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/diff?src=pr&el=tree#diff-c3JjL2VpbnN0ZWlucHkvY29vcmRpbmF0ZXMvdmVsb2NpdHkucHk=) | `100% <ø> (ø)` | :arrow_up: |
| [src/einsteinpy/coordinates/core.py](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203/diff?src=pr&el=tree#diff-c3JjL2VpbnN0ZWlucHkvY29vcmRpbmF0ZXMvY29yZS5weQ==) | `94.59% <ø> (ø)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=footer). Last update [20a2979...54f4888](https://codecov.io/gh/einsteinpy/einsteinpy/pull/203?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

pep8speaks: Hello @ritzvik! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`src/einsteinpy/coordinates/core.py`](https://github.com/einsteinpy/einsteinpy/blob/dcce1149dd54a0c53f52627bf219863b02f7a1ae/src/einsteinpy/coordinates/core.py):

> [Line 212:1](https://github.com/einsteinpy/einsteinpy/blob/dcce1149dd54a0c53f52627bf219863b02f7a1ae/src/einsteinpy/coordinates/core.py#L212): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace

* In the file [`src/einsteinpy/coordinates/velocity.py`](https://github.com/einsteinpy/einsteinpy/blob/dcce1149dd54a0c53f52627bf219863b02f7a1ae/src/einsteinpy/coordinates/velocity.py):

> [Line 262:1](https://github.com/einsteinpy/einsteinpy/blob/dcce1149dd54a0c53f52627bf219863b02f7a1ae/src/einsteinpy/coordinates/velocity.py#L262): [W293](https://duckduckgo.com/?q=pep8%20W293) blank line contains whitespace


ritzvik: @shreyasbapat I refactored the schwarzschild module. See whether it is satisfactory or requires changes, as I would the same in Kerr & Kerr-Newman also.
ritzvik: Would fix the other modules in subsequent PRs. Merge this if you feel like doing so. Rebasing it!
ritzvik: @shreyasbapat We should not make exceptions in CodeClimate. What do you think? Maybe open an issue to tackle the code smells.

Patch:
diff --git a/src/einsteinpy/coordinates/core.py b/src/einsteinpy/coordinates/core.py
index a07d2b4..320acf5 100644
--- a/src/einsteinpy/coordinates/core.py
+++ b/src/einsteinpy/coordinates/core.py
@@ -14,9 +14,10 @@ class Cartesian:
 
         Parameters
         ----------
-        x : ~astropy.units
-        y : ~astropy.units
-        z : ~astropy.units
+        x : ~astropy.units.quantity.Quantity
+        y : ~astropy.units.quantity.Quantity
+        z : ~astropy.units.quantity.Qauntity
+
         """
         self.x = x
         self.y = y
@@ -47,7 +48,7 @@ class Cartesian:
 
         Returns
         -------
-        ~astropy.units
+        ~astropy.units.quantity.Quantity
             Euclidean norm with units.
 
         """
@@ -63,7 +64,7 @@ class Cartesian:
 
         Returns
         -------
-        ~astropy.units
+        ~astropy.units.quantity.Quantity
             Dot product with units
 
         """
@@ -97,7 +98,7 @@ class Cartesian:
 
         Parameters
         ----------
-        a : astropy.units
+        a : ~astropy.units.quantity.Quantity
             a = J/Mc , the angular momentum per unit mass of the black hole per speed of light.
 
         Returns
@@ -126,9 +127,10 @@ class Spherical:
 
         Parameters
         ----------
-        r : ~astropy.units
-        theta : ~astropy.units
-        phi : ~astropy.units
+        r : ~astropy.units.quantity.Quantity
+        theta : ~astropy.units.quantity.Quantity
+        phi : ~astropy.units.quantity.Quantity
+
         """
         self.r = r
         self.theta = theta
@@ -178,7 +180,7 @@ class Spherical:
 
         Parameters
         ----------
-        a : astropy.units
+        a : ~astropy.units.quantity.Quantity
             a = J/Mc , the angular momentum per unit mass of the black hole per speed of light.
 
         Returns
@@ -203,10 +205,11 @@ class BoyerLindquist:
 
         Parameters
         ----------
-        r : ~astropy.units
-        theta : ~astropy.units
-        phi : ~astropy.units
-        a : ~astropy.units
+        r : ~astropy.units.quantity.Quantity
+        theta : ~astropy.units.quantity.Quantity
+        phi : ~astropy.units.quantity.Quantity
+        a : ~astropy.units.quantity.Quantity
+        
         """
         self.r = r
         self.theta = theta
diff --git a/src/einsteinpy/coordinates/velocity.py b/src/einsteinpy/coordinates/velocity.py
index f9f939f..872e446 100644
--- a/src/einsteinpy/coordinates/velocity.py
+++ b/src/einsteinpy/coordinates/velocity.py
@@ -18,12 +18,13 @@ class CartesianDifferential(Cartesian):
 
         Parameters
         ----------
-        x : ~astropy.units
-        y : ~astropy.units
-        z : ~astropy.units
-        v_x : ~astropy.units
-        v_y : ~astropy.units
-        v_z : ~astropy.units
+        x : ~astropy.units.quantity.Quantity
+        y : ~astropy.units.quantity.Quantity
+        z : ~astropy.units.quantity.Quantity
+        v_x : ~astropy.units.quantity.Quantity
+        v_y : ~astropy.units.quantity.Quantity
+        v_z : ~astropy.units.quantity.Quantity
+
         """
         super(CartesianDifferential, self).__init__(x, y, z)
         self.v_x = v_x
@@ -66,7 +67,7 @@ class CartesianDifferential(Cartesian):
 
         Returns
         -------
-        SphericalDifferential : ~einsteinpy.coordinates.velocity.SphericalDifferential
+        ~einsteinpy.coordinates.velocity.SphericalDifferential
             Spherical representation of the velocity in Cartesian Coordinates.
 
         """
@@ -93,12 +94,12 @@ class CartesianDifferential(Cartesian):
 
         Parameters
         ----------
-        a : astropy.units
+        a : ~astropy.units.quantity.Quantity
             a = J/Mc , the angular momentum per unit mass of the black hole per speed of light.
 
         Returns
         -------
-        BoyerLindquistDifferential : ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
+        ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
             Boyer-Lindquist representation of the velocity in Cartesian Coordinates.
 
         """
@@ -141,12 +142,13 @@ class SphericalDifferential(Spherical):
 
         Parameters
         ----------
-        r : ~astropy.units
-        theta : ~astropy.units
-        phi : ~astropy.units
-        v_r : ~astropy.units
-        v_t : ~astropy.units
-        v_p : ~astropy.units
+        r : ~astropy.units.quantity.Quantity
+        theta : ~astropy.units.quantity.Quantity
+        phi : ~astropy.units.quantity.Quantity
+        v_r : ~astropy.units.quantity.Quantity
+        v_t : ~astropy.units.quantity.Quantity
+        v_p : ~astropy.units.quantity.Quantity
+
         """
         super(SphericalDifferential, self).__init__(r, theta, phi)
         self.v_r = v_r
@@ -189,7 +191,7 @@ class SphericalDifferential(Spherical):
 
         Returns
         -------
-        CartesianDifferential : ~einsteinpy.coordinates.velocity.CartesianDifferential
+        ~einsteinpy.coordinates.velocity.CartesianDifferential
             Cartesian representation of the velocity in Spherical Coordinates.
 
         """
@@ -217,12 +219,12 @@ class SphericalDifferential(Spherical):
 
         Parameters
         ----------
-        a : astropy.units
+        a : ~astropy.units.quantity.Quantity
             a = J/Mc , the angular momentum per unit mass of the black hole per speed of light.
 
         Returns
         -------
-        BoyerLindquistDifferential : ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
+        ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
             Boyer-Lindquist representation of the velocity in Spherical Coordinates.
 
         """
@@ -250,13 +252,14 @@ class BoyerLindquistDifferential(BoyerLindquist):
 
         Parameters
         ----------
-        r : ~astropy.units
-        theta : ~astropy.units
-        phi : ~astropy.units
-        v_r : ~astropy.units
-        v_t : ~astropy.units
-        v_p : ~astropy.units
-        a : ~astropy.units
+        r : ~astropy.units.quantity.Quantity
+        theta : ~astropy.units.quantity.Quantity
+        phi : ~astropy.units.quantity.Quantity
+        v_r : ~astropy.units.quantity.Quantity
+        v_t : ~astropy.units.quantity.Quantity
+        v_p : ~astropy.units.quantity.Quantity
+        a : ~astropy.units.quantity.Quantity
+
         """
         super(BoyerLindquistDifferential, self).__init__(r, theta, phi, a)
         self.v_r = v_r
@@ -303,7 +306,7 @@ class BoyerLindquistDifferential(BoyerLindquist):
 
         Returns
         -------
-        CartesianDifferential : ~einsteinpy.coordinates.velocity.CartesianDifferential
+        ~einsteinpy.coordinates.velocity.CartesianDifferential
             Cartesian representation of the velocity in Boyer-Lindquist Coordinates.
 
         """
@@ -330,7 +333,7 @@ class BoyerLindquistDifferential(BoyerLindquist):
 
         Returns
         -------
-        SphericalDifferential : ~einsteinpy.coordinates.velocity.SphericalDifferential
+        ~einsteinpy.coordinates.velocity.SphericalDifferential
             Spherical representation of the velocity in Boyer-Lindquist Coordinates.
 
         """
diff --git a/src/einsteinpy/metric/kerr.py b/src/einsteinpy/metric/kerr.py
index 6a81b84..6f98b0b 100644
--- a/src/einsteinpy/metric/kerr.py
+++ b/src/einsteinpy/metric/kerr.py
@@ -4,9 +4,9 @@ import astropy.units as u
 import numpy as np
 
 from einsteinpy import constant
+from einsteinpy.coordinates import BoyerLindquistDifferential, CartesianDifferential
 from einsteinpy.integrators import RK45
-from einsteinpy.utils import *
-from einsteinpy.utils import kerr_utils, schwarzschild_radius as scr
+from einsteinpy.utils import kerr_utils, schwarzschild_radius
 
 _G = constant.G.value
 _c = constant.c.value
@@ -18,104 +18,68 @@ class Kerr:
     """
 
     @u.quantity_input(time=u.s, M=u.kg)
-    def __init__(self, pos_vec, vel_vec, time, M, a):
+    def __init__(self, bl_coords, M, time):
+        self.input_coords = bl_coords
         self.M = M
-        self.a = a
-        self.pos_vec = pos_vec
-        self.vel_vec = vel_vec
+        self.a = self.input_coords.a.to(u.m)
         self.time = time
-        self.time_vel = kerr_utils.kerr_time_velocity(
-            self.pos_vec, self.vel_vec, self.M, self.a
-        )
-        self.initial_vec = np.hstack(
-            (self.time.value, self.pos_vec, self.time_vel.value, self.vel_vec)
+
+        pos_vec, vel_vec = (
+            self.input_coords.si_values()[:3],
+            self.input_coords.si_values()[3:],
         )
-        self.schwarzschild_r = scr(M)
+        time_vel = kerr_utils.kerr_time_velocity(pos_vec, vel_vec, self.M, self.a.value)
 
-    @classmethod
-    def _classmethod_handler(cls, pos_vec, vel_vec, time, M, a):
-        cls.units_list = [
-            u.s,
-            u.m,
-            u.rad,
-            u.rad,
-            u.one,
-            u.m / u.s,
-            u.rad / u.s,
-            u.rad / u.s,
-        ]
-        pos_vec_vals = [
-            pos_vec[i].to(cls.units_list[i + 1]).value for i in range(len(pos_vec))
-        ]
-        vel_vec_vals = [
-            vel_vec[i].to(cls.units_list[i + 5]).value for i in range(len(vel_vec))
-        ]
-        return cls(
-            np.array(pos_vec_vals), np.array(vel_vec_vals), time.to(u.s), M.to(u.kg), a
+        self.initial_vec = np.hstack(
+            (self.time.value, pos_vec, time_vel.value, vel_vec)
         )
+        self.scr = schwarzschild_radius(M)
 
     @classmethod
     @u.quantity_input(time=u.s, M=u.kg)
-    def from_BL(cls, pos_vec, vel_vec, time, M, a):
+    def from_BL(cls, coords, M, time=0 * u.s):
         """
         Constructor
         Initialize from Boyer-Lindquist Coordinates
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
+        coords : ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
+            Initial positions and velocities of particle in BL Coordinates, and spin factor of massive body.
+        M : ~astropy.units.quantity.Quantity
             Mass of the body
-        a : float
-            Spin factor of massive body
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
         cls.input_coord_system = "Boyer-Lindquist"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        return cls._classmethod_handler(pos_vec, vel_vec, time, M, a)
+        return cls(coords, M, time)
 
     @classmethod
-    def from_cartesian(cls, pos_vec, vel_vec, time, M, a):
+    @u.quantity_input(time=u.s, M=u.kg, a=u.m)
+    def from_cartesian(cls, coords, M, a, time=0 * u.s):
         """
         Constructor
         Initialize from Cartesian Coordinates
 
         Parameters
         ----------
-        pos_vec : list
-            list of x, y and z components along with ~astropy.units
-        vel_vec : list
-            list of velocities of x, y, and z components along with ~astropy.units
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
+        coords : ~einsteinpy.coordinates.velocity.CartesianDifferential
+            Object having both initial positions and velocities of particle in Cartesian Coordinates
+        M : ~astropy.units.quantity.Quantity
             Mass of the body
-        a : float
-            Spin factor of massive body
+        a : ~astropy.units.quantity.Quantity
+            Spin factor of the massive body. Angular momentum divided by mass divided by speed of light.
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
         cls.input_coord_system = "Cartesian"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        bl_pos_vec, bl_vel_vec = C2BL_units(pos_vec, vel_vec, a)
-        return cls._classmethod_handler(bl_pos_vec, bl_vel_vec, time, M, a)
+        bl_coords = coords.bl_differential(a)
+        return cls(bl_coords, M, time)
 
     def f_vec(self, ld, vec):
-        chl = kerr_utils.christoffels(vec[1], vec[2], self.M.value, self.a)
+        chl = kerr_utils.christoffels(vec[1], vec[2], self.M.value, self.a.value)
         vals = np.zeros(shape=(8,), dtype=float)
         for i in range(4):
             vals[i] = vec[i + 4]
@@ -163,9 +127,9 @@ class Kerr:
         Parameters
         ----------
         start_lambda : float
-            Starting lambda, defaults to 0.0, (lambda ~= t)
+            Starting lambda(proper time), defaults to 0, (lambda ~= t)
         end_lamdba : float
-            Lambda where iteartions will stop, defaults to 100000
+            Lambda(proper time) where iteartions will stop, defaults to 100000
         stop_on_singularity : bool
             Whether to stop further computation on reaching singularity, defaults to True
         OdeMethodKwargs : dict
@@ -176,12 +140,14 @@ class Kerr:
 
         Returns
         -------
-        tuple
-            (~numpy.array of lambda, (n,8) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        ~numpy.ndarray
+            N-element array containing proper time.
+        ~numpy.ndarray
+            (n,8) shape array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        vec_list = list()
-        lambda_list = list()
+        vecs = list()
+        lambdas = list()
         singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
@@ -190,40 +156,42 @@ class Kerr:
             t_bound=end_lambda,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
-        _event_hor = kerr_utils.event_horizon(self.M.value, self.a)[0] * 1.001
+
+        _event_hor = kerr_utils.event_horizon(self.M.value, self.a.value)[0] * 1.001
+
         while ODE.t < end_lambda:
-            vec_list.append(ODE.y)
-            lambda_list.append(ODE.t)
+            vecs.append(ODE.y)
+            lambdas.append(ODE.t)
             ODE.step()
             if (not singularity_reached) and (ODE.y[1] <= _event_hor):
-                warnings.warn(
-                    "r component of position vector reached event horizon. ",
-                    RuntimeWarning,
-                )
+                warnings.warn("particle reached event horizon. ", RuntimeWarning)
                 if stop_on_singularity:
                     break
                 else:
                     singularity_reached = True
 
-        def _not_cartesian():
-            return (np.array(lambda_list), np.array(vec_list))
-
-        def _cartesian():
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-            return (np.array(lambda_list), BL2C_8dim(np.array(vec_list), self.a))
+        vecs, lambdas = np.array(vecs), np.array(lambdas)
 
-        choice_dict = {0: _not_cartesian, 1: _cartesian}
-        return choice_dict[int(return_cartesian)]()
+        if not return_cartesian:
+            return lambdas, vecs
+        else:
+            cart_vecs = list()
+            for v in vecs:
+                si_vals = (
+                    BoyerLindquistDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
+                        self.a,
+                    )
+                    .cartesian_differential()
+                    .si_values()
+                )
+                cart_vecs.append(np.hstack((v[0], si_vals[:3], v[4], si_vals[3:])))
+            return lambdas, np.array(cart_vecs)
 
     def calculate_trajectory_iterator(
         self,
@@ -233,8 +201,8 @@ class Kerr:
         return_cartesian=False,
     ):
         """
-        Calculate trajectory in manifold according to geodesic equation
-        Yields an iterator
+        Calculate trajectory in manifold according to geodesic equation.
+        Yields an iterator.
 
         Parameters
         ----------
@@ -250,11 +218,12 @@ class Kerr:
 
         Yields
         ------
-        tuple
-            (lambda, (8,) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        float
+            proper time
+        ~numpy.ndarray
+            array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
             t0=start_lambda,
@@ -262,39 +231,32 @@ class Kerr:
             t_bound=1e300,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
-        _event_hor = kerr_utils.event_horizon(self.M.value, self.a)[0] * 1.001
+        singularity_reached = False
+        _event_hor = kerr_utils.event_horizon(self.M.value, self.a.value)[0] * 1.001
 
-        def yielder_func():
-            nonlocal singularity_reached
-            while True:
-                if not return_cartesian:
-                    yield (ODE.t, ODE.y)
-                else:
-                    temp = np.copy(ODE.y)
-                    temp[1:4] = BLToCartesian_pos(ODE.y[1:4], self.a)
-                    temp[5:8] = BLToCartesian_vel(ODE.y[1:4], ODE.y[5:8], self.a)
-                    yield (ODE.t, temp)
-                ODE.step()
-                if (not singularity_reached) and (ODE.y[1] <= _event_hor):
-                    warnings.warn(
-                        "r component of position vector reached event horizon. ",
-                        RuntimeWarning,
+        while True:
+            if not return_cartesian:
+                yield ODE.t, ODE.y
+            else:
+                v = ODE.y
+                si_vals = (
+                    BoyerLindquistDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
+                        self.a,
                     )
-                    if stop_on_singularity:
-                        break
-                    else:
-                        singularity_reached = True
-
-        if return_cartesian:
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-        return yielder_func()
+                    .cartesian_differential()
+                    .si_values()
+                )
+                yield ODE.t, np.hstack((v[0], si_vals[:3], v[4], si_vals[3:]))
+            ODE.step()
+            if (not singularity_reached) and (ODE.y[1] <= _event_hor):
+                warnings.warn("particle reached event horizon. ", RuntimeWarning)
+                if stop_on_singularity:
+                    break
+                else:
+                    singularity_reached = True
diff --git a/src/einsteinpy/metric/kerrnewman.py b/src/einsteinpy/metric/kerrnewman.py
index 2e041b9..8b774ea 100644
--- a/src/einsteinpy/metric/kerrnewman.py
+++ b/src/einsteinpy/metric/kerrnewman.py
@@ -4,8 +4,8 @@ import astropy.units as u
 import numpy as np
 
 from einsteinpy import constant
+from einsteinpy.coordinates import BoyerLindquistDifferential, CartesianDifferential
 from einsteinpy.integrators import RK45
-from einsteinpy.utils import *
 from einsteinpy.utils import kerrnewman_utils, schwarzschild_radius
 
 _G = constant.G.value
@@ -20,129 +20,86 @@ class KerrNewman:
     """
 
     @u.quantity_input(time=u.s, M=u.kg, Q=u.C)
-    def __init__(self, pos_vec, vel_vec, q, time, M, a, Q):
+    def __init__(self, bl_coords, q, M, Q, time):
+        self.input_coords = bl_coords
         self.M = M
-        self.a = a
+        self.a = self.input_coords.a.to(u.m)
         self.Q = Q
         self.q = q
-        self.pos_vec = pos_vec
-        self.vel_vec = vel_vec
+
+        pos_vec, vel_vec = (
+            self.input_coords.si_values()[:3],
+            self.input_coords.si_values()[3:],
+        )
         self.time = time
-        self.time_vel = kerrnewman_utils.kerrnewman_time_velocity(
-            self.pos_vec, self.vel_vec, self.M, self.a, Q
+        time_vel = kerrnewman_utils.kerrnewman_time_velocity(
+            pos_vec, vel_vec, self.M, self.a.value, Q
         )
         self.initial_vec = np.hstack(
-            (self.time.value, self.pos_vec, self.time_vel.value, self.vel_vec)
-        )
-        self.schwarzschild_r = schwarzschild_radius(M)
-
-    @classmethod
-    def _classmethod_handler(cls, pos_vec, vel_vec, q, time, M, a, Q):
-        cls.units_list = [
-            u.s,
-            u.m,
-            u.rad,
-            u.rad,
-            u.one,
-            u.m / u.s,
-            u.rad / u.s,
-            u.rad / u.s,
-        ]
-        pos_vec_vals = [
-            pos_vec[i].to(cls.units_list[i + 1]).value for i in range(len(pos_vec))
-        ]
-        vel_vec_vals = [
-            vel_vec[i].to(cls.units_list[i + 5]).value for i in range(len(vel_vec))
-        ]
-        return cls(
-            np.array(pos_vec_vals),
-            np.array(vel_vec_vals),
-            q.to(u.C / u.kg),
-            time.to(u.s),
-            M.to(u.kg),
-            a,
-            Q.to(u.C),
+            (self.time.value, pos_vec, time_vel.value, vel_vec)
         )
+        self.scr = schwarzschild_radius(M)
 
     @classmethod
     @u.quantity_input(q=u.C / u.kg, time=u.s, M=u.kg, Q=u.C)
-    def from_BL(cls, pos_vec, vel_vec, q, time, M, a, Q):
+    def from_BL(cls, coords, q, M, Q, time=0 * u.s):
         """
-        Constructor
-        Initialize from Boyer-Lindquist Coordinates
+        Constructor.
+        Initialize from Boyer-Lindquist Coordinates.
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units
-        q : ~astropy.units.C/astropy.units.kg
-            Charge per mass of test particle
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
-            Mass of the body
-        a : float
-            Spin factor of massive body
-        Q : ~astropy.units.C
+        coords : ~einsteinpy.coordinates.velocity.BoyerLindquistDifferential
+            Initial positions and velocities of particle in BL Coordinates, and spin factor of massive body.
+        q : ~astropy.units.quantity.Quantity
+            Charge per unit mass of test particle
+        M : ~astropy.units.quantity.Quantity
+            Mass of the massive body
+        Q : ~astropy.units.quantity.Quantity
             Charge on the massive body
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
         cls.input_coord_system = "Boyer-Lindquist"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        return cls._classmethod_handler(pos_vec, vel_vec, q, time, M, a, Q)
+        return cls(coords, q, M, Q, time)
 
     @classmethod
-    @u.quantity_input(q=u.C / u.kg, time=u.s, M=u.kg, Q=u.C)
-    def from_cartesian(cls, pos_vec, vel_vec, q, time, M, a, Q):
+    @u.quantity_input(q=u.C / u.kg, time=u.s, M=u.kg, Q=u.C, a=u.m)
+    def from_cartesian(cls, coords, q, M, a, Q, time=0 * u.s):
         """
-        Constructor
-        Initialize from Cartesian Coordinates
+        Constructor.
+        Initialize from Cartesian Coordinates.
 
         Parameters
         ----------
-        pos_vec : list
-            list of x, y and z components along with ~astropy.units
-        vel_vec : list
-            list of velocities of x, y, and z components along with ~astropy.units
-        q : ~astropy.units.C/astropy.units.kg
-            Charge per mass of test particle
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
-            Mass of the body
-        a : float
-            Spin factor of massive body
-        Q : ~astropy.units.C
+        coords : ~einsteinpy.coordinates.velocity.CartesianDifferential
+            Initial positions and velocities of particle in Cartesian Coordinates.
+        q : ~astropy.units.quantity.Quantity
+            Charge per unit mass of test particle
+        M : ~astropy.units.quantity.Quantity
+            Mass of the massive body
+        a : ~astropy.units.quantity.Quantity
+            Spin factor of the massive body(Angular Momentum per unit mass per speed of light)
+        Q : ~astropy.units.quantity.Quantity
             Charge on the massive body
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
         cls.input_coord_system = "Cartesian"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        bl_pos_vec, bl_vel_vec = C2BL_units(pos_vec, vel_vec, a)
-        return cls._classmethod_handler(bl_pos_vec, bl_vel_vec, q, time, M, a, Q)
+        bl_coords = coords.bl_differential(a)
+        return cls(bl_coords, q, M, Q, time)
 
     def f_vec(self, ld, vec):
-        _scr = self.schwarzschild_r.value
         chl = kerrnewman_utils.christoffels(
-            vec[1], vec[2], self.M.value, self.a, self.Q.value
+            vec[1], vec[2], self.M.value, self.a.value, self.Q.value
         )
         maxwell = kerrnewman_utils.maxwell_tensor_contravariant(
-            vec[1], vec[2], self.a, self.Q.value, self.M.value
+            vec[1], vec[2], self.a.value, self.Q.value, self.M.value
         )
         metric = kerrnewman_utils.metric(
-            vec[1], vec[2], self.M.value, self.a, self.Q.value
+            vec[1], vec[2], self.M.value, self.a.value, self.Q.value
         )
         vals = np.zeros(shape=(8,), dtype=float)
         for i in range(4):
@@ -194,9 +151,9 @@ class KerrNewman:
         Parameters
         ----------
         start_lambda : float
-            Starting lambda, defaults to 0.0, (lambda ~= t)
+            Starting lambda(proper time), defaults to 0.0, (lambda ~= t)
         end_lamdba : float
-            Lambda where iteartions will stop, defaults to 100000
+            Lambda(proper time) where iteartions will stop, defaults to 100000
         stop_on_singularity : bool
             Whether to stop further computation on reaching singularity, defaults to True
         OdeMethodKwargs : dict
@@ -207,12 +164,14 @@ class KerrNewman:
 
         Returns
         -------
-        tuple
-            (~numpy.array of lambda, (n,8) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        ~numpy.ndarray
+            N-element array containing proper time.
+        ~numpy.ndarray
+            (n,8) shape array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        vec_list = list()
-        lambda_list = list()
+        vecs = list()
+        lambdas = list()
         singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
@@ -221,39 +180,42 @@ class KerrNewman:
             t_bound=end_lambda,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
+
+        _scr = self.scr.value * 1.001
+
         while ODE.t < end_lambda:
-            vec_list.append(ODE.y)
-            lambda_list.append(ODE.t)
+            vecs.append(ODE.y)
+            lambdas.append(ODE.t)
             ODE.step()
             if (not singularity_reached) and (ODE.y[1] <= _scr):
-                warnings.warn(
-                    "r component of position vector reached Schwarzchild Radius. ",
-                    RuntimeWarning,
-                )
+                warnings.warn("particle reached Schwarzchild Radius. ", RuntimeWarning)
                 if stop_on_singularity:
                     break
                 else:
                     singularity_reached = True
 
-        def _not_cartesian():
-            return (np.array(lambda_list), np.array(vec_list))
-
-        def _cartesian():
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-            return (np.array(lambda_list), BL2C_8dim(np.array(vec_list), self.a))
+        vecs, lambdas = np.array(vecs), np.array(lambdas)
 
-        choice_dict = {0: _not_cartesian, 1: _cartesian}
-        return choice_dict[int(return_cartesian)]()
+        if not return_cartesian:
+            return lambdas, vecs
+        else:
+            cart_vecs = list()
+            for v in vecs:
+                si_vals = (
+                    BoyerLindquistDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
+                        self.a,
+                    )
+                    .cartesian_differential()
+                    .si_values()
+                )
+                cart_vecs.append(np.hstack((v[0], si_vals[:3], v[4], si_vals[3:])))
+            return lambdas, np.array(cart_vecs)
 
     def calculate_trajectory_iterator(
         self,
@@ -263,8 +225,8 @@ class KerrNewman:
         return_cartesian=False,
     ):
         """
-        Calculate trajectory in manifold according to geodesic equation
-        Yields an iterator
+        Calculate trajectory in manifold according to geodesic equation.
+        Yields an iterator.
 
         Parameters
         ----------
@@ -280,11 +242,12 @@ class KerrNewman:
 
         Yields
         ------
-        tuple
-            (lambda, (8,) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        float
+            proper time
+        ~numpy.ndarray
+            array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
             t0=start_lambda,
@@ -292,38 +255,32 @@ class KerrNewman:
             t_bound=1e300,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
+        singularity_reached = False
+        _scr = self.scr.value * 1.001
 
-        def yielder_func():
-            nonlocal singularity_reached
-            while True:
-                if not return_cartesian:
-                    yield (ODE.t, ODE.y)
-                else:
-                    temp = np.copy(ODE.y)
-                    temp[1:4] = BLToCartesian_pos(ODE.y[1:4], self.a)
-                    temp[5:8] = BLToCartesian_vel(ODE.y[1:4], ODE.y[5:8], self.a)
-                    yield (ODE.t, temp)
-                ODE.step()
-                if (not singularity_reached) and (ODE.y[1] <= _scr):
-                    warnings.warn(
-                        "r component of position vector reached Schwarzchild Radius. ",
-                        RuntimeWarning,
+        while True:
+            if not return_cartesian:
+                yield ODE.t, ODE.y
+            else:
+                v = ODE.y
+                si_vals = (
+                    BoyerLindquistDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
+                        self.a,
                     )
-                    if stop_on_singularity:
-                        break
-                    else:
-                        singularity_reached = True
-
-        if return_cartesian:
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-        return yielder_func()
+                    .cartesian_differential()
+                    .si_values()
+                )
+                yield ODE.t, np.hstack((v[0], si_vals[:3], v[4], si_vals[3:]))
+            ODE.step()
+            if (not singularity_reached) and (ODE.y[1] <= _scr):
+                warnings.warn("particle reached Schwarzschild Radius. ", RuntimeWarning)
+                if stop_on_singularity:
+                    break
+                else:
+                    singularity_reached = True
diff --git a/src/einsteinpy/metric/schwarzschild.py b/src/einsteinpy/metric/schwarzschild.py
index 7bac69b..fdcf16f 100644
--- a/src/einsteinpy/metric/schwarzschild.py
+++ b/src/einsteinpy/metric/schwarzschild.py
@@ -4,8 +4,8 @@ import astropy.units as u
 import numpy as np
 
 from einsteinpy import constant
+from einsteinpy.coordinates import CartesianDifferential, SphericalDifferential
 from einsteinpy.integrators import RK45
-from einsteinpy.utils import *
 from einsteinpy.utils import schwarzschild_radius, schwarzschild_utils
 
 _G = constant.G.value
@@ -18,93 +18,56 @@ class Schwarzschild:
     """
 
     @u.quantity_input(time=u.s, M=u.kg)
-    def __init__(self, pos_vec, vel_vec, time, M):
+    def __init__(self, sph_coords, M, time):
         self.M = M
-        self.pos_vec = pos_vec
-        self.vel_vec = vel_vec
+        self.input_coords = sph_coords
         self.time = time
-        self.time_vel = time_velocity(pos_vec, vel_vec, M)
-        self.initial_vec = np.hstack(
-            (self.time.value, self.pos_vec, self.time_vel.value, self.vel_vec)
+        pos_vec, vel_vec = (
+            self.input_coords.si_values()[:3],
+            self.input_coords.si_values()[3:],
         )
-        self.schwarzschild_r = schwarzschild_radius(M)
-
-    @classmethod
-    def _classmethod_handler(cls, pos_vec, vel_vec, time, M):
-        cls.units_list = [
-            u.s,
-            u.m,
-            u.rad,
-            u.rad,
-            u.one,
-            u.m / u.s,
-            u.rad / u.s,
-            u.rad / u.s,
-        ]
-        pos_vec_vals = [
-            pos_vec[i].to(cls.units_list[i + 1]).value for i in range(len(pos_vec))
-        ]
-        vel_vec_vals = [
-            vel_vec[i].to(cls.units_list[i + 5]).value for i in range(len(vel_vec))
-        ]
-        return cls(
-            np.array(pos_vec_vals), np.array(vel_vec_vals), time.to(u.s), M.to(u.kg)
+        time_vel = schwarzschild_utils.time_velocity(pos_vec, vel_vec, M)
+        self.initial_vec = np.hstack(
+            (self.time.value, pos_vec, time_vel.value, vel_vec)
         )
+        self.scr = schwarzschild_radius(M)
 
     @classmethod
     @u.quantity_input(time=u.s, M=u.kg)
-    def from_spherical(cls, pos_vec, vel_vec, time, M):
+    def from_spherical(cls, coords, M, time=0 * u.s):
         """
         Constructor
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Object having both initial positions and velocities of particle in Spherical Coordinates
+        M : ~astropy.units.quantity.Quantity
             Mass of the body
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
-        cls.input_coord_system = "Spherical"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        return cls._classmethod_handler(pos_vec, vel_vec, time, M)
+        return cls(coords, M, time)
 
     @classmethod
     @u.quantity_input(time=u.s, M=u.kg)
-    def from_cartesian(cls, pos_vec, vel_vec, time, M):
+    def from_cartesian(cls, coords, M, time=0 * u.s):
         """
         Constructor
 
         Parameters
         ----------
-        pos_vec : list
-            list of x, y and z components along with ~astropy.units
-        vel_vec : list
-            list of velocities of x, y, and z components along with ~astropy.units
-        time : ~astropy.units.s
-            Time of start
-        M : ~astropy.units.kg
+        coords : ~einsteinpy.coordinates.velocity.CartesianDifferential
+            Object having both initial positions and velocities of particle in Cartesian Coordinates
+        M : ~astropy.units.quantity.Quantity
             Mass of the body
+        time : ~astropy.units.quantity.Quantity
+            Time of start, defaults to 0 seconds.
 
         """
-        cls.input_coord_system = "Cartesian"
-        cls.input_units_list = (
-            [time.unit]
-            + [pos_vec[i].unit for i in range(len(pos_vec))]
-            + [u.one]
-            + [vel_vec[i].unit for i in range(len(vel_vec))]
-        )
-        sp_pos_vec, sp_vel_vec = C2S_units(pos_vec, vel_vec)
-        return cls._classmethod_handler(sp_pos_vec, sp_vel_vec, time, M)
+        sph_coords = coords.spherical_differential()
+        return cls(sph_coords, M, time)
 
     def f_vec(self, ld, vec):
         vals = np.zeros(shape=vec.shape, dtype=vec.dtype)
@@ -135,9 +98,9 @@ class Schwarzschild:
         Parameters
         ----------
         start_lambda : float
-            Starting lambda, defaults to 0.0, (lambda ~= t)
+            Starting lambda(proper time), defaults to 0, (lambda ~= t)
         end_lamdba : float
-            Lambda where iteartions will stop, defaults to 100000
+            Lambda(proper time) where iteartions will stop, defaults to 100000
         stop_on_singularity : bool
             Whether to stop further computation on reaching singularity, defaults to True
         OdeMethodKwargs : dict
@@ -148,13 +111,12 @@ class Schwarzschild:
 
         Returns
         -------
-        tuple
-            (~numpy.array of lambda, (n,8) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        ~numpy.ndarray
+            N-element array containing proper time.
+        ~numpy.ndarray
+            (n,8) shape array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        vec_list = list()
-        lambda_list = list()
-        singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
             t0=start_lambda,
@@ -162,39 +124,44 @@ class Schwarzschild:
             t_bound=end_lambda,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
+
+        vecs = list()
+        lambdas = list()
+        singularity_reached = False
+        _scr = self.scr.value * 1.001
+
         while ODE.t < end_lambda:
-            vec_list.append(ODE.y)
-            lambda_list.append(ODE.t)
+            vecs.append(ODE.y)
+            lambdas.append(ODE.t)
             ODE.step()
             if (not singularity_reached) and (ODE.y[1] <= _scr):
-                warnings.warn(
-                    "r component of position vector reached Schwarzchild Radius. ",
-                    RuntimeWarning,
-                )
+                warnings.warn("particle reached Schwarzchild Radius. ", RuntimeWarning)
                 if stop_on_singularity:
                     break
                 else:
                     singularity_reached = True
 
-        def _not_cartesian():
-            return (np.array(lambda_list), np.array(vec_list))
-
-        def _cartesian():
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-            return (np.array(lambda_list), S2C_8dim(np.array(vec_list)))
-
-        choice_dict = {0: _not_cartesian, 1: _cartesian}
-        return choice_dict[int(return_cartesian)]()
+        vecs, lambdas = np.array(vecs), np.array(lambdas)
+
+        if not return_cartesian:
+            return lambdas, vecs
+        else:
+            cart_vecs = list()
+            for v in vecs:
+                si_vals = (
+                    SphericalDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
+                    )
+                    .cartesian_differential()
+                    .si_values()
+                )
+                cart_vecs.append(np.hstack((v[0], si_vals[:3], v[4], si_vals[3:])))
+            return lambdas, np.array(cart_vecs)
 
     def calculate_trajectory_iterator(
         self,
@@ -221,11 +188,12 @@ class Schwarzschild:
 
         Yields
         ------
-        tuple
-            (lambda, (8,) shape ~numpy.array of [t, pos1, pos2, pos3, velocity_of_time, vel1, vel2, vel3])
+        float
+            proper time
+        ~numpy.ndarray
+            array of [t, x1, x2, x3, velocity_of_time, v1, v2, v3] for each proper time(lambda).
 
         """
-        singularity_reached = False
         ODE = RK45(
             fun=self.f_vec,
             t0=start_lambda,
@@ -233,38 +201,32 @@ class Schwarzschild:
             t_bound=1e300,
             **OdeMethodKwargs
         )
-        _scr = self.schwarzschild_r.value * 1.001
 
-        def yielder_func():
-            nonlocal singularity_reached
-            while True:
-                if not return_cartesian:
-                    yield (ODE.t, ODE.y)
-                else:
-                    temp = np.copy(ODE.y)
-                    temp[1:4] = SphericalToCartesian_pos(ODE.y[1:4])
-                    temp[5:8] = SphericalToCartesian_vel(ODE.y[1:4], ODE.y[5:8])
-                    yield (ODE.t, temp)
-                ODE.step()
-                if (not singularity_reached) and (ODE.y[1] <= _scr):
-                    warnings.warn(
-                        "r component of position vector reached Schwarzchild Radius. ",
-                        RuntimeWarning,
+        singularity_reached = False
+        _scr = self.scr.value * 1.001
+
+        while True:
+            if not return_cartesian:
+                yield ODE.t, ODE.y
+            else:
+                v = ODE.y
+                si_vals = (
+                    SphericalDifferential(
+                        v[1] * u.m,
+                        v[2] * u.rad,
+                        v[3] * u.rad,
+                        v[5] * u.m / u.s,
+                        v[6] * u.rad / u.s,
+                        v[7] * u.rad / u.s,
                     )
-                    if stop_on_singularity:
-                        break
-                    else:
-                        singularity_reached = True
-
-        if return_cartesian:
-            self.units_list = [
-                u.s,
-                u.m,
-                u.m,
-                u.m,
-                u.one,
-                u.m / u.s,
-                u.m / u.s,
-                u.m / u.s,
-            ]
-        return yielder_func()
+                    .cartesian_differential()
+                    .si_values()
+                )
+                yield ODE.t, np.hstack((v[0], si_vals[:3], v[4], si_vals[3:]))
+            ODE.step()
+            if (not singularity_reached) and (ODE.y[1] <= _scr):
+                warnings.warn("particle reached Schwarzchild Radius. ", RuntimeWarning)
+                if stop_on_singularity:
+                    break
+                else:
+                    singularity_reached = True
diff --git a/src/einsteinpy/plotting/geodesics_scatter.py b/src/einsteinpy/plotting/geodesics_scatter.py
index d2a89d0..ac4a6c6 100644
--- a/src/einsteinpy/plotting/geodesics_scatter.py
+++ b/src/einsteinpy/plotting/geodesics_scatter.py
@@ -34,15 +34,13 @@ class ScatterGeodesicPlotter:
         self._attractor_present = True
         plt.scatter(0, 0, color=self.attractor_color)
 
-    def plot(self, pos_vec, vel_vec, end_lambda=10, step_size=1e-3):
+    def plot(self, coords, end_lambda=10, step_size=1e-3):
         """
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units.
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units.
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Position and velocity components of particle in Spherical Coordinates.
         end_lambda : float, optional
             Lambda where iteartions will stop.
         step_size : float, optional
@@ -50,7 +48,7 @@ class ScatterGeodesicPlotter:
 
         """
 
-        swc = Schwarzschild.from_spherical(pos_vec, vel_vec, self.time, self.mass)
+        swc = Schwarzschild.from_spherical(coords, self.mass, self.time)
 
         vals = swc.calculate_trajectory(
             end_lambda=end_lambda, OdeMethodKwargs={"stepsize": step_size}
@@ -70,16 +68,14 @@ class ScatterGeodesicPlotter:
         if not self._attractor_present:
             self._plot_attractor()
 
-    def animate(self, pos_vec, vel_vec, end_lambda=10, step_size=1e-3, interval=50):
+    def animate(self, coords, end_lambda=10, step_size=1e-3, interval=50):
         """
         Function to generate animated plots of geodesics.
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units.
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units.
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Position and velocity components of particle in Spherical Coordinates.
         end_lambda : float, optional
             Lambda where iteartions will stop.
         step_size : float, optional
@@ -89,7 +85,7 @@ class ScatterGeodesicPlotter:
 
         """
 
-        swc = Schwarzschild.from_spherical(pos_vec, vel_vec, self.time, self.mass)
+        swc = Schwarzschild.from_spherical(coords, self.mass, self.time)
 
         vals = swc.calculate_trajectory(
             end_lambda=end_lambda, OdeMethodKwargs={"stepsize": step_size}
diff --git a/src/einsteinpy/plotting/geodesics_static.py b/src/einsteinpy/plotting/geodesics_static.py
index a1d476d..476b91c 100644
--- a/src/einsteinpy/plotting/geodesics_static.py
+++ b/src/einsteinpy/plotting/geodesics_static.py
@@ -45,15 +45,13 @@ class StaticGeodesicPlotter:
         self.__yarr = np.array([])
         self.get_curr_plot_radius = 0
 
-    def plot_trajectory(self, pos_vec, vel_vec, end_lambda, step_size, color):
+    def plot_trajectory(self, coords, end_lambda, step_size, color):
         """
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units.
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units.
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Initial position and velocity of particle in Spherical coordinates.
         end_lambda : float, optional
             Lambda where iteartions will stop.
         step_size : float, optional
@@ -62,15 +60,15 @@ class StaticGeodesicPlotter:
             Color of the Geodesic
 
         """
-        swc = Schwarzschild.from_spherical(pos_vec, vel_vec, self.time, self.mass)
+        swc = Schwarzschild.from_spherical(coords, self.mass, self.time)
 
         vals = swc.calculate_trajectory(
             end_lambda=end_lambda, OdeMethodKwargs={"stepsize": step_size}
         )[1]
 
-        time = np.array([coord[0] for coord in vals])
+        # time = np.array([coord[0] for coord in vals])
         r = np.array([coord[1] for coord in vals])
-        theta = np.array([coord[2] for coord in vals])
+        # theta = np.array([coord[2] for coord in vals])
         phi = np.array([coord[3] for coord in vals])
 
         x = r * np.cos(phi)
@@ -84,30 +82,28 @@ class StaticGeodesicPlotter:
         if not self._attractor_present:
             self._draw_attractor()
 
-    def __get_x_y(self, pos_vec, vel_vec, end_lambda, step_size):
+    def __get_x_y(self, coords, end_lambda, step_size):
         """
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units.
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units.
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Initial position and velocity of particle in Spherical coordinates.
         end_lambda : float, optional
             Lambda where iteartions will stop.
         step_size : float, optional
             Step size for the ODE.
 
         """
-        swc = Schwarzschild.from_spherical(pos_vec, vel_vec, self.time, self.mass)
+        swc = Schwarzschild.from_spherical(coords, self.mass, self.time)
 
         vals = swc.calculate_trajectory(
             end_lambda=end_lambda, OdeMethodKwargs={"stepsize": step_size}
         )[1]
 
-        time = np.array([coord[0] for coord in vals])
+        # time = np.array([coord[0] for coord in vals])
         r = np.array([coord[1] for coord in vals])
-        theta = np.array([coord[2] for coord in vals])
+        # theta = np.array([coord[2] for coord in vals])
         phi = np.array([coord[3] for coord in vals])
 
         x = r * np.cos(phi)
@@ -150,8 +146,7 @@ class StaticGeodesicPlotter:
 
     def plot(
         self,
-        pos_vec,
-        vel_vec,
+        coords,
         end_lambda=10,
         step_size=1e-3,
         color="#{:06x}".format(random.randint(0, 0xFFFFFF)),
@@ -160,10 +155,8 @@ class StaticGeodesicPlotter:
 
         Parameters
         ----------
-        pos_vec : list
-            list of r, theta & phi components along with ~astropy.units.
-        vel_vec : list
-            list of velocities of r, theta & phi components along with ~astropy.units.
+        coords : ~einsteinpy.coordinates.velocity.SphericalDifferential
+            Initial position and velocity of particle in Spherical coordinates.
         end_lambda : float, optional
             Lambda where iteartions will stop.
         step_size : float, optional
@@ -172,15 +165,11 @@ class StaticGeodesicPlotter:
             Color of the Geodesic
 
         """
-        self.__xarr, self.__yarr = self.__get_x_y(
-            pos_vec, vel_vec, end_lambda, step_size
-        )
+        self.__xarr, self.__yarr = self.__get_x_y(coords, end_lambda, step_size)
         self.plot_attractor()
         self._attractor_present = True
 
-        lines, x0, y0 = self.plot_trajectory(
-            pos_vec, vel_vec, end_lambda, step_size, color
-        )
+        lines, x0, y0 = self.plot_trajectory(coords, end_lambda, step_size, color)
 
         l, = self.ax.plot(x0, y0, "o", mew=0, color=lines[0].get_color())
         lines.append(l)


Test patch:
diff --git a/src/einsteinpy/tests/test_metric/test_kerr.py b/src/einsteinpy/tests/test_metric/test_kerr.py
index 4c0b755..d24ec80 100644
--- a/src/einsteinpy/tests/test_metric/test_kerr.py
+++ b/src/einsteinpy/tests/test_metric/test_kerr.py
@@ -6,6 +6,7 @@ from astropy import units as u
 from numpy.testing import assert_allclose
 
 from einsteinpy import constant
+from einsteinpy.coordinates import BoyerLindquistDifferential, CartesianDifferential
 from einsteinpy.metric import Kerr
 from einsteinpy.utils import kerr_utils, schwarzschild_radius_dimensionless
 
@@ -13,49 +14,63 @@ _c = constant.c.value
 
 
 @pytest.mark.parametrize(
-    "pos_vec, vel_vec, time, M, a, start_lambda, end_lambda, OdeMethodKwargs",
+    "coords, time, M, start_lambda, end_lambda, OdeMethodKwargs",
     [
         (
-            [306 * u.m, np.pi / 2.05 * u.rad, np.pi / 2 * u.rad],
-            [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s],
+            BoyerLindquistDifferential(
+                306 * u.m,
+                np.pi / 2.05 * u.rad,
+                np.pi / 2 * u.rad,
+                0 * u.m / u.s,
+                0 * u.rad / u.s,
+                951.0 * u.rad / u.s,
+                2e-3 * u.m,
+            ),
             0 * u.s,
             4e24 * u.kg,
-            2e-3,
             0.0,
             0.001,
             {"stepsize": 0.5e-6},
         ),
         (
-            [1 * u.km, 0.15 * u.rad, np.pi / 2 * u.rad],
-            [
+            BoyerLindquistDifferential(
+                1 * u.km,
+                0.15 * u.rad,
+                np.pi / 2 * u.rad,
                 0.1 * _c * u.m / u.s,
                 0.5e-5 * _c * u.rad / u.s,
                 0.5e-4 * _c * u.rad / u.s,
-            ],
+                2e-3 * u.m,
+            ),
             0 * u.s,
             5.972e24 * u.kg,
-            2e-3,
             0.0,
             0.0001,
             {"stepsize": 0.5e-6},
         ),
         (
-            [50 * u.km, np.pi / 2 * u.rad, np.pi / 2 * u.rad],
-            [0.1 * _c * u.m / u.s, 2e-7 * _c * u.rad / u.s, 1e-5 * u.rad / u.s],
+            BoyerLindquistDifferential(
+                50 * u.km,
+                np.pi / 2 * u.rad,
+                np.pi / 2 * u.rad,
+                0.1 * _c * u.m / u.s,
+                2e-7 * _c * u.rad / u.s,
+                1e-5 * u.rad / u.s,
+                0.0 * u.km,
+            ),
             0 * u.s,
             5.972e24 * u.g,
             0.0,
-            0.0,
             0.001,
             {"stepsize": 5e-6},
         ),
     ],
 )
 def test_calculate_trajectory(
-    pos_vec, vel_vec, time, M, a, start_lambda, end_lambda, OdeMethodKwargs
+    coords, time, M, start_lambda, end_lambda, OdeMethodKwargs
 ):
     _scr = schwarzschild_radius_dimensionless(M)
-    obj = Kerr.from_BL(pos_vec, vel_vec, time, M, a)
+    obj = Kerr.from_BL(coords, M, time)
     ans = obj.calculate_trajectory(
         start_lambda=start_lambda,
         end_lambda=end_lambda,
@@ -64,7 +79,7 @@ def test_calculate_trajectory(
     ans = ans[1]
     testarray = list()
     for i in ans:
-        g = kerr_utils.metric(i[1], i[2], M.value, a)
+        g = kerr_utils.metric(i[1], i[2], M.value, coords.a.to(u.m).value)
         testarray.append(
             g[0][0] * (i[4] ** 2)
             + g[1][1] * (i[5] ** 2)
@@ -85,18 +100,17 @@ def test_calculate_trajectory3():
     M = 1.989e30 * u.kg
     distance_at_perihelion = 147.10e6 * u.km
     speed_at_perihelion = 30.29 * u.km / u.s
-    pos_vec = [
+    cart_obj = CartesianDifferential(
         distance_at_perihelion / np.sqrt(2),
         distance_at_perihelion / np.sqrt(2),
         0 * u.km,
-    ]
-    vel_vec = [
         -1 * speed_at_perihelion / np.sqrt(2),
         speed_at_perihelion / np.sqrt(2),
         0 * u.km / u.h,
-    ]
+    )
     end_lambda = ((1 * u.year).to(u.s)).value
-    cl = Kerr.from_cartesian(pos_vec, vel_vec, 0 * u.min, M, 0.0)
+    a = 0 * u.km
+    cl = Kerr.from_cartesian(cart_obj, M, a)
     ans = cl.calculate_trajectory(
         start_lambda=0.0,
         end_lambda=end_lambda,
@@ -115,26 +129,38 @@ def test_calculate_trajectory3():
 
 
 @pytest.mark.parametrize(
-    "pos_vec, vel_vec, time, M, a, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian",
+    "coords, time, M, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian",
     [
         (
-            [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad],
-            [0 * u.m / u.s, 0.1 * u.rad / u.s, 951.0 * u.rad / u.s],
+            BoyerLindquistDifferential(
+                306 * u.m,
+                np.pi / 2 * u.rad,
+                np.pi / 2 * u.rad,
+                0 * u.m / u.s,
+                0.1 * u.rad / u.s,
+                951.0 * u.rad / u.s,
+                2e-3 * u.m,
+            ),
             0 * u.s,
             4e24 * u.kg,
-            2e-3,
             0.0,
             0.0003,
             {"stepsize": 0.3e-6},
             True,
         ),
         (
-            [1 * u.km, 0.15 * u.rad, np.pi / 2 * u.rad],
-            [_c * u.m / u.s, 0.5e-5 * _c * u.rad / u.s, 1e-4 * _c * u.rad / u.s],
+            BoyerLindquistDifferential(
+                1 * u.km,
+                0.15 * u.rad,
+                np.pi / 2 * u.rad,
+                0.2 * _c * u.m / u.s,
+                0.5e-5 * _c * u.rad / u.s,
+                1e-4 * _c * u.rad / u.s,
+                0.0 * u.km,
+            ),
             0 * u.s,
             5.972e24 * u.kg,
             0.0,
-            0.0,
             0.0004,
             {"stepsize": 0.5e-6},
             False,
@@ -142,24 +168,16 @@ def test_calculate_trajectory3():
     ],
 )
 def test_calculate_trajectory_iterator(
-    pos_vec,
-    vel_vec,
-    time,
-    M,
-    a,
-    start_lambda,
-    end_lambda,
-    OdeMethodKwargs,
-    return_cartesian,
+    coords, time, M, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian
 ):
-    cl1 = Kerr.from_BL(pos_vec, vel_vec, time, M, a)
+    cl1 = Kerr.from_BL(coords, M, time)
     arr1 = cl1.calculate_trajectory(
         start_lambda=start_lambda,
         end_lambda=end_lambda,
         OdeMethodKwargs=OdeMethodKwargs,
         return_cartesian=return_cartesian,
     )[1]
-    cl2 = Kerr.from_BL(pos_vec, vel_vec, time, M, a)
+    cl2 = Kerr.from_BL(coords, M, time)
     it = cl2.calculate_trajectory_iterator(
         start_lambda=start_lambda,
         OdeMethodKwargs=OdeMethodKwargs,
@@ -170,3 +188,28 @@ def test_calculate_trajectory_iterator(
         arr2_list.append(val[1])
     arr2 = np.array(arr2_list)
     assert_allclose(arr1[:100, :], arr2, rtol=1e-10)
+
+
+def test_calculate_trajectory_iterator_RuntimeWarning():
+    bl_obj = BoyerLindquistDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 2 * u.rad,
+        0 * u.m / u.s,
+        0.01 * u.rad / u.s,
+        10 * u.rad / u.s,
+        0 * u.m,
+    )
+    M = 1e25 * u.kg
+    start_lambda = 0.0
+    OdeMethodKwargs = {"stepsize": 0.4e-6}
+    cl = Kerr.from_BL(bl_obj, M)
+    with warnings.catch_warnings(record=True) as w:
+        it = cl.calculate_trajectory_iterator(
+            start_lambda=start_lambda,
+            OdeMethodKwargs=OdeMethodKwargs,
+            stop_on_singularity=True,
+        )
+        for _, _ in zip(range(1000), it):
+            pass
+        assert len(w) >= 1
diff --git a/src/einsteinpy/tests/test_metric/test_kerrnewman.py b/src/einsteinpy/tests/test_metric/test_kerrnewman.py
index d89cd7f..858a68d 100644
--- a/src/einsteinpy/tests/test_metric/test_kerrnewman.py
+++ b/src/einsteinpy/tests/test_metric/test_kerrnewman.py
@@ -6,6 +6,7 @@ from astropy import units as u
 from numpy.testing import assert_allclose
 
 from einsteinpy import constant
+from einsteinpy.coordinates import BoyerLindquistDifferential, CartesianDifferential
 from einsteinpy.metric import KerrNewman
 from einsteinpy.utils import kerrnewman_utils, schwarzschild_radius
 
@@ -21,20 +22,20 @@ def test_calculate_trajectory0():
     # Function returning cartesian coordinates
     M = 1.989e30 * u.kg
     q = 0 * u.C / u.kg
+    Q = 0 * u.C
     distance_at_perihelion = 147.10e6 * u.km
     speed_at_perihelion = 30.29 * u.km / u.s
-    pos_vec = [
+    cart_obj = CartesianDifferential(
         distance_at_perihelion / np.sqrt(2),
         distance_at_perihelion / np.sqrt(2),
         0 * u.km,
-    ]
-    vel_vec = [
         -1 * speed_at_perihelion / np.sqrt(2),
         speed_at_perihelion / np.sqrt(2),
         0 * u.km / u.h,
-    ]
+    )
+    a = 0 * u.m
     end_lambda = ((1 * u.year).to(u.s)).value
-    cl = KerrNewman.from_cartesian(pos_vec, vel_vec, q, 0 * u.min, M, 0.0, 0 * u.C)
+    cl = KerrNewman.from_cartesian(cart_obj, q, M, a, Q)
     ans = cl.calculate_trajectory(
         start_lambda=0.0,
         end_lambda=end_lambda,
@@ -60,11 +61,18 @@ def test_calculate_trajectory1():
     stepsize = 0.5
     tmp = _G * M.value / _cc
     q = tmp * u.C / u.kg
-    print(tmp, 5.900455 * tmp ** 3)
+    a = 0 * u.km
     Q = 11604461683.91822052001953125 * u.C
-    pos_vec = [r, 0.5 * np.pi * u.rad, 0 * u.rad]
-    vel_vec = [0 * u.m / u.s, 0 * u.rad / u.s, 0.0 * u.rad / u.s]
-    cl = KerrNewman.from_BL(pos_vec, vel_vec, q, 0 * u.s, M, 0.0, Q)
+    bl_obj = BoyerLindquistDifferential(
+        r,
+        0.5 * np.pi * u.rad,
+        0 * u.rad,
+        0 * u.m / u.s,
+        0 * u.rad / u.s,
+        0.0 * u.rad / u.s,
+        a,
+    )
+    cl = KerrNewman.from_BL(bl_obj, q, M, Q)
     ans = cl.calculate_trajectory(
         end_lambda=end_lambda, OdeMethodKwargs={"stepsize": stepsize}
     )
@@ -74,7 +82,7 @@ def test_calculate_trajectory1():
 @pytest.fixture()
 def test_input():
     q = 1 * u.C / u.g
-    a = 1e-6
+    a = 1e-6 * u.m
     Q = 100 * u.C
     el = 200.0
     ss = 1.0
@@ -82,12 +90,19 @@ def test_input():
 
 
 def test_compare_calculate_trajectory_iterator_bl(test_input):
-    pos_vec = [1000.0 * u.km, 0.6 * np.pi * u.rad, np.pi / 8 * u.rad]
-    vel_vec = [10000 * u.m / u.s, -0.01 * u.rad / u.s, 0.0 * u.rad / u.s]
-    M = 0.5 * 5.972e24 * u.kg
     q, a, Q, el, ss = test_input
-    cl1 = KerrNewman.from_BL(pos_vec, vel_vec, q, 0 * u.s, M, a, Q)
-    cl2 = KerrNewman.from_BL(pos_vec, vel_vec, q, 0 * u.s, M, a, Q)
+    bl_obj = BoyerLindquistDifferential(
+        1000.0 * u.km,
+        0.6 * np.pi * u.rad,
+        np.pi / 8 * u.rad,
+        10000 * u.m / u.s,
+        -0.01 * u.rad / u.s,
+        0.0 * u.rad / u.s,
+        a,
+    )
+    M = 0.5 * 5.972e24 * u.kg
+    cl1 = KerrNewman.from_BL(bl_obj, q, M, Q)
+    cl2 = KerrNewman.from_BL(bl_obj, q, M, Q)
     ans1 = cl1.calculate_trajectory(end_lambda=el, OdeMethodKwargs={"stepsize": ss})[1]
     it = cl2.calculate_trajectory_iterator(OdeMethodKwargs={"stepsize": ss})
     ans2 = list()
@@ -99,12 +114,18 @@ def test_compare_calculate_trajectory_iterator_bl(test_input):
 
 
 def test_compare_calculate_trajectory_iterator_cartesians(test_input):
-    pos_vec = [1000000 * u.m, 1000000 * u.m, 20.5 * u.m]
-    vel_vec = [10000 * u.m / u.s, 10000 * u.m / u.s, -30 * u.m / u.s]
+    cart_obj = CartesianDifferential(
+        1000000 * u.m,
+        1000000 * u.m,
+        20.5 * u.m,
+        10000 * u.m / u.s,
+        10000 * u.m / u.s,
+        -30 * u.m / u.s,
+    )
     M = 2e24 * u.kg
     q, a, Q, el, ss = test_input
-    cl1 = KerrNewman.from_cartesian(pos_vec, vel_vec, q, 0 * u.s, M, a, Q)
-    cl2 = KerrNewman.from_cartesian(pos_vec, vel_vec, q, 0 * u.s, M, a, Q)
+    cl1 = KerrNewman.from_cartesian(cart_obj, q, M, a, Q)
+    cl2 = KerrNewman.from_cartesian(cart_obj, q, M, a, Q)
     ans1 = cl1.calculate_trajectory(
         end_lambda=el, OdeMethodKwargs={"stepsize": ss}, return_cartesian=True
     )[1]
@@ -117,3 +138,29 @@ def test_compare_calculate_trajectory_iterator_cartesians(test_input):
     ans2 = np.array(ans2)
     print(ans1)
     assert_allclose(ans1[:20], ans2)
+
+
+def test_calculate_trajectory_iterator_RuntimeWarning():
+    bl_obj = BoyerLindquistDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 2 * u.rad,
+        0 * u.m / u.s,
+        0.01 * u.rad / u.s,
+        10 * u.rad / u.s,
+        0 * u.m,
+    )
+    M = 1e25 * u.kg
+    start_lambda = 0.0
+    q, Q = 0 * u.C / u.kg, 0 * u.C
+    OdeMethodKwargs = {"stepsize": 0.4e-6}
+    cl = KerrNewman.from_BL(bl_obj, q, M, Q)
+    with warnings.catch_warnings(record=True) as w:
+        it = cl.calculate_trajectory_iterator(
+            start_lambda=start_lambda,
+            OdeMethodKwargs=OdeMethodKwargs,
+            stop_on_singularity=True,
+        )
+        for _, _ in zip(range(1000), it):
+            pass
+        assert len(w) >= 1
diff --git a/src/einsteinpy/tests/test_metric/test_schwarzschild.py b/src/einsteinpy/tests/test_metric/test_schwarzschild.py
index 06acdcf..a15d706 100644
--- a/src/einsteinpy/tests/test_metric/test_schwarzschild.py
+++ b/src/einsteinpy/tests/test_metric/test_schwarzschild.py
@@ -6,6 +6,7 @@ from astropy import units as u
 from numpy.testing import assert_allclose
 
 from einsteinpy import constant
+from einsteinpy.coordinates import CartesianDifferential, SphericalDifferential
 from einsteinpy.metric import Schwarzschild
 from einsteinpy.utils import schwarzschild_radius, schwarzschild_utils
 
@@ -13,11 +14,17 @@ _c = constant.c.value
 
 
 @pytest.mark.parametrize(
-    "pos_vec, vel_vec, time, M, start_lambda, end_lambda, OdeMethodKwargs",
+    "coords, time, M, start_lambda, end_lambda, OdeMethodKwargs",
     [
         (
-            [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad],
-            [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s],
+            SphericalDifferential(
+                306 * u.m,
+                np.pi / 2 * u.rad,
+                np.pi / 2 * u.rad,
+                0 * u.m / u.s,
+                0 * u.rad / u.s,
+                951.0 * u.rad / u.s,
+            ),
             0 * u.s,
             4e24 * u.kg,
             0.0,
@@ -25,12 +32,14 @@ _c = constant.c.value
             {"stepsize": 0.5e-6},
         ),
         (
-            [1 * u.km, 0.15 * u.rad, np.pi / 2 * u.rad],
-            [
+            SphericalDifferential(
+                1 * u.km,
+                0.15 * u.rad,
+                90 * u.deg,
                 0.1 * _c * u.m / u.s,
                 0.5e-5 * _c * u.rad / u.s,
                 0.5e-4 * _c * u.rad / u.s,
-            ],
+            ),
             0 * u.s,
             5.972e24 * u.kg,
             0.0,
@@ -38,8 +47,14 @@ _c = constant.c.value
             {"stepsize": 0.5e-6},
         ),
         (
-            [50 * u.km, np.pi / 2 * u.rad, np.pi / 2 * u.rad],
-            [0.1 * _c * u.m / u.s, 2e-7 * _c * u.rad / u.s, 1e-5 * u.rad / u.s],
+            SphericalDifferential(
+                50 * u.km,
+                np.pi / 2 * u.rad,
+                np.pi / 2 * u.rad,
+                0.1 * _c * u.m / u.s,
+                2e-7 * _c * u.rad / u.s,
+                1e-5 * u.rad / u.s,
+            ),
             0 * u.s,
             5.972e24 * u.g,
             0.0,
@@ -49,9 +64,9 @@ _c = constant.c.value
     ],
 )
 def test_calculate_trajectory(
-    pos_vec, vel_vec, time, M, start_lambda, end_lambda, OdeMethodKwargs
+    coords, time, M, start_lambda, end_lambda, OdeMethodKwargs
 ):
-    cl = Schwarzschild.from_spherical(pos_vec, vel_vec, time, M)
+    cl = Schwarzschild.from_spherical(coords, M, time)
     ans = cl.calculate_trajectory(
         start_lambda=start_lambda,
         end_lambda=end_lambda,
@@ -79,10 +94,16 @@ def test_calculate_trajectory2():
     distance_at_perihelion = 147.10e6 * u.km
     speed_at_perihelion = 30.29 * u.km / u.s
     angular_vel = (speed_at_perihelion / distance_at_perihelion) * u.rad
-    pos_vec = [distance_at_perihelion, np.pi / 2 * u.rad, 0 * u.rad]
-    vel_vec = [0 * u.km / u.s, 0 * u.rad / u.s, angular_vel]
+    sph_obj = SphericalDifferential(
+        distance_at_perihelion,
+        np.pi / 2 * u.rad,
+        0 * u.rad,
+        0 * u.km / u.s,
+        0 * u.rad / u.s,
+        angular_vel,
+    )
     end_lambda = ((1 * u.year).to(u.s)).value
-    cl = Schwarzschild.from_spherical(pos_vec, vel_vec, 0 * u.s, M)
+    cl = Schwarzschild.from_spherical(sph_obj, M)
     ans = cl.calculate_trajectory(
         start_lambda=0.0,
         end_lambda=end_lambda,
@@ -101,18 +122,16 @@ def test_calculate_trajectory3():
     M = 1.989e30 * u.kg
     distance_at_perihelion = 147.10e6 * u.km
     speed_at_perihelion = 30.29 * u.km / u.s
-    pos_vec = [
+    cart_obj = CartesianDifferential(
         distance_at_perihelion / np.sqrt(2),
         distance_at_perihelion / np.sqrt(2),
         0 * u.km,
-    ]
-    vel_vec = [
         -1 * speed_at_perihelion / np.sqrt(2),
         speed_at_perihelion / np.sqrt(2),
         0 * u.km / u.h,
-    ]
+    )
     end_lambda = ((1 * u.year).to(u.s)).value
-    cl = Schwarzschild.from_cartesian(pos_vec, vel_vec, 0 * u.min, M)
+    cl = Schwarzschild.from_cartesian(cart_obj, M)
     ans = cl.calculate_trajectory(
         start_lambda=0.0,
         end_lambda=end_lambda,
@@ -131,68 +150,78 @@ def test_calculate_trajectory3():
 
 
 @pytest.mark.parametrize(
-    "pos_vec, vel_vec, time, M, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian",
+    "coords, time, M, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian",
     [
         (
-            [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad],
-            [0 * u.m / u.s, 0.1 * u.rad / u.s, 951.0 * u.rad / u.s],
+            SphericalDifferential(
+                306 * u.m,
+                np.pi / 2 * u.rad,
+                np.pi / 2 * u.rad,
+                0 * u.m / u.s,
+                0.1 * u.rad / u.s,
+                951.0 * u.rad / u.s,
+            ),
             0 * u.s,
             4e24 * u.kg,
             0.0,
-            0.0003,
+            0.0002,
             {"stepsize": 0.3e-6},
             True,
         ),
         (
-            [1 * u.km, 0.15 * u.rad, np.pi / 2 * u.rad],
-            [_c * u.m / u.s, 0.5e-5 * _c * u.rad / u.s, 1e-4 * _c * u.rad / u.s],
+            SphericalDifferential(
+                1 * u.km,
+                0.15 * u.rad,
+                np.pi / 2 * u.rad,
+                _c * u.m / u.s,
+                0.5e-5 * _c * u.rad / u.s,
+                1e-4 * _c * u.rad / u.s,
+            ),
             0 * u.s,
             5.972e24 * u.kg,
             0.0,
-            0.0004,
+            0.0002,
             {"stepsize": 0.5e-6},
             False,
         ),
     ],
 )
 def test_calculate_trajectory_iterator(
-    pos_vec,
-    vel_vec,
-    time,
-    M,
-    start_lambda,
-    end_lambda,
-    OdeMethodKwargs,
-    return_cartesian,
+    coords, time, M, start_lambda, end_lambda, OdeMethodKwargs, return_cartesian
 ):
-    cl1 = Schwarzschild.from_spherical(pos_vec, vel_vec, time, M)
+    cl1 = Schwarzschild.from_spherical(coords, M, time)
     arr1 = cl1.calculate_trajectory(
         start_lambda=start_lambda,
         end_lambda=end_lambda,
         OdeMethodKwargs=OdeMethodKwargs,
         return_cartesian=return_cartesian,
     )[1]
-    cl2 = Schwarzschild.from_spherical(pos_vec, vel_vec, time, M)
+    cl2 = Schwarzschild.from_spherical(coords, M, time)
     it = cl2.calculate_trajectory_iterator(
         start_lambda=start_lambda,
         OdeMethodKwargs=OdeMethodKwargs,
         return_cartesian=return_cartesian,
     )
     arr2_list = list()
-    for _, val in zip(range(100), it):
+    for _, val in zip(range(50), it):
         arr2_list.append(val[1])
     arr2 = np.array(arr2_list)
-    assert_allclose(arr1[:100, :], arr2, rtol=1e-10)
+    assert_allclose(arr1[:50, :], arr2, rtol=1e-10)
 
 
 def test_calculate_trajectory_iterator_RuntimeWarning():
-    pos_vec = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    vel_vec = [0 * u.m / u.s, 0.01 * u.rad / u.s, 10 * u.rad / u.s]
-    time = 0 * u.s
+    sph_obj = SphericalDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 2 * u.rad,
+        0 * u.m / u.s,
+        0.01 * u.rad / u.s,
+        10 * u.rad / u.s,
+    )
     M = 1e25 * u.kg
     start_lambda = 0.0
     OdeMethodKwargs = {"stepsize": 0.4e-6}
-    cl = Schwarzschild.from_spherical(pos_vec, vel_vec, time, M)
+    cl = Schwarzschild.from_spherical(sph_obj, M)
     with warnings.catch_warnings(record=True) as w:
         it = cl.calculate_trajectory_iterator(
             start_lambda=start_lambda,
@@ -205,13 +234,18 @@ def test_calculate_trajectory_iterator_RuntimeWarning():
 
 
 def test_calculate_trajectory_iterator_RuntimeWarning2():
-    pos_vec = [306 * u.m, np.pi / 2 * u.rad, np.pi / 3 * u.rad]
-    vel_vec = [0 * u.m / u.s, 0.01 * u.rad / u.s, 10 * u.rad / u.s]
-    time = 0 * u.s
+    sph_obj = SphericalDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 3 * u.rad,
+        0 * u.m / u.s,
+        0.01 * u.rad / u.s,
+        10 * u.rad / u.s,
+    )
     M = 1e25 * u.kg
     start_lambda = 0.0
     OdeMethodKwargs = {"stepsize": 0.4e-6}
-    cl = Schwarzschild.from_spherical(pos_vec, vel_vec, time, M)
+    cl = Schwarzschild.from_spherical(sph_obj, M)
     with warnings.catch_warnings(record=True) as w:
         it = cl.calculate_trajectory_iterator(
             start_lambda=start_lambda,
diff --git a/src/einsteinpy/tests/test_plotting/test_scattergeodesicplotter.py b/src/einsteinpy/tests/test_plotting/test_scattergeodesicplotter.py
index 053c57b..049df62 100644
--- a/src/einsteinpy/tests/test_plotting/test_scattergeodesicplotter.py
+++ b/src/einsteinpy/tests/test_plotting/test_scattergeodesicplotter.py
@@ -4,77 +4,68 @@ import astropy.units as u
 import numpy as np
 import pytest
 
+from einsteinpy.coordinates import SphericalDifferential
 from einsteinpy.plotting import ScatterGeodesicPlotter
 
 
 @pytest.fixture()
 def dummy_data():
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
+    sph_obj = SphericalDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 2 * u.rad,
+        0 * u.m / u.s,
+        0 * u.rad / u.s,
+        951.0 * u.rad / u.s,
+    )
     t = 0 * u.s
     m = 4e24 * u.kg
     start_lambda = 0.0
     end_lambda = 0.002
     step_size = 0.5e-6
-    return r, v, t, m, start_lambda, end_lambda, step_size
+    return sph_obj, t, m, start_lambda, end_lambda, step_size
 
 
 def test_plot_attractor_is_called_only_once(dummy_data):
-    r, v, _, m, _, el, ss = dummy_data
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = ScatterGeodesicPlotter(m)
     assert not cl._attractor_present
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     assert cl._attractor_present
 
 
 @mock.patch(
     "einsteinpy.plotting.geodesics_scatter.ScatterGeodesicPlotter._plot_attractor"
 )
-def test_plot_calls_plot_attractor(mock_plot_attractor):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_plot_calls_plot_attractor(mock_plot_attractor, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = ScatterGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     mock_plot_attractor.assert_called_with()
 
 
 @mock.patch("einsteinpy.plotting.geodesics_scatter.plt.figure")
-def test_animate_calls_figure(mock_figure):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_animate_calls_figure(mock_figure, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = ScatterGeodesicPlotter(m)
-    cl.animate(r, v, el, ss)
+    cl.animate(sph_obj, el, ss)
     mock_figure.assert_called_with()
 
 
 @mock.patch("einsteinpy.plotting.geodesics_scatter.plt.show")
-def test_plot_show_shows_plot(mock_show):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_plot_show_shows_plot(mock_show, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = ScatterGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     cl.show()
     mock_show.assert_called_with()
 
 
 @mock.patch("einsteinpy.plotting.geodesics_scatter.plt.savefig")
-def test_plot_save_saves_plot(mock_save):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_plot_save_saves_plot(mock_save, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = ScatterGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     name = "test_plot.png"
     cl.save(name)
     mock_save.assert_called_with(name)
diff --git a/src/einsteinpy/tests/test_plotting/test_staticgeodesicplotter.py b/src/einsteinpy/tests/test_plotting/test_staticgeodesicplotter.py
index 88f2fc0..3a11c77 100644
--- a/src/einsteinpy/tests/test_plotting/test_staticgeodesicplotter.py
+++ b/src/einsteinpy/tests/test_plotting/test_staticgeodesicplotter.py
@@ -6,23 +6,30 @@ import matplotlib.pyplot as plt
 import numpy as np
 import pytest
 
+from einsteinpy.coordinates import SphericalDifferential
 from einsteinpy.plotting import StaticGeodesicPlotter
 
 
 @pytest.fixture()
 def dummy_data():
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
+    sph_obj = SphericalDifferential(
+        306 * u.m,
+        np.pi / 2 * u.rad,
+        np.pi / 2 * u.rad,
+        0 * u.m / u.s,
+        0 * u.rad / u.s,
+        951.0 * u.rad / u.s,
+    )
     t = 0 * u.s
     m = 4e24 * u.kg
     start_lambda = 0.0
     end_lambda = 0.002
     step_size = 0.5e-6
-    return r, v, t, m, start_lambda, end_lambda, step_size
+    return sph_obj, t, m, start_lambda, end_lambda, step_size
 
 
 def test_staticgeodesicplotter_has_axes(dummy_data):
-    r, v, _, m, _, el, ss = dummy_data
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = StaticGeodesicPlotter(m)
     assert isinstance(cl.ax, mpl.axes.SubplotBase)
     assert cl.time.value == 0.0
@@ -30,45 +37,37 @@ def test_staticgeodesicplotter_has_axes(dummy_data):
 
 
 @mock.patch("einsteinpy.plotting.geodesics_static.plt.show")
-def test_plot_calls_plt_show(mock_show):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_plot_calls_plt_show(mock_show, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = StaticGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     cl.show()
     mock_show.assert_called_with()
     assert cl._attractor_present
 
 
 @mock.patch("einsteinpy.plotting.geodesics_static.plt.savefig")
-def test_plot_save_saves_plot(mock_save):
-    r = [306 * u.m, np.pi / 2 * u.rad, np.pi / 2 * u.rad]
-    v = [0 * u.m / u.s, 0 * u.rad / u.s, 951.0 * u.rad / u.s]
-    m = 4e24 * u.kg
-    el = 0.002
-    ss = 0.5e-6
+def test_plot_save_saves_plot(mock_save, dummy_data):
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = StaticGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     name = "test_plot.png"
     cl.save(name)
     mock_save.assert_called_with(name)
 
 
 def test_plot_calls_draw_attractor_Manualscale(dummy_data):
-    r, v, _, m, _, el, ss = dummy_data
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = StaticGeodesicPlotter(m, attractor_radius_scale=1500)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     assert cl._attractor_present
     assert cl.attractor_radius_scale == 1500
     assert cl.get_curr_plot_radius != -1
 
 
 def test_plot_calls_draw_attractor_AutoScale(dummy_data):
-    r, v, _, m, _, el, ss = dummy_data
+    sph_obj, _, m, _, el, ss = dummy_data
     cl = StaticGeodesicPlotter(m)
-    cl.plot(r, v, el, ss)
+    cl.plot(sph_obj, el, ss)
     assert cl._attractor_present
     assert cl.get_curr_plot_radius != -1