diff --git a/examples/gravity/satellite_velocity.py b/examples/gravity/satellite_velocity.py
index 65ce97787..aa2b974d0 100644
--- a/examples/gravity/satellite_velocity.py
+++ b/examples/gravity/satellite_velocity.py
@@ -17,14 +17,14 @@
 
 solution = Eq(centripetal_acceleration_law.law.rhs, free_fall_law.law.rhs)
 solution_applied = solution.subs(centripetal_acceleration_law.radius_of_curvature,
-    free_fall_law.planet_radius + free_fall_law.height_above_surface)
+    free_fall_law.planet_radius + free_fall_law.elevation)
 
 # first solution is negative and corresponds to the backwards direction of velocity - ignore it
 satellite_velocity = solve(solution_applied, centripetal_acceleration_law.speed,
     dict=True)[1][centripetal_acceleration_law.speed]
 
 print(
-    f"The formula for satellite linear velocity is: {print_expression(simplify(satellite_velocity))}"
+    f"The formula for satellite linear velocity is:\n{print_expression(simplify(satellite_velocity))}"
 )
 
 ## As a curve radius we are having radius of the planet plus desired height of the orbit. Let's take Earth as an example and 100km height.
@@ -35,7 +35,7 @@
 required_velocity_expression = satellite_velocity.subs({
     free_fall_law.planet_mass: planet_mass,
     free_fall_law.planet_radius: planet_radius_,
-    free_fall_law.height_above_surface: height_above_surface_
+    free_fall_law.elevation: height_above_surface_
 })
 
 result_velocity = Quantity(required_velocity_expression)
diff --git a/plots/gravity/free_fall_acceleration.py b/plots/gravity/free_fall_acceleration.py
index 1e3a84979..fd40277a0 100644
--- a/plots/gravity/free_fall_acceleration.py
+++ b/plots/gravity/free_fall_acceleration.py
@@ -1,7 +1,7 @@
 from sympy import solve, symbols
 from sympy.plotting import plot
 from sympy.plotting.plot import MatplotlibBackend
-from symplyphysics import (print_expression, units, convert_to, Quantity)
+from symplyphysics import (print_expression, units, convert_to, Quantity, quantities)
 from symplyphysics.laws.gravity import free_fall_acceleration_from_height as acceleration
 
 print(f"Formula is:\n{print_expression(acceleration.law)}")
@@ -17,8 +17,8 @@
 result_acceleration = solved.subs({
     acceleration.planet_mass: EARTH_MASS,
     acceleration.planet_radius: EARTH_RADIUS,
-    acceleration.height_above_surface: height_above_ground,
-    acceleration.units.gravitational_constant: gravity_constant_value
+    acceleration.elevation: height_above_ground,
+    quantities.gravitational_constant: gravity_constant_value
 })
 
 print(
diff --git a/symplyphysics/laws/gravity/area_rate_of_change_is_proportional_to_angular_momentum.py b/symplyphysics/laws/gravity/area_rate_of_change_is_proportional_to_angular_momentum.py
index 7dddd8ea6..c927834d5 100644
--- a/symplyphysics/laws/gravity/area_rate_of_change_is_proportional_to_angular_momentum.py
+++ b/symplyphysics/laws/gravity/area_rate_of_change_is_proportional_to_angular_momentum.py
@@ -1,41 +1,54 @@
+"""
+Area rate of change is proportional to angular momentum
+=======================================================
+
+The **law of areas**, also known as **Kepler's second law of planetary motion**,
+states that a line the connects a planet to the attracting body (the Sun) sweeps
+out equal areas in the plane of the planet's orbit in equal time intervals, and
+its rate of change is proportional to the planet's angular momentum. It is
+equivalent to saying that the planet's angular momentum is conserved.
+
+#. Sivukhin D.V. (1979), *Obshchiy kurs fiziki* [General course of Physics], vol. 1, pp. 312—314.
+"""
+
 from sympy import Eq, Derivative
 from symplyphysics import (
     units,
     Quantity,
-    Symbol,
-    Function,
-    print_expression,
     validate_input,
     validate_output,
     symbols,
+    clone_as_function,
 )
 
-# Description
-## The law of areas, also known as Kepler's second law of planetary motion, states that
-## a line the connects a planet to the attracting body (the Sun) sweeps out equal areas
-## in the plane of the planet's orbit in equal time intervals, and its rate of change
-## is proportional to the planets' angular momentum.
+time = symbols.time
+"""
+:symbols:`time`.
+"""
 
-# Note: this law is equivalent to saying that the planet's angular momentum is conserved.
+area_swept = clone_as_function(symbols.area, [time])
+"""
+:symbols:`area` swept by the planet.
+"""
 
-# Law: dA/dt = L/(2*m)
-## A - sweeping area
-## d/dt - derivative w.r.t. time
-## L - planet's angular momentum
-## m - planet's mass
+planet_angular_momentum = symbols.angular_momentum
+"""
+:symbols:`angular_momentum` of the planet.
+"""
 
-time = Symbol("time", units.time)
-area_swept = Function("area_swept", units.area)
-planet_angular_momentum = Symbol("planet_angular_momentum", units.length * units.momentum)
 planet_mass = symbols.mass
+"""
+:symbols:`mass` of the planet.
+"""
 
 law = Eq(Derivative(area_swept(time), time), planet_angular_momentum / (2 * planet_mass))
+"""
+:laws:symbol::
 
-# TODO: derive law from expression of area and definition of angular momentum
-
+:laws:latex::
+"""
 
-def print_law() -> str:
-    return print_expression(law)
+# TODO: derive law from expression of area and definition of angular momentum
 
 
 @validate_input(planet_angular_momentum_=planet_angular_momentum, planet_mass_=planet_mass)
diff --git a/symplyphysics/laws/gravity/corrected_planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py b/symplyphysics/laws/gravity/corrected_planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
index 28d812cbf..010980cc7 100644
--- a/symplyphysics/laws/gravity/corrected_planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
+++ b/symplyphysics/laws/gravity/corrected_planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
@@ -7,34 +7,34 @@
 the orbiting planet. If this assumption isn't made, the first and the second law stay the same, but
 the third law requires a revision. As a result, the period of the planet's rotation depends not only
 on the attracting mass, but also on the mass of the planet itself.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. Sivukhin D.V. (1979), *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 321, (59.3).
 """
 
 from sympy import Eq, pi, solve
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
     symbols,
     clone_as_symbol,
+    quantities,
 )
 
-rotation_period = Symbol("rotation_period", units.time)
+rotation_period = symbols.period
 """
-The planet's period of rotation.
-
-Symbol:
-    :code:`T`
+The planet's :symbols:`period` of rotation.
 """
 
-semimajor_axis = Symbol("semimajor_axis", units.length)
+semimajor_axis = symbols.semimajor_axis
 """
-The semi-major axis of the planet's orbit.
-
-Symbol:
-    :code:`a`
+The :symbols:`semimajor_axis` of the planet's orbit.
 """
 
 attracting_mass = clone_as_symbol(symbols.mass, display_symbol="M", display_latex="M")
@@ -49,14 +49,12 @@
 
 law = Eq(
     rotation_period**2,
-    (4 * pi**2 * semimajor_axis**3) / (gravitational_constant * (attracting_mass + planetary_mass)),
+    (4 * pi**2 * semimajor_axis**3) / (quantities.gravitational_constant * (attracting_mass + planetary_mass)),
 )
-r"""
-:code:`T^2 = (4 * pi^2 * a^3) / (G * (M + m))`
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        T^2 = \frac{4 \pi^2 a^3}{G \left( M + m \right)}
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/easterly_deviation_from_plumbline_of_falling_bodies.py b/symplyphysics/laws/gravity/easterly_deviation_from_plumbline_of_falling_bodies.py
index 4ad7aabd7..829fe3d13 100644
--- a/symplyphysics/laws/gravity/easterly_deviation_from_plumbline_of_falling_bodies.py
+++ b/symplyphysics/laws/gravity/easterly_deviation_from_plumbline_of_falling_bodies.py
@@ -1,39 +1,69 @@
+"""
+Easterly deviation from plumbline of falling bodies
+===================================================
+
+Suppose a body is falling freely in Earth's gravity field with its initial velocity being zero.
+Then the effect of the Coriolis force on the falling body can be found in the fact that it deflects
+from plumbline in the easterly and southerly (equatorial) directions.
+
+**Conditions:**
+
+#. The vector of free fall acceleration is constant.
+
+**Links:**
+
+#. Sivukhin D.V. (1979), *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 355, (67.10).
+"""
+
 from sympy import Eq, cos, pi
 from symplyphysics import (
-    units,
-    angle_type,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
+    symbols,
+    clone_as_symbol,
 )
 from symplyphysics.core.symbols.quantities import scale_factor
 
-# Description
-## Suppose a body is falling freely in Earth's gravity field with its initial velocity being zero.
-## Then the effect of the Coriolis force on the falling body can be found in the fact that it deflects
-## from plumbline in the easterly and southerly (equatorial) directions.
+easterly_deviation_from_plumbline = clone_as_symbol(
+    symbols.euclidean_distance,
+    display_symbol="s_east",
+    display_latex="s_\\text{east}",
+)
+"""
+Easterly deviation of falling body from plumbline due to Earth's rotation.
+See :symbols:`euclidean_distance`
+"""
+
+fall_time = symbols.time
+"""
+:symbols:`time` of the body's fall.
+"""
 
-# Conditions
-## - The vector of free fall acceleration `g` is considered constant.
+rotation_period = symbols.period
+"""
+:symbols:`period` of the Earth's rotation.
+"""
 
-# Law: s_east = (4 * pi / 3) * (t / T) * h * cos(theta)
-## s_east - easterly deviation of falling body from plumbline due to Earth's rotation
-## t - time of the body's fall
-## T - period of Earth's rotation
-## h - initial elevation of the body from Earth's surface
-## theta - latitude of the place the body is located in
+initial_elevation = symbols.height
+"""
+Initial elevation (:symbols:`height`) of the body from the Earth's surface.
+"""
 
-easterly_deviation_from_plumbline = Symbol("easterly_deviation_from_plumbline", units.length)
-fall_time = Symbol("fall_time", units.time)
-rotation_period = Symbol("rotation_period", units.time)
-initial_elevation = Symbol("initial_elevation", units.length)
-latitude = Symbol("latitude", angle_type)
+latitude = symbols.latitude
+"""
+:symbols:`latitude` of the location of the body.
+"""
 
 law = Eq(
     easterly_deviation_from_plumbline,
     (4 * pi / 3) * (fall_time / rotation_period) * initial_elevation * cos(latitude),
 )
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # TODO: derive from the solution of the relative motion equation in Earth's gravitational field.
 
diff --git a/symplyphysics/laws/gravity/escape_velocity.py b/symplyphysics/laws/gravity/escape_velocity.py
index 74a9f5906..26d9e578c 100644
--- a/symplyphysics/laws/gravity/escape_velocity.py
+++ b/symplyphysics/laws/gravity/escape_velocity.py
@@ -1,73 +1,96 @@
+"""
+First escape speed
+==================
+
+**First escape speed** is the minimum speed needed for an object to escape from contact with or orbit of a primary body.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Escape_velocity#From_an_orbiting_body>`__.
+
+..
+    TODO rename file
+"""
+
 from sympy import Eq, solve, sqrt
-from sympy.physics.units import gravitational_constant
-from symplyphysics import (units, Quantity, Symbol, print_expression, validate_input,
-    validate_output, symbols)
+from symplyphysics import Quantity, validate_input, validate_output, symbols, quantities
 from symplyphysics.core.expr_comparisons import expr_equals
-
 from symplyphysics.laws.gravity import gravity_force_from_mass_and_distance as gravity_force_law
 from symplyphysics.laws.dynamics import acceleration_is_force_over_mass as acceleration_law
 from symplyphysics.laws.kinematics import centripetal_acceleration_via_linear_speed_and_radius as centripetal_law
 
-# Law: V = √(G * M / (R + h))
-# Where:
-# V - initial velocity
-# M - planet_mass of planet
-# G - gravitational constant
-# h - height above the planet surface
-# R - radius of planet
-
-velocity = Symbol("initial_velocity", units.velocity)
-radius = Symbol("radius", units.length)
-height = Symbol("height", units.length)
+speed = symbols.speed
+"""
+Escape :symbols:`speed` of the body.
+"""
+
+radius = symbols.radius
+"""
+:symbols:`radius` of the planet.
+"""
+
+height = symbols.height
+"""
+Elevation (:symbols:`height`) of the body from the surface of the planet
+"""
+
 planet_mass = symbols.mass
+"""
+:symbols:`mass` of the planet.
+"""
 
-law = Eq(velocity, sqrt(gravitational_constant * planet_mass / (radius + height)))
+law = Eq(speed, sqrt(quantities.gravitational_constant * planet_mass / (radius + height)))
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # This law might be derived via "gravity_force_from_mass_and_distance" law, "acceleration_from_force" law,
 # "centripetal_acceleration_is_squared_velocity_by_radius" law.
 
 # The radius of the orbit consists of the radius of the planet and the height above its surface.
-centripetal_law_applied = centripetal_law.law.subs({
-    centripetal_law.speed: velocity,
+_centripetal_law_applied = centripetal_law.law.subs({
+    centripetal_law.speed: speed,
     centripetal_law.radius_of_curvature: radius + height,
 })
-acceleration_derived = solve(centripetal_law_applied,
+_acceleration_derived = solve(_centripetal_law_applied,
     centripetal_law.centripetal_acceleration,
     dict=True)[0][centripetal_law.centripetal_acceleration]
 
-acceleration_law_applied = acceleration_law.law.subs({
-    acceleration_law.symbols.acceleration: acceleration_derived,
+_acceleration_law_applied = acceleration_law.law.subs({
+    acceleration_law.symbols.acceleration: _acceleration_derived,
 })
-force_derived = solve(acceleration_law_applied, acceleration_law.symbols.force,
+_force_derived = solve(_acceleration_law_applied, acceleration_law.symbols.force,
     dict=True)[0][acceleration_law.symbols.force]
 
 # Let's write down Newton's second law: ma = F. F is, in this case, the force of gravity. And in the general case,
-# when a body moves along a circle with a constant velocity in modulus, its acceleration is equal to the centripetal
+# when a body moves along a circle with a constant speed in modulus, its acceleration is equal to the centripetal
 # acceleration.
-gravity_force_law_applied = gravity_force_law.law.subs({
+_gravity_force_law_applied = gravity_force_law.law.subs({
     gravity_force_law.first_mass: planet_mass,
-    gravity_force_law.gravitational_force: force_derived,
+    gravity_force_law.gravitational_force: _force_derived,
     gravity_force_law.second_mass: acceleration_law.symbols.mass,
     gravity_force_law.distance_between_mass_centers: radius + height,
 })
-# The first cosmic velocity is the minimum horizontal velocity that must be given to an object so that it moves in
-# a circular orbit around the planet. Based on this definition, the velocity could be left negative, implying a different
+# The first cosmic speed is the minimum horizontal speed that must be given to an object so that it moves in
+# a circular orbit around the planet. Based on this definition, the speed could be left negative, implying a different
 # direction. But most often they are written with a plus sign, implying the modulus of the horizontal component of the
-# velocity tangent to the orbit. Therefore, the first solution that returns a minus is ignored.
-velocity_derived = solve(gravity_force_law_applied, velocity, dict=True)[1][velocity]
-
-# Check if derived velocity is same as declared.
-assert expr_equals(velocity_derived, law.rhs)
-
+# speed tangent to the orbit. Therefore, the first solution that returns a minus is ignored.
+_velocity_derived = solve(_gravity_force_law_applied, speed, dict=True)[1][speed]
 
-def print_law() -> str:
-    return print_expression(law)
+# Check if derived speed is same as declared.
+assert expr_equals(_velocity_derived, law.rhs)
 
 
 @validate_input(planet_mass_=planet_mass, radius_=radius, height_=height)
-@validate_output(velocity)
+@validate_output(speed)
 def calculate_velocity(planet_mass_: Quantity, radius_: Quantity, height_: Quantity) -> Quantity:
-    result_velocity_expr = solve(law, velocity, dict=True)[0][velocity]
+    result_velocity_expr = solve(law, speed, dict=True)[0][speed]
     result_expr = result_velocity_expr.subs({
         planet_mass: planet_mass_,
         radius: radius_,
diff --git a/symplyphysics/laws/gravity/free_fall_acceleration_from_height.py b/symplyphysics/laws/gravity/free_fall_acceleration_from_height.py
index 140fb20b5..834e1c99a 100644
--- a/symplyphysics/laws/gravity/free_fall_acceleration_from_height.py
+++ b/symplyphysics/laws/gravity/free_fall_acceleration_from_height.py
@@ -1,52 +1,68 @@
-from sympy import (Eq, solve)
-from sympy.physics.units import gravitational_constant
-from symplyphysics import (units, Quantity, Symbol, print_expression, validate_input,
-    validate_output, symbols)
+"""
+Free fall acceleration from height
+==================================
+
+**Free fall** is any motion of a body where gravity is the only force acting upon it.
+**Free fall acceleration** is the acceleration the body experiences during the free fall.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+"""
+
+from sympy import Eq, solve
+from symplyphysics import Quantity, validate_input, validate_output, symbols, quantities
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.laws.gravity import gravity_force_from_mass_and_distance as gravity_law
 from symplyphysics.laws.dynamics import acceleration_is_force_over_mass as newton2_law
 
-# Description
-## Every planet generates gravity field which causes free falling. Free fall acceleration depends on height above the planet surface.
-## Law: g = G * M / (R + h)**2
-## Where:
-## g is free fall acceleration
-## G is gravitational constant
-## M is mass of the planet
-## R is radius of the planet
-## h is height above the planet surface
-
 free_fall_acceleration = symbols.acceleration
-planet_radius = Symbol("planet_radius", units.length)
-height_above_surface = Symbol("height_above_surface", units.length)
+"""
+Free fall :symbols:`acceleration` of the body.
+"""
+
+planet_radius = symbols.radius
+"""
+:symbols:`radius` of the planet.
+"""
+
+elevation = symbols.height
+"""
+Elevation (:symbols:`height`) of the body from the planet's surface.
+"""
+
 planet_mass = symbols.mass
+"""
+:symbols:`mass` of the planet.
+"""
 
 law = Eq(free_fall_acceleration,
-    gravitational_constant * planet_mass / (planet_radius + height_above_surface)**2)
+    quantities.gravitational_constant * planet_mass / (planet_radius + elevation)**2)
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # This law might be easily derived from gravitational law via Newton's law #2
 ## Distance between mass centers is radius of the planet plus height above it's surface.
-gravitational_force = gravity_law.law.rhs.subs({
+_gravitational_force = gravity_law.law.rhs.subs({
     gravity_law.first_mass: planet_mass,
-    gravity_law.distance_between_mass_centers: planet_radius + height_above_surface
+    gravity_law.distance_between_mass_centers: planet_radius + elevation
 })
 
 # Substitute mass first
-derived_free_fall_acceleration = newton2_law.law.rhs.subs(newton2_law.mass, gravity_law.second_mass)
-derived_free_fall_acceleration = derived_free_fall_acceleration.subs(newton2_law.force,
-    gravitational_force)
+_derived_free_fall_acceleration = newton2_law.law.rhs.subs(newton2_law.mass, gravity_law.second_mass)
+_derived_free_fall_acceleration = _derived_free_fall_acceleration.subs(newton2_law.force,
+    _gravitational_force)
 
 # Check if derived acceleration is same as declared
-assert expr_equals(derived_free_fall_acceleration, law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_derived_free_fall_acceleration, law.rhs)
 
 
 @validate_input(planet_mass_=planet_mass,
     planet_radius_=planet_radius,
-    height_above_surface_=height_above_surface)
+    height_above_surface_=elevation)
 @validate_output(free_fall_acceleration)
 def calculate_acceleration(planet_mass_: Quantity, planet_radius_: Quantity,
     height_above_surface_: Quantity) -> Quantity:
@@ -54,6 +70,6 @@ def calculate_acceleration(planet_mass_: Quantity, planet_radius_: Quantity,
     result_expr = result_accel_expr.subs({
         planet_mass: planet_mass_,
         planet_radius: planet_radius_,
-        height_above_surface: height_above_surface_
+        elevation: height_above_surface_
     })
     return Quantity(result_expr)
diff --git a/symplyphysics/laws/gravity/gravitational_potential_energy.py b/symplyphysics/laws/gravity/gravitational_potential_energy.py
index 9eea083ca..e0f0362d7 100644
--- a/symplyphysics/laws/gravity/gravitational_potential_energy.py
+++ b/symplyphysics/laws/gravity/gravitational_potential_energy.py
@@ -1,39 +1,58 @@
+"""
+Gravitational potential energy
+==============================
+
+The **gravitational potential energy** of a system of two particles is the negative of the work
+that would be done by the gravitational force of either particle acting on the other if the
+particles were brought together from infinity to the given distance.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. `Physics LibreTexts. Gravitational Potential Energy beyond Earth (13.4.6) <https://phys.libretexts.org/Workbench/PH_245_Textbook_V2/13%3A_Gravitation/13.04%3A_Gravitational_Potential_Energy_and_Total_Energy>`__.
+"""
+
 from sympy import Eq
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
-    print_expression,
     validate_input,
     validate_output,
     clone_as_symbol,
     symbols,
+    quantities,
 )
 
-# Description
-## The gravitational potential energy U(r) of a system of two particles, with masses m1 and m2
-## and separated by a distance r, is the negative of the work that would be done by the gravitational
-## force of either particle acting on the other if the separation between the particles were changed
-## from infinite to r.
+gravitational_potential_energy = symbols.potential_energy
+"""
+Gravitational :symbols:`potential_energy`.
+"""
 
-# Law: U = -G * m1 * m2 / r
-## U - gravitational potential energy
-## G - gravitational constant
-## m1, m2 - masses of interacting particles
-## r - distance between mass centers of particles
+first_mass = clone_as_symbol(symbols.mass, subscript="1")
+"""
+:symbols:`mass` of the first particle.
+"""
 
-gravitational_potential_energy = Symbol("gravitational_potential_energy", units.energy)
-first_mass = clone_as_symbol(symbols.mass, display_symbol="m_1", display_latex="m_1")
-second_mass = clone_as_symbol(symbols.mass, display_symbol="m_2", display_latex="m_2")
-distance_between_mass_centers = Symbol("distance_between_mass_centers", units.length)
+second_mass = clone_as_symbol(symbols.mass, subscript="2")
+"""
+:symbols:`mass` of the second particle.
+"""
 
-law = Eq(gravitational_potential_energy,
-    -1 * gravitational_constant * first_mass * second_mass / distance_between_mass_centers)
+distance_between_mass_centers = symbols.euclidean_distance
+"""
+:symbols:`euclidean_distance` between mass centers of the particles.
+"""
 
 
-def print_law() -> str:
-    return print_expression(law)
+law = Eq(gravitational_potential_energy,
+    -1 * quantities.gravitational_constant * first_mass * second_mass / distance_between_mass_centers)
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 
 @validate_input(
diff --git a/symplyphysics/laws/gravity/gravitational_radius_of_body_mass.py b/symplyphysics/laws/gravity/gravitational_radius_of_body_mass.py
index f646543c7..6dd730e40 100644
--- a/symplyphysics/laws/gravity/gravitational_radius_of_body_mass.py
+++ b/symplyphysics/laws/gravity/gravitational_radius_of_body_mass.py
@@ -1,38 +1,53 @@
-from sympy import (
-    Eq,
-    solve,
-)
-from sympy.physics.units import gravitational_constant, speed_of_light
+"""
+Gravitational radius of massive body
+====================================
+
+The **gravitational radius**, or **Schwarzschild radius**, is a characteristic radius
+defined for any physical body with mass. This is the radius of the sphere on which the
+event horizon created by this mass would be located (from the point of view of general
+theory of relativity) if it were distributed spherically symmetrically, would be
+stationary (in particular, it would not rotate, but radial movements are permissible)
+and would lie entirely inside this sphere.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+#. :quantity_notation:`speed_of_light`.
+
+**Links:**
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Schwarzschild_radius>`__.
+
+..
+    TODO rename file
+"""
+
+from sympy import Eq, solve
 from symplyphysics import (
     symbols,
-    units,
     Quantity,
-    Symbol,
-    print_expression,
     validate_input,
     validate_output,
+    quantities,
 )
 
-# Description
-## The gravitational radius is a characteristic radius defined for any physical body with mass. This is the radius of the sphere
-## on which the event horizon created by this mass would be located (from the point of view of general theory of relativity) if it were distributed spherically
-## symmetrically, would be stationary (in particular, it would not rotate, but radial movements are permissible) and would lie entirely
-## inside this sphere.
+gravitational_radius = symbols.radius
+"""
+Gravitational :symbols:`radius` of the body.
+"""
 
-## Law is: R = 2 * G * M / c^2, where
-## R - gravitational radius of the body,
-## G - gravitational constant,
-## M - mass of the body,
-## c - speed of light.
-
-gravitational_radius = Symbol("gravitational_radius", units.length)
 body_mass = symbols.mass
+"""
+:symbols:`mass` of the body.
+"""
 
-law = Eq(gravitational_radius, 2 * gravitational_constant * body_mass / speed_of_light**2)
-
+law = Eq(gravitational_radius, 2 * quantities.gravitational_constant * body_mass / quantities.speed_of_light**2)
+"""
+:laws:symbol::
 
-def print_law() -> str:
-    return print_expression(law)
+:laws:latex::
+"""
 
 
 @validate_input(body_mass_=body_mass)
diff --git a/symplyphysics/laws/gravity/gravity_force_from_mass_and_distance.py b/symplyphysics/laws/gravity/gravity_force_from_mass_and_distance.py
index 65eb394d5..4642b0d1c 100644
--- a/symplyphysics/laws/gravity/gravity_force_from_mass_and_distance.py
+++ b/symplyphysics/laws/gravity/gravity_force_from_mass_and_distance.py
@@ -1,15 +1,29 @@
+"""
+Gravitational force from mass and distance
+==========================================
+
+In classical mechanics, the **gravitational force** is a fundamental attractive force that
+exists between any two massive bodies. Its magnitude is proportional to the mass of the
+bodies and inversely proportional to the distance squared between the bodies.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. `Physics LibreTexts. Newton's Law of Universal Gravitation (13.2.1) <https://phys.libretexts.org/Workbench/PH_245_Textbook_V2/13%3A_Gravitation/13.02%3A_Newton's_Law_of_Universal_Gravitation>`__.
+"""
+
 from sympy import Eq, solve, sqrt
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
-    print_expression,
     validate_input,
     validate_output,
     vector_magnitude,
     clone_as_symbol,
     symbols,
+    quantities,
 )
 from symplyphysics.core.dimensions import ScalarValue
 from symplyphysics.core.points.cartesian_point import CartesianPoint
@@ -19,28 +33,38 @@
     conservative_force_is_gradient_of_potential_energy as gradient_law,)
 from symplyphysics.laws.gravity import gravitational_potential_energy
 
-# Description
-## Every object generates gravity field around it. Any other object in this field is pulled toward generator.
-## Gravitational force between two object if proportional to masses of objects and counter-proportional to distance between their mass centers.
-## Law: F = G * m1 * m2 / R**2
-## Where:
-## F - gravitational force
-## m1 and m2 - masses of objects
-## R - distance between mass centers of objects
-## G - gravitational constant
-
 gravitational_force = symbols.force
-first_mass = clone_as_symbol(symbols.mass, display_symbol="m_1", display_latex="m_1")
-second_mass = clone_as_symbol(symbols.mass, display_symbol="m_2", display_latex="m_2")
-distance_between_mass_centers = Symbol("distance_between_mass_centers", units.length)
+"""
+Gravitational :symbols:`force`.
+"""
+
+first_mass = clone_as_symbol(symbols.mass, subscript="1")
+"""
+:symbols:`mass` of the first body.
+"""
+
+second_mass = clone_as_symbol(symbols.mass, subscript="2")
+"""
+:symbols:`mass` of the second body.
+"""
+
+distance_between_mass_centers = symbols.euclidean_distance
+"""
+:symbols:`euclidean_distance` between the centers of the bodies..
+"""
 
 law = Eq(gravitational_force,
-    gravitational_constant * first_mass * second_mass / distance_between_mass_centers**2)
+    quantities.gravitational_constant * first_mass * second_mass / distance_between_mass_centers**2)
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # Derive law from the gravitational potential energy
 # Condition: space must be 3-dimensional and flat
 
-potential = gravitational_potential_energy.law.rhs.subs({
+_potential = gravitational_potential_energy.law.rhs.subs({
     gravitational_potential_energy.first_mass: first_mass,
     gravitational_potential_energy.second_mass: second_mass,
     gravitational_potential_energy.distance_between_mass_centers: distance_between_mass_centers,
@@ -48,27 +72,23 @@
 
 
 def potential_field_function(point: CartesianPoint) -> ScalarValue:
-    return potential.subs(
+    return _potential.subs(
         distance_between_mass_centers,
         sqrt(point.x**2 + point.y**2 + point.z**2),
     )
 
 
-potential_field = ScalarField(potential_field_function)
+_potential_field = ScalarField(potential_field_function)
 
-gravitational_force_vector = gradient_law.law(potential_field)
-gravitational_force_derived = vector_magnitude(gravitational_force_vector).simplify()
+_gravitational_force_vector = gradient_law.law(_potential_field)
+_gravitational_force_derived = vector_magnitude(_gravitational_force_vector).simplify()
 
-x, y, z = gravitational_force_vector.coordinate_system.coord_system.base_scalars()
-gravitational_force_from_law = law.rhs.subs(distance_between_mass_centers, sqrt(x**2 + y**2 + z**2))
+_x, _y, _z = _gravitational_force_vector.coordinate_system.coord_system.base_scalars()
+_gravitational_force_from_law = law.rhs.subs(distance_between_mass_centers, sqrt(_x**2 + _y**2 + _z**2))
 
 # sympy avoids oversimplifications in case of square roots without certain assumptions,
 # therefore we resort to squaring both sides to make it work
-assert expr_equals(gravitational_force_derived**2, gravitational_force_from_law**2)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_gravitational_force_derived**2, _gravitational_force_from_law**2)
 
 
 @validate_input(first_object_mass_=first_mass,
diff --git a/symplyphysics/laws/gravity/keplers_constant_via_attracting_body_mass.py b/symplyphysics/laws/gravity/keplers_constant_via_attracting_body_mass.py
index 285663240..e8090eb25 100644
--- a/symplyphysics/laws/gravity/keplers_constant_via_attracting_body_mass.py
+++ b/symplyphysics/laws/gravity/keplers_constant_via_attracting_body_mass.py
@@ -5,27 +5,25 @@
 *Kepler's constant* is a physical quantity that only depends on the gravitational constant and the mass of
 the orbited body, such as the Sun. It is constant in the sense that all planets that orbit the same body
 approximately have the same value of the Kepler's constant.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
 """
 
 from sympy import Eq, pi
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
     symbols,
     clone_as_symbol,
+    quantities,
 )
 
-keplers_constant = Symbol("keplers_constant", units.length**3 / units.time**2)
-r"""
-Symbol:
-    K
-
-Latex:
-    :math:`\mathfrak{K}`
+kepler_constant = symbols.kepler_constant
+"""
+:symbols:`kepler_constant` of the system.
 """
 
 attracting_mass = clone_as_symbol(symbols.mass, display_symbol="M", display_latex="M")
@@ -33,17 +31,16 @@
 The :symbols:`mass` of the attracting body.
 """
 
-law = Eq(keplers_constant, gravitational_constant * attracting_mass / (4 * pi**2))
-r"""
-K = G * M / (4 * pi^2)
+law = Eq(kepler_constant, quantities.gravitational_constant * attracting_mass / (4 * pi**2))
+"""
+:laws:symbol::
 
-Latex:
-    :math:`\mathfrak{K} = \frac{G M}{4 \pi^2}`
+:laws:latex::
 """
 
 
 @validate_input(attracting_body_mass_=attracting_mass)
-@validate_output(keplers_constant)
+@validate_output(kepler_constant)
 def calculate_keplers_constant(attracting_body_mass_: Quantity) -> Quantity:
     result = law.rhs.subs({
         attracting_mass: attracting_body_mass_,
diff --git a/symplyphysics/laws/gravity/maximum_height_of_a_body_thrown_at_an_angle_to_horizon.py b/symplyphysics/laws/gravity/maximum_height_of_a_body_thrown_at_an_angle_to_horizon.py
index 1d1b46c8b..3ab419e2e 100644
--- a/symplyphysics/laws/gravity/maximum_height_of_a_body_thrown_at_an_angle_to_horizon.py
+++ b/symplyphysics/laws/gravity/maximum_height_of_a_body_thrown_at_an_angle_to_horizon.py
@@ -1,74 +1,106 @@
-from sympy import (Eq, solve, sin, pi)
-from sympy.physics.units import acceleration_due_to_gravity as earth_free_fall_acceleration
-from symplyphysics import (units, Quantity, Symbol, print_expression, validate_input,
-    validate_output, angle_type)
+r"""
+Maximum height of body thrown at angle to horizon
+=================================================
+
+The maximum height of a projectile above its launch position is a function of
+the vertical component of the initial velocity, the angle between the initial
+velocity and the horizon, and the magnitude of the free fall acceleration.
+
+**Notation:**
+
+#. :quantity_notation:`acceleration_due_to_gravity`.
+
+**Conditions:**
+
+#. :math:`\mathbf{v}_0 \cdot \mathbf{g} < 0`, i.e. the vector of initial velocity
+   :math:`\mathbf{v}_0` must have a non-zero component that is antiparallel to the
+   vector of free fall acceleration :math:`\mathbf{g}`.
+
+#. The height is measured with respect to the horizontal plane where the projectile
+   was located at initial time.
+
+**Links:**
+
+#. `Physics LibreTexts. Projectile Motion, Maximum Height (3.3.14) <https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion>`__.
+"""
+
+from sympy import Eq, solve, sin, pi
+from symplyphysics import (
+    Quantity,
+    validate_input,
+    validate_output,
+    symbols,
+    quantities,
+    convert_to_si,
+    clone_as_symbol,
+)
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.laws.kinematics import position_via_constant_acceleration_and_time as distance_law
 from symplyphysics.laws.kinematics import speed_via_constant_acceleration_and_time as velocity_law
 from symplyphysics.laws.geometry import planar_projection_is_cosine as projection_law
 
-# Description
-## Let's say we throw the body at an angle to the horizon with some initial velocity.
-## Then the height of throw of the body will depend only on the initial velocity,
-## the angle of the throw and the acceleration of free fall.
-## The angle of the throw is the angle between the initial velocity vector and the horizontal axis.
+height = symbols.height
+"""
+Maximum :symbols:`height` of the projectile.
+"""
 
-## Law is: h = v0^2 * sin(a)^2 / (2 * g), where
-## h - height of throw,
-## v0 - initial velocity,
-## a - angle of throw,
-## g - acceleration of free fall.
+initial_speed = clone_as_symbol(symbols.speed, subscript="0")
+"""
+Initial :symbols:`speed` of the body.
+"""
 
-height = Symbol("height", units.length)
+angle = symbols.angle
+"""
+:symbols:`angle` between the initial velocity and the horizon.
+"""
 
-initial_velocity = Symbol("initial_velocity", units.velocity)
-angle = Symbol("angle", angle_type)
+law = Eq(height, initial_speed**2 * sin(angle)**2 / (2 * quantities.acceleration_due_to_gravity))
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
-law = Eq(height, initial_velocity**2 * sin(angle)**2 / (2 * earth_free_fall_acceleration))
 
 # This law might be derived via "constant_acceleration_movement_is_parabolic" law, "planar_projection_is_cosine" law
 # and "accelerated_velocity_from_time" law.
 
 # The law seeks a projection on the horizontal axis, but a projection on the vertical axis is necessary,
 # so the angle is represented as a "pi/2 - angle".
-projection_law_applied = projection_law.law.subs({
-    projection_law.vector_length: initial_velocity,
+_projection_law_applied = projection_law.law.subs({
+    projection_law.vector_length: initial_speed,
     projection_law.vector_angle: (pi / 2) - angle,
 })
-vertical_projection_derived = solve(projection_law_applied, projection_law.projection,
+_vertical_projection_derived = solve(_projection_law_applied, projection_law.projection,
     dict=True)[0][projection_law.projection]
 
 # Vertical velocity is zero in the highest point of trajectory.
-velocity_law_applied = velocity_law.law.subs({
-    velocity_law.initial_speed: vertical_projection_derived,
+_velocity_law_applied = velocity_law.law.subs({
+    velocity_law.initial_speed: _vertical_projection_derived,
     velocity_law.final_speed: 0,
-    velocity_law.acceleration: -earth_free_fall_acceleration,
+    velocity_law.acceleration: -quantities.acceleration_due_to_gravity,
 })
-time_derived = solve(velocity_law_applied, velocity_law.time, dict=True)[0][velocity_law.time]
+_time_derived = solve(_velocity_law_applied, velocity_law.time, dict=True)[0][velocity_law.time]
 
 # The acceleration of gravity is directed opposite to the vertical coordinate axis,
 ## so there is a minus sign before the acceleration.
-height_law_applied = distance_law.law.subs({
-    distance_law.initial_speed: vertical_projection_derived,
-    distance_law.time: time_derived,
-    distance_law.acceleration: -earth_free_fall_acceleration,
+_height_law_applied = distance_law.law.subs({
+    distance_law.initial_speed: _vertical_projection_derived,
+    distance_law.time: _time_derived,
+    distance_law.acceleration: -quantities.acceleration_due_to_gravity,
     distance_law.initial_position: 0,
 })
 
 # Check if derived height is same as declared.
-assert expr_equals(height_law_applied.rhs, law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_height_law_applied.rhs, law.rhs)
 
 
-@validate_input(initial_velocity_=initial_velocity, angle_=angle)
+@validate_input(initial_velocity_=initial_speed, angle_=angle)
 @validate_output(height)
 def calculate_height(initial_velocity_: Quantity, angle_: float | Quantity) -> Quantity:
     result_expr = solve(law, height, dict=True)[0][height]
     result_expr = result_expr.subs({
-        initial_velocity: initial_velocity_,
-        angle: angle_,
+        initial_speed: initial_velocity_,
+        angle: convert_to_si(angle_),
     })
     return Quantity(result_expr)
diff --git a/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon.py b/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon.py
index ef602f92d..0bae633a5 100644
--- a/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon.py
+++ b/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon.py
@@ -1,66 +1,87 @@
-from sympy import (Eq, solve, sin, pi)
-from sympy.physics.units import acceleration_due_to_gravity as earth_free_fall_acceleration
-from symplyphysics import (units, Quantity, Symbol, print_expression, validate_input,
-    validate_output, angle_type)
+"""
+Time of flight of a projectile via initial velocity
+===================================================
+
+Let the body be thrown at an angle to the horizon with some initial velocity. Then the time of the motion
+of this body until it reaches the initial elevation depends only on initial velocity and acceleration
+of free fall.
+
+**Conditions:**
+
+#. The acceleration of gravity stays constant throughout the movement of the body.
+
+#. The initial and final height of the body with respect to the "ground" are the same,
+   i.e. the body lands on a point with the same vertical coordinate as at the start
+   of the flight.
+
+**Links:**
+
+#. `Physics LibreTexts. Projectile Motion, Maximum Height (3.3.13) <https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion>`__.
+
+..
+    TODO rename file
+"""
+
+from sympy import Eq, solve, sin, pi
+from symplyphysics import Quantity, validate_input, validate_output, symbols, quantities, clone_as_symbol
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.laws.kinematics import position_via_constant_acceleration_and_time as distance_law
 from symplyphysics.laws.geometry import planar_projection_is_cosine as projection_law
 
-# Description
-## Let's say we throw the body at an angle to the horizon with some initial velocity.
-## Then the time of movement (time of fall) of the body will depend only on the initial velocity,
-## the angle of the throw and the acceleration of free fall.
-## The angle of the throw is the angle between the initial velocity vector and the horizontal axis.
+time = symbols.time
+"""
+:symbols:`time` of flight of the projectile.
+"""
 
-## Law is: t = 2 * v0 * sin(a) / g, where
-## t - movement time,
-## v0 - initial velocity,
-## a - angle of throw,
-## g - acceleration of free fall.
+initial_speed = clone_as_symbol(symbols.speed, subscript="0")
+"""
+Initial :symbols:`speed` of the projectile.
+"""
 
-movement_time = Symbol("movement_time", units.time)
+angle = symbols.angle
+"""
+:symbols:`angle`
+"""
 
-initial_velocity = Symbol("initial_velocity", units.velocity)
-angle = Symbol("angle", angle_type)
+law = Eq(time, 2 * initial_speed * sin(angle) / quantities.acceleration_due_to_gravity)
+"""
+:laws:symbol::
 
-law = Eq(movement_time, 2 * initial_velocity * sin(angle) / earth_free_fall_acceleration)
+:laws:latex::
+"""
 
 # This law might be derived via "constant_acceleration_movement_is_parabolic" law
 # and "planar_projection_is_cosine" law.
 
 # The law seeks a projection on the horizontal axis, but a projection on the vertical axis is necessary,
 # so the angle is represented as a "pi/2 - angle".
-projection_law_applied = projection_law.law.subs({
-    projection_law.vector_length: initial_velocity,
+_projection_law_applied = projection_law.law.subs({
+    projection_law.vector_length: initial_speed,
     projection_law.vector_angle: (pi / 2) - angle,
 })
-projection_derived = solve(projection_law_applied, projection_law.projection,
+_projection_derived = solve(_projection_law_applied, projection_law.projection,
     dict=True)[0][projection_law.projection]
 
 # The acceleration of gravity is directed opposite to the vertical coordinate axis,
 ## so there is a minus sign before the acceleration.
-distance_law_applied = distance_law.law.subs({
-    distance_law.initial_speed: projection_derived,
-    distance_law.acceleration: -earth_free_fall_acceleration,
+_distance_law_applied = distance_law.law.subs({
+    distance_law.initial_speed: _projection_derived,
+    distance_law.acceleration: -quantities.acceleration_due_to_gravity,
     distance_law.initial_position: 0,
     distance_law.final_position: 0,
 })
-time_derived = solve(distance_law_applied, distance_law.time, dict=True)[1][distance_law.time]
+_time_derived = solve(_distance_law_applied, distance_law.time, dict=True)[1][distance_law.time]
 
 # Check if derived movement time is same as declared.
-assert expr_equals(time_derived, law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_time_derived, law.rhs)
 
 
-@validate_input(initial_velocity_=initial_velocity, angle_=angle)
-@validate_output(movement_time)
+@validate_input(initial_velocity_=initial_speed, angle_=angle)
+@validate_output(time)
 def calculate_movement_time(initial_velocity_: Quantity, angle_: float | Quantity) -> Quantity:
-    result_expr = solve(law, movement_time, dict=True)[0][movement_time]
+    result_expr = solve(law, time, dict=True)[0][time]
     result_expr = result_expr.subs({
-        initial_velocity: initial_velocity_,
+        initial_speed: initial_velocity_,
         angle: angle_,
     })
     return Quantity(result_expr)
diff --git a/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_horizontally.py b/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_horizontally.py
index a59ba66ac..45183d27b 100644
--- a/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_horizontally.py
+++ b/symplyphysics/laws/gravity/maximum_movement_time_of_a_body_thrown_horizontally.py
@@ -1,54 +1,74 @@
+"""
+Time of flight of a projectile via maximum height
+=================================================
+
+Let the body be thrown horizontally with some initial velocity. Then the time of the motion
+of this body until it reaches the initial elevation depends only on the height and acceleration
+of free fall.
+
+**Conditions:**
+
+#. The acceleration of gravity stays constant throughout the movement of the body.
+
+#. The initial and final height of the body with respect to the "ground" are the same,
+   i.e. the body lands on a point with the same vertical coordinate as at the start
+   of the flight.
+
+**Links:**
+
+#. `Physics LibreTexts. Zero Launch Angle, Duration of Flight (3.3.22) <https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion>`__.
+
+..
+    TODO rename file
+"""
+
 from sympy import (Eq, solve, sqrt)
-from sympy.physics.units import acceleration_due_to_gravity as earth_free_fall_acceleration
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
-    print_expression,
     validate_input,
     validate_output,
+    quantities,
+    symbols,
 )
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.laws.kinematics import position_via_constant_acceleration_and_time as distance_law
 
-# Description
-## Let the body be thrown horizontally with some initial velocity. Then the time of motion (time of fall)
-## of this body depends only on the height and acceleration of free fall.
-
-## Law is: t = sqrt(2 * h / g), where
-## t - movement time,
-## h - height,
-## g - acceleration of free fall.
+time = symbols.time
+"""
+:symbols:`time` of flight of the projectile.
+"""
 
-movement_time = Symbol("movement_time", units.time)
+height = symbols.height
+"""
+Maximum :symbols:`height` which the projectile reaches during its motion.
+"""
 
-height = Symbol("height", units.length)
+law = Eq(time, sqrt(2 * height / quantities.acceleration_due_to_gravity))
+"""
+:laws:symbol::
 
-law = Eq(movement_time, sqrt(2 * height / earth_free_fall_acceleration))
+:laws:latex::
+"""
 
 # This law might be derived via "constant_acceleration_movement_is_parabolic" law.
 # Horizontal vector of movement does not change falling time.
 
-distance_law_applied = distance_law.law.subs({
+_distance_law_applied = distance_law.law.subs({
     distance_law.initial_speed: 0,
-    distance_law.acceleration: earth_free_fall_acceleration,
+    distance_law.acceleration: quantities.acceleration_due_to_gravity,
     distance_law.initial_position: 0,
     distance_law.final_position: height,
 })
-time_derived = solve(distance_law_applied, distance_law.time, dict=True)[1][distance_law.time]
+_time_derived = solve(_distance_law_applied, distance_law.time, dict=True)[1][distance_law.time]
 
 # Check if derived movement time is same as declared.
-assert expr_equals(time_derived, law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_time_derived, law.rhs)
 
 
 @validate_input(height_=height)
-@validate_output(movement_time)
+@validate_output(time)
 def calculate_movement_time(height_: Quantity) -> Quantity:
-    result_expr = solve(law, movement_time, dict=True)[0][movement_time]
+    result_expr = solve(law, time, dict=True)[0][time]
     result_expr = result_expr.subs({
         height: height_,
     })
diff --git a/symplyphysics/laws/gravity/planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py b/symplyphysics/laws/gravity/planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
index 34eac5b6a..e056d1b76 100644
--- a/symplyphysics/laws/gravity/planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
+++ b/symplyphysics/laws/gravity/planet_period_squared_is_proportional_to_cube_of_semimajor_axis.py
@@ -1,14 +1,23 @@
+"""
+Planet period squared is proportional to cube of semimajor axis
+===============================================================
+
+Also known as **Kepler's third law** of planetary motion, the law of periods relates
+the period of rotation of any planet to the semi-major axis of its orbit.
+
+**Links:**
+
+#. `Physics LibreTexts. Kepler's Third Law, Derivation of Kepler's Third Law (5.6.23) <https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/5%3A_Uniform_Circular_Motion_and_Gravitation/5.6%3A_Keplers_Laws>`__.
+"""
+
 from sympy import Eq, solve, pi, Symbol as SymSymbol
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
-    print_expression,
     validate_input,
     validate_output,
     clone_as_symbol,
     symbols,
+    quantities,
 )
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.definitions import (
@@ -20,24 +29,31 @@
 from symplyphysics.laws.kinematics import (
     centripetal_acceleration_via_angular_speed_and_radius as centripetal_law,)
 
-# Description
-## Also known as Kepler's third law of planetary motion, the law of periods relates
-## the period of rotation of any planet to the semi-major axis of its orbit.
-
-# Law: T**2 = 4*pi**2 / (G*M) * a**3
-## T - planet's period of rotation
-## G - gravitational constant
-## M - mass of the attracting body, the Sun in case of the solar system
-## a - semi-major axis of the planet's orbit, its radius in case of a round orbit
+rotation_period = clone_as_symbol(symbols.period, positive=True)
+"""
+The planet's :symbols:`period` of rotation.
+"""
 
-rotation_period = Symbol("rotation_period", units.time, positive=True)
 attracting_mass = clone_as_symbol(symbols.mass, positive=True)
-semimajor_axis = Symbol("semimajor_axis", units.length, positive=True)
+"""
+:symbols:`mass` of the attracting body, e.g. the Sun in case of the solar system.
+"""
+
+semimajor_axis = clone_as_symbol(symbols.semimajor_axis, positive=True)
+"""
+:symbols:`semimajor_axis` of the planet's orbit. It is equal to the :symbols:`radius` in case of
+a round orbit.
+"""
 
 law = Eq(
     rotation_period**2,
-    4 * pi**2 / (gravitational_constant * attracting_mass) * semimajor_axis**3,
+    4 * pi**2 / (quantities.gravitational_constant * attracting_mass) * semimajor_axis**3,
 )
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # Derive law from Newton's second law of motion
 
@@ -75,10 +91,6 @@
 assert expr_equals(_period_derived, _period_from_law)
 
 
-def print_law() -> str:
-    return print_expression(law)
-
-
 @validate_input(attracting_mass_=attracting_mass, semimajor_axis_=semimajor_axis)
 @validate_output(rotation_period)
 def calculate_rotation_period(
diff --git a/symplyphysics/laws/gravity/radial_motion/average_potential_energy_via_average_kinetic_energy.py b/symplyphysics/laws/gravity/radial_motion/average_potential_energy_via_average_kinetic_energy.py
index 2f363dfd9..a31d66e00 100644
--- a/symplyphysics/laws/gravity/radial_motion/average_potential_energy_via_average_kinetic_energy.py
+++ b/symplyphysics/laws/gravity/radial_motion/average_potential_energy_via_average_kinetic_energy.py
@@ -1,4 +1,4 @@
-r"""
+"""
 Average potential energy via average kinetic energy
 ===================================================
 
@@ -7,47 +7,45 @@
 
 **Conditions:**
 
-#. Works for elliptical orbits, i.e. the total energy :math:`E` of the planet is negative.
+#. Works for elliptical orbits, in which case the total energy :math:`E` of the planet is negative.
+
+**Links:**
+
+#. `Physics LibreTexts. Energy in Circular Orbits <https://phys.libretexts.org/Workbench/PH_245_Textbook_V2/13%3A_Gravitation/13.05%3A_Satellite_Orbits_and_Energy>`__.
 """
 
 from sympy import Eq
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
+    symbols,
+    clone_as_symbol,
 )
 
-average_potential_energy = Symbol("average_potential_energy", units.energy)
-r"""
-The potential energy of the planet averaged over time.
-
-Symbol:
-    :code:`avg(U)`
-
-Latex:
-    :math:`\langle U \rangle`
+average_potential_energy = clone_as_symbol(
+    symbols.potential_energy,
+    display_symbol="avg(U)",
+    display_latex="\\langle U \\rangle",
+)
+"""
+The :symbols:`potential_energy` of the planet averaged over :symbols:`time`.
 """
 
-average_kinetic_energy = Symbol("average_kinetic_energy", units.energy)
-r"""
-The kinetic energy of the planet averaged over time.
-
-Symbol:
-    :code:`avg(K)`
-
-Latex:
-    :math:`\langle K \rangle`
+average_kinetic_energy = clone_as_symbol(
+    symbols.kinetic_energy,
+    display_symbol="avg(K)",
+    display_latex="\\langle K \\rangle",
+)
+"""
+The :symbols:`kinetic_energy` of the planet averaged over :symbols:`time`.
 """
 
 law = Eq(average_potential_energy, -2 * average_kinetic_energy)
-r"""
-:code:`avg(U) = -2 * avg(K)`
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        \langle U \rangle = -2 \langle K \rangle
+:laws:latex::
 """
 
 # TODO Prove using the virial's theorem.
diff --git a/symplyphysics/laws/gravity/radial_motion/potential_energy_of_planetary_motion.py b/symplyphysics/laws/gravity/radial_motion/potential_energy_of_planetary_motion.py
index fa6d27196..76620d754 100644
--- a/symplyphysics/laws/gravity/radial_motion/potential_energy_of_planetary_motion.py
+++ b/symplyphysics/laws/gravity/radial_motion/potential_energy_of_planetary_motion.py
@@ -6,46 +6,42 @@
 The potential energy is in turn the sum of the potential energy due to the gravitational interaction
 between the planet and the Sun, and the the energy of the tangential motion, which depends on the
 planet's angular momentum.
+
+**Links:**
+
+#. Sivukhin, D.V. (1979). *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 315.
 """
 
 from sympy import Eq
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
     symbols,
+    clone_as_symbol,
 )
 
-total_potential_energy = Symbol("total_potential_energy", units.energy)
-r"""
-The total potential energy of the planet.
-
-Symbol:
-    :code:`U_tot`
-
-Latex:
-    :math:`U_\text{tot}`
+total_potential_energy = clone_as_symbol(
+    symbols.potential_energy,
+    display_symbol="U_tot",
+    display_latex="U_\\text{tot}",
+)
 """
-
-gravitational_potential_energy = Symbol("gravitational_potential_energy", units.energy)
-r"""
-The potential energy of the planet due to the gravitational interaction of the planet and the star.
-
-Symbol:
-    :code:`U_gr`
-
-Latex:
-    :math:`U_\text{gr}`
+The total :symbols:`potential_energy` of the planet.
 """
 
-angular_momentum = Symbol("angular_momentum", units.length * units.momentum)
+gravitational_potential_energy = clone_as_symbol(
+    symbols.potential_energy,
+    display_symbol="U_gr",
+    display_latex="U_\\text{gr}",
+)
+"""
+The :symbols:`potential_energy` of the planet due to the gravitational interaction of the planet and the star.
 """
-The angular momentum of the planet.
 
-Symbol:
-    :code:`L`
+angular_momentum = symbols.angular_momentum
+"""
+The :symbols:`angular_momentum` of the planet.
 """
 
 planetary_mass = symbols.mass
@@ -53,24 +49,19 @@
 The :symbols:`mass` of the planet.
 """
 
-distance = Symbol("distance", units.length)
+distance = symbols.euclidean_distance
 """
-The distance between the star and the planet.
-
-Symbol:
-    :code:`r`
+The :symbols:`euclidean_distance` between the star and the planet.
 """
 
 law = Eq(
     total_potential_energy,
     gravitational_potential_energy + angular_momentum**2 / (2 * planetary_mass * distance**2),
 )
-r"""
-:code:`U_tot = U_gr + L^2 / (2 * m * r^2)`
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        U_\text{tot} = U_\text{gr} + \frac{L^2}{2 m r^2}
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/radial_motion/radial_kinetic_energy_plus_potential_energy_is_constant.py b/symplyphysics/laws/gravity/radial_motion/radial_kinetic_energy_plus_potential_energy_is_constant.py
index 4ba0d5a62..acadc5e18 100644
--- a/symplyphysics/laws/gravity/radial_motion/radial_kinetic_energy_plus_potential_energy_is_constant.py
+++ b/symplyphysics/laws/gravity/radial_motion/radial_kinetic_energy_plus_potential_energy_is_constant.py
@@ -8,16 +8,19 @@
 #. If :math:`E < 0`, the orbit is elliptical.
 #. If :math:`E = 0`, the orbit is parabolical.
 #. If :math:`E > 0`, the orbit is hyperbolical.
+
+**Links:**
+
+#. `Physics LibreTexts <https://phys.libretexts.org/Bookshelves/Classical_Mechanics/Classical_Mechanics_(Dourmashkin)/25%3A_Celestial_Mechanics/25.03%3A_Energy_and_Angular_Momentum_Constants_of_the_Motion>`__.
 """
 
 from sympy import Eq
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
     symbols,
+    clone_as_symbol,
 )
 
 planetary_mass = symbols.mass
@@ -25,39 +28,26 @@
 The :symbols:`mass` of the planet.
 """
 
-radial_speed = Symbol("radial_speed", units.velocity)
+radial_speed = clone_as_symbol(symbols.speed, subscript="r")
 """
-The projection of the velocity vector in the radial direction.
-
-Symbol:
-    :code:`v_r`
-
-Latex:
-    :math:`v_r`
+The projection of the velocity vector in the radial direction. See :symbols:`speed`.
 """
 
-potential_energy = Symbol("potential_energy", units.energy)
+potential_energy = symbols.potential_energy
 """
 The :doc:`potential energy <laws.gravity.radial_motion.potential_energy_of_planetary_motion>` of the planet.
-
-Symbol:
-    :code:`U`
+See :symbols:`potential_energy`.
 """
 
-total_mechanical_energy = Symbol("total_mechanical_energy", units.energy)
-"""Total mechanical energy of the planet, assumed to be constant.
-
-Symbol:
-    :code:`E`
+total_mechanical_energy = symbols.mechanical_energy
+"""Total :symbols:`mechanical_energy` of the planet, assumed to be constant.
 """
 
 law = Eq(planetary_mass * radial_speed**2 / 2 + potential_energy, total_mechanical_energy)
-r"""
-:code:`m * v_r**2 / 2 + U = E`
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        \frac{1}{2} m v_r^2 + U = E
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/radial_motion/semimajor_axis_via_kepler_constant_and_total_energy.py b/symplyphysics/laws/gravity/radial_motion/semimajor_axis_via_kepler_constant_and_total_energy.py
index 8f01421a7..fa9399173 100644
--- a/symplyphysics/laws/gravity/radial_motion/semimajor_axis_via_kepler_constant_and_total_energy.py
+++ b/symplyphysics/laws/gravity/radial_motion/semimajor_axis_via_kepler_constant_and_total_energy.py
@@ -3,62 +3,48 @@
 =====================================================
 
 The semi-major axis of an orbiting planet depends on the Kepler's constant of
-the star---planet system and the total energy of the planet per unit of its mass.
+the star—planet system and the total energy of the planet per unit of its mass.
 
 **Notes:**
 
 #. Works for both elliptical (:math:`\varepsilon < 0`) and hyperbolical
    (:math:`\varepsilon > 0`) orbits.
+
+**Links:**
+
+#. Sivukhin D.V. (1979), *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 317, (58.2).
 """
 
 from sympy import Eq, pi
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
+    symbols,
 )
 
-semimajor_axis = Symbol("semimajor_axis", units.length)
+semimajor_axis = symbols.semimajor_axis
 """
-The semi-major axis of the planet's orbit.
-
-Symbol:
-    a
+The :symbols:`semimajor_axis` of the planet's orbit.
 """
 
-kepler_constant = Symbol("kepler_constant", units.length**3 / units.time**2)
-r"""
-The Kepler's constant, whose value is determined by the :symbols:`mass`
+kepler_constant = symbols.kepler_constant
+"""
+The :symbols:`kepler_constant`, whose value is determined by the :symbols:`mass`
 of the orbited star.
-
-Symbol:
-    K
-
-Latex:
-    :math:`\mathfrak{K}`
 """
 
-specific_energy = Symbol("specific_energy", units.energy / units.mass)
-r"""
-The total energy of the planet per unit of its :symbols:`mass`.
-Can be negative or positive depending on the sign of the planet's energy.
-
-Symbol:
-    epsilon
-
-Latex:
-    :math:`\varepsilon`
+specific_energy = symbols.specific_energy
+"""
+:symbols:`specific_energy` of the planet. Can be negative or positive depending on the sign of
+the planet's energy.
 """
 
 law = Eq(semimajor_axis, 2 * pi**2 * kepler_constant / abs(specific_energy))
-r"""
-a = 2 * pi^2 * K / abs(epsilon)
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        a = \frac{2 \pi^2 \mathfrak{K}}{|\varepsilon|}
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/radial_motion/semiminor_axis_of_elliptical_orbit_via_orbit_parameters.py b/symplyphysics/laws/gravity/radial_motion/semiminor_axis_of_elliptical_orbit_via_orbit_parameters.py
index 8ff150026..8713f421f 100644
--- a/symplyphysics/laws/gravity/radial_motion/semiminor_axis_of_elliptical_orbit_via_orbit_parameters.py
+++ b/symplyphysics/laws/gravity/radial_motion/semiminor_axis_of_elliptical_orbit_via_orbit_parameters.py
@@ -4,45 +4,39 @@
 
 The minor semiaxis can be found as a function of the sector speed of the planet,
 the major semiaxis of its orbit, and the mass of body that attracts it, such as the Sun.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. Sivukhin, D.V. (1979). *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 318, (58.5)
 """
 
 from sympy import Eq, sqrt
-from sympy.physics.units import gravitational_constant
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
     clone_as_symbol,
     symbols,
+    quantities,
 )
 
-semiminor_axis = Symbol("semiminor_axis", units.length)
+semiminor_axis = symbols.semiminor_axis
 """
-The semi-minor axis of the planet's orbit.
-
-Symbol:
-    b
+The :symbols:`semiminor_axis` of the planet's orbit.
 """
 
-sector_speed = Symbol("sector_speed", units.area / units.time)
-r"""
-The sector speed of the planet, i.e. the area swept by the planet per unit time.
-
-Symbol:
-    sigma
-
-Latex:
-    :math:`\sigma`
+sector_speed = symbols.sector_speed
 """
-
-semimajor_axis = Symbol("semimajor_axis", units.length)
+The :symbols:`sector_speed` of the planet, i.e. the area swept by the planet per unit time.
 """
-The semi-major axis of the planet's orbit.
 
-Symbol:
-    a
+semimajor_axis = symbols.semimajor_axis
+"""
+The :symbols:`semimajor_axis` of the planet's orbit.
 """
 
 attracting_mass = clone_as_symbol(symbols.mass, display_symbol="M", display_latex="M")
@@ -52,13 +46,12 @@
 
 law = Eq(
     semiminor_axis,
-    2 * sector_speed * sqrt(semimajor_axis / (gravitational_constant * attracting_mass)),
+    2 * sector_speed * sqrt(semimajor_axis / (quantities.gravitational_constant * attracting_mass)),
 )
-r"""
-b = 2 * sigma * sqrt(a / (G * M))
+"""
+:laws:symbol::
 
-Latex:
-    :math:`b = 2 \sigma \sqrt{\frac{a}{G M}}`
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/radial_motion/total_energy_is_negative_average_kinetic_energy.py b/symplyphysics/laws/gravity/radial_motion/total_energy_is_negative_average_kinetic_energy.py
index 3318a517c..3c49e647e 100644
--- a/symplyphysics/laws/gravity/radial_motion/total_energy_is_negative_average_kinetic_energy.py
+++ b/symplyphysics/laws/gravity/radial_motion/total_energy_is_negative_average_kinetic_energy.py
@@ -7,43 +7,40 @@
 **Conditions:**
 
 #. Works for elliptical (:math:`E < 0`) orbits.
+
+**Links:**
+
+#. `Physics LibreTexts. Section "Energy in Circular Orbits" <https://phys.libretexts.org/Workbench/PH_245_Textbook_V2/13%3A_Gravitation/13.05%3A_Satellite_Orbits_and_Energy>`__.
 """
 
 from sympy import Eq
 from symplyphysics import (
-    units,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
+    symbols,
+    clone_as_symbol,
 )
 
-total_mechanical_energy = Symbol("total_mechanical_energy", units.energy)
+total_mechanical_energy = symbols.mechanical_energy
 """
-The total mechanical energy of the planet.
-
-Symbol:
-    :code:`E`
+The total :symbols:`mechanical_energy` of the planet.
 """
 
-average_kinetic_energy = Symbol("average_kinetic_energy", units.energy)
-r"""
-The kinetic energy of the planet averaged over time.
-
-Symbol:
-    :code:`avg(K)`
-
-Latex:
-    :math:`\langle K \rangle`
+average_kinetic_energy = clone_as_symbol(
+    symbols.kinetic_energy,
+    display_symbol="avg(K)",
+    display_latex="\\langle K \\rangle",
+)
+"""
+The :symbols:`kinetic_energy` of the planet averaged over :symbols:`time`.
 """
 
 law = Eq(total_mechanical_energy, -1 * average_kinetic_energy)
-r"""
-:code:`E = -1 * avg(K)`
+"""
+:laws:symbol::
 
-Latex:
-    .. math::
-        E = - \langle K \rangle
+:laws:latex::
 """
 
 
diff --git a/symplyphysics/laws/gravity/range_of_flight_of_a_body_thrown_at_an_angle_to_horizon.py b/symplyphysics/laws/gravity/range_of_flight_of_a_body_thrown_at_an_angle_to_horizon.py
index 791fe289a..6f2a655b0 100644
--- a/symplyphysics/laws/gravity/range_of_flight_of_a_body_thrown_at_an_angle_to_horizon.py
+++ b/symplyphysics/laws/gravity/range_of_flight_of_a_body_thrown_at_an_angle_to_horizon.py
@@ -1,69 +1,88 @@
-from sympy import (Eq, solve, sin, simplify)
-from sympy.physics.units import acceleration_due_to_gravity as earth_free_fall_acceleration
-from symplyphysics import (units, Quantity, Symbol, print_expression, validate_input,
-    validate_output, angle_type)
+"""
+Horizontal displacement of projectile
+=====================================
+
+Let's say we throw the body at an angle to the horizon with some initial velocity.
+Then the horizontal displacement, or **range**, of the body depends on the initial
+speed, the angle of the throw and the acceleration of free fall.
+
+**Conditions:**
+
+#. The acceleration due to gravity is constant throughout the movement
+   of the projectile.
+
+#. The elevation of the initial and final points of the projectile are
+   the same.
+
+**Links:**
+
+#. `Physics LibreTexts. Projectile Motion, Range (3.3.15) <https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion>`__.
+"""
+
+from sympy import Eq, solve, sin
+from symplyphysics import Quantity, validate_input, validate_output, symbols, quantities, clone_as_symbol
 from symplyphysics.core.expr_comparisons import expr_equals
 from symplyphysics.laws.kinematics import position_via_constant_speed_and_time as distance_law
 from symplyphysics.laws.geometry import planar_projection_is_cosine as projection_law
 from symplyphysics.laws.gravity import maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon as time_law
 
-# Description
-## Let's say we throw the body at an angle to the horizon with some initial velocity.
-## Then the range of the throw depends on the initial velocity, the angle of the throw and
-## the acceleration of free fall.
+horizontal_displacement = symbols.euclidean_distance
+"""
+Horizontal displacement, or **range**, of the projectile. See :symbols:`euclidean_distance`.
+"""
 
-## Law is: L = v0^2 * sin(2 * a) / g, where
-## L - throw range,
-## v0 - initial velocity,
-## a - angle of throw,
-## g - acceleration of free fall.
+initial_speed = clone_as_symbol(symbols.speed, subscript="0")
+"""
+Initial :symbols:`speed` of the projectile.
+"""
 
-throw_range = Symbol("throw_range", units.length)
+angle = symbols.angle
+"""
+:symbols:`angle` between the initial velocity and the horizon.
+"""
 
-initial_velocity = Symbol("initial_velocity", units.velocity)
-angle = Symbol("angle", angle_type)
+law = Eq(horizontal_displacement, initial_speed**2 * sin(2 * angle) / quantities.acceleration_due_to_gravity)
+"""
+:laws:symbol::
 
-law = Eq(throw_range, initial_velocity**2 * sin(2 * angle) / earth_free_fall_acceleration)
+:laws:latex::
+"""
 
 # This law might be derived via "distance_from_constant_velocity" law, "planar_projection_is_cosine" law
 # and "maximum_movement_time_of_a_body_thrown_at_an_angle_to_horizon" law.
 
-projection_law_applied = projection_law.law.subs({
-    projection_law.vector_length: initial_velocity,
+_projection_law_applied = projection_law.law.subs({
+    projection_law.vector_length: initial_speed,
     projection_law.vector_angle: angle,
 })
-horizontal_projection_derived = solve(projection_law_applied, projection_law.projection,
+_horizontal_projection_derived = solve(_projection_law_applied, projection_law.projection,
     dict=True)[0][projection_law.projection]
 
-time_law_applied = time_law.law.subs({
-    time_law.initial_velocity: initial_velocity,
+_time_law_applied = time_law.law.subs({
+    time_law.initial_speed: initial_speed,
     time_law.angle: angle,
 })
-time_derived = solve(time_law_applied, time_law.movement_time, dict=True)[0][time_law.movement_time]
+_time_derived = solve(_time_law_applied, time_law.time, dict=True)[0][time_law.time]
 
-range_law_applied = distance_law.law.subs({
+_range_law_applied = distance_law.law.subs({
     distance_law.initial_position: 0,
-    distance_law.speed: horizontal_projection_derived,
-    distance_law.time: time_derived,
+    distance_law.speed: _horizontal_projection_derived,
+    distance_law.time: _time_derived,
 })
 
 # Check if derived range is same as declared.
-assert expr_equals(simplify(range_law_applied.rhs), law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_range_law_applied.rhs, law.rhs)
 
 
-@validate_input(initial_velocity_=initial_velocity, angle_=angle)
-@validate_output(throw_range)
+@validate_input(initial_velocity_=initial_speed, angle_=angle)
+@validate_output(horizontal_displacement)
 def calculate_range(
     initial_velocity_: Quantity,
     angle_: float | Quantity,
 ) -> Quantity:
-    result_expr = solve(law, throw_range, dict=True)[0][throw_range]
+    result_expr = solve(law, horizontal_displacement, dict=True)[0][horizontal_displacement]
     result_expr = result_expr.subs({
-        initial_velocity: initial_velocity_,
+        initial_speed: initial_velocity_,
         angle: angle_,
     })
     return Quantity(result_expr)
diff --git a/symplyphysics/laws/gravity/second_escape_velocity.py b/symplyphysics/laws/gravity/second_escape_velocity.py
index f6a1fe339..1d974eb04 100644
--- a/symplyphysics/laws/gravity/second_escape_velocity.py
+++ b/symplyphysics/laws/gravity/second_escape_velocity.py
@@ -1,49 +1,71 @@
+"""
+Second escape velocity
+======================
+
+The second cosmic velocity is the lowest velocity that must be given to an object starting from the surface of
+a celestial body, the mass of which is negligible compared to the mass of a celestial body, in order to overcome
+the gravitational attraction of this celestial body and leave a closed orbit around it.
+
+**Notation:**
+
+#. :quantity_notation:`gravitational_constant`.
+
+**Links:**
+
+#. `Physics LibreTexts. Conservation of Energy, Escape velocity (13.4.10) <https://phys.libretexts.org/Bookshelves/University_Physics/University_Physics_(OpenStax)/Book%3A_University_Physics_I_-_Mechanics_Sound_Oscillations_and_Waves_(OpenStax)/13%3A_Gravitation/13.04%3A_Gravitational_Potential_Energy_and_Total_Energy>`__.
+"""
+
 from sympy import Eq, solve, sqrt
-from sympy.physics.units import gravitational_constant
-from symplyphysics import (symbols, units, Quantity, Symbol, print_expression, validate_input,
-    validate_output)
+from symplyphysics import symbols, Quantity, validate_input, validate_output, quantities
 from symplyphysics.core.expr_comparisons import expr_equals
-
 from symplyphysics.laws.gravity import gravitational_potential_energy as potential_energy_law
 from symplyphysics.laws.dynamics import kinetic_energy_from_mass_and_speed as kinetic_energy_law
 from symplyphysics.laws.conservation import mechanical_energy_after_equals_to_mechanical_energy_before as conservation_law
 
-# The second cosmic velocity is the lowest velocity that must be given to an object starting from the surface of
-# a celestial body, the mass of which is negligible compared to the mass of a celestial body, in order to overcome
-# the gravitational attraction of this celestial body and leave a closed orbit around it.
+speed = symbols.speed
+"""
+Second escape :symbols:`speed`.
+"""
 
-# Law: V = √(2 * G * M / (R + h))
-# Where:
-# V - initial velocity
-# M - planet_mass of planet
-# G - gravitational constant
-# h - height above the planet surface
-# R - radius of planet
-
-velocity = Symbol("initial_velocity", units.velocity)
 planet_mass = symbols.mass
-radius = Symbol("radius", units.length)
-height = Symbol("height", units.length)
+"""
+Planetary :symbols:`mass`.
+"""
+
+planet_radius = symbols.radius
+"""
+Planetary :symbols:`radius`.
+"""
+
+height = symbols.height
+"""
+Elevation (see :symbols:`height`) of the body above the planet's surface.
+"""
 
-law = Eq(velocity, sqrt(2 * gravitational_constant * planet_mass / (radius + height)))
+law = Eq(speed, sqrt(2 * quantities.gravitational_constant * planet_mass / (planet_radius + height)))
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # This law might be derived via "gravitational_potential_energy" law, "kinetic_energy_from_mass_and_velocity" law,
 # "mechanical_energy_after_equals_to_mechanical_energy_before" law.
 
 # The radius of the orbit consists of the radius of the planet and the height above its surface.
-potential_energy_law_applied = potential_energy_law.law.subs({
+_potential_energy_law_applied = potential_energy_law.law.subs({
     potential_energy_law.first_mass: planet_mass,
-    potential_energy_law.distance_between_mass_centers: radius + height,
+    potential_energy_law.distance_between_mass_centers: planet_radius + height,
 })
-potential_energy_derived = solve(potential_energy_law_applied,
+_potential_energy_derived = solve(_potential_energy_law_applied,
     potential_energy_law.gravitational_potential_energy,
     dict=True)[0][potential_energy_law.gravitational_potential_energy]
 
-kinetic_energy_law_applied = kinetic_energy_law.law.subs({
+_kinetic_energy_law_applied = kinetic_energy_law.law.subs({
     kinetic_energy_law.mass: potential_energy_law.second_mass,
-    kinetic_energy_law.speed: velocity,
+    kinetic_energy_law.speed: speed,
 })
-kinetic_energy_derived = solve(kinetic_energy_law_applied,
+_kinetic_energy_derived = solve(_kinetic_energy_law_applied,
     kinetic_energy_law.kinetic_energy,
     dict=True)[0][kinetic_energy_law.kinetic_energy]
 
@@ -51,30 +73,26 @@
 # look like: kinetic energy minus potential. But in the law "gravitational_potential_energy", the potential energy has
 # already been introduced with a minus. And in the law "mechanical_energy_after_equals_to_mechanical_energy_before",
 # the conservation law looks like: E1=E2. Therefore, a minus sign is placed before the potential energy.
-conservation_law_applied = conservation_law.law.subs({
-    conservation_law.mechanical_energy(conservation_law.time_after): kinetic_energy_derived,
-    conservation_law.mechanical_energy(conservation_law.time_before): -1 * potential_energy_derived,
+_conservation_law_applied = conservation_law.law.subs({
+    conservation_law.mechanical_energy(conservation_law.time_after): _kinetic_energy_derived,
+    conservation_law.mechanical_energy(conservation_law.time_before): -1 * _potential_energy_derived,
 })
 # The proof considers the fall of a body from infinity to the earth. But at the same time, in most cases, I consider
 # this speed as directed from the planet from which the body starts, which corresponds to the plus sign. Therefore,
 # the first solution is ignored.
-velocity_derived = solve(conservation_law_applied, velocity, dict=True)[1][velocity]
+_velocity_derived = solve(_conservation_law_applied, speed, dict=True)[1][speed]
 
 # Check if derived velocity is same as declared.
-assert expr_equals(velocity_derived, law.rhs)
-
-
-def print_law() -> str:
-    return print_expression(law)
+assert expr_equals(_velocity_derived, law.rhs)
 
 
-@validate_input(planet_mass_=planet_mass, radius_=radius, height_=height)
-@validate_output(velocity)
+@validate_input(planet_mass_=planet_mass, radius_=planet_radius, height_=height)
+@validate_output(speed)
 def calculate_velocity(planet_mass_: Quantity, radius_: Quantity, height_: Quantity) -> Quantity:
-    result_velocity_expr = solve(law, velocity, dict=True)[0][velocity]
+    result_velocity_expr = solve(law, speed, dict=True)[0][speed]
     result_expr = result_velocity_expr.subs({
         planet_mass: planet_mass_,
-        radius: radius_,
+        planet_radius: radius_,
         height: height_
     })
     return Quantity(result_expr)
diff --git a/symplyphysics/laws/gravity/southerly_deviation_from_plumbline_of_falling_bodies.py b/symplyphysics/laws/gravity/southerly_deviation_from_plumbline_of_falling_bodies.py
index 415e9f742..c6cf6f093 100644
--- a/symplyphysics/laws/gravity/southerly_deviation_from_plumbline_of_falling_bodies.py
+++ b/symplyphysics/laws/gravity/southerly_deviation_from_plumbline_of_falling_bodies.py
@@ -1,42 +1,79 @@
+"""
+Southerly deviation from plumbline of falling bodies
+====================================================
+
+Suppose a body is falling freely in Earth's gravity field with its initial velocity being zero.
+Then the effect of the Coriolis force on the falling body can be found in the fact that it deflects
+from plumbline in the easterly and southerly (equatorial) directions.
+
+**Notes:**
+
+#. The southerly deviation is extremely small and almost unobservable due to the :math:`t / T` factor.
+
+#. Also see :ref:`Easterly deviation from plumbline of falling bodies`.
+
+**Conditions:**
+
+#. The vector of free fall acceleration is considered constant.
+
+**Links:**
+
+#. Sivukhin D.V. (1979), *Obshchiy kurs fiziki* [General course of Physics], vol. 1, p. 355, (67.11).
+"""
+
 from sympy import Eq, sin, pi
 from symplyphysics import (
-    units,
-    angle_type,
     Quantity,
-    Symbol,
     validate_input,
     validate_output,
+    symbols,
+    clone_as_symbol,
 )
 from symplyphysics.core.symbols.quantities import scale_factor
 
-# Description
-## Suppose a body is falling freely in Earth's gravity field with its initial velocity being zero.
-## Then the effect of the Coriolis force on the falling body can be found in the fact that it deflects
-## from plumbline in the easterly and southerly (equatorial) directions.
+southerly_deviation_from_plumbline = clone_as_symbol(
+    symbols.euclidean_distance,
+    display_symbol="s_south",
+    display_latex="s_\\text{south}",
+)
+"""
+Southerly (equatorial) deviation of falling body from plumbline due to Earth's rotation.
+See :symbols:`euclidean_distance`.
+"""
 
-# Conditions
-## - The vector of free fall acceleration `g` is considered constant.
+fall_time = symbols.time
+"""
+:symbols:`time` elapsed during the body's fall.
+"""
 
-# Notes
-## - The southerly deviation is extremely small and almost unobservable due to the `t / T` factor.
+rotation_period = symbols.period
+"""
+:symbols:`period` of Earth's rotation.
+"""
 
-# Law: s_south = pi * (t / T) * s_east * sin(theta)
-## s_sourth - southerly (equatorial) deviation of falling body from plumbline due to Earth's rotation
-## t - time of the body's fall
-## T - period of Earth's rotation
-## s_east - easterly deviation of falling body from plumbline
-## theta - latitude of the place the body is located in
+easterly_deviation_from_plumbline = clone_as_symbol(
+    symbols.euclidean_distance,
+    display_symbol="s_east",
+    display_latex="s_\\text{east}",
+)
+"""
+Easterly deviation of falling body from plumbline. See :symbols:`euclidean_distance`.
+"""
 
-southerly_deviation_from_plumbline = Symbol("southerly_deviation_from_plumbline", units.length)
-fall_time = Symbol("fall_time", units.time)
-rotation_period = Symbol("rotation_period", units.time)
-easterly_deviation_from_plumbline = Symbol("easterly_deviation_from_plumbline", units.length)
-latitude = Symbol("latitude", angle_type)
+latitude = symbols.latitude
+"""
+:symbols:`latitude` of the place the body is located in.
+"""
 
 law = Eq(
-    easterly_deviation_from_plumbline,
+    southerly_deviation_from_plumbline,
     pi * (fall_time / rotation_period) * easterly_deviation_from_plumbline * sin(latitude),
 )
+"""
+:laws:symbol::
+
+:laws:latex::
+"""
 
 # TODO: derive from the solution of the relative motion equation in Earth's gravitational field.
 
diff --git a/symplyphysics/laws/gravity/relative_acceleration_from_force.py b/symplyphysics/laws/gravity/vector/relative_acceleration_from_force.py
similarity index 100%
rename from symplyphysics/laws/gravity/relative_acceleration_from_force.py
rename to symplyphysics/laws/gravity/vector/relative_acceleration_from_force.py
diff --git a/symplyphysics/quantities/__init__.py b/symplyphysics/quantities/__init__.py
index 5bd85597d..28723edfd 100644
--- a/symplyphysics/quantities/__init__.py
+++ b/symplyphysics/quantities/__init__.py
@@ -166,6 +166,16 @@
 #. `Wien's displacement law <https://en.wikipedia.org/wiki/Wien%27s_displacement_law>`__.
 """
 
+gravitational_constant = Quantity(units.gravitational_constant, display_symbol="G")
+"""
+The **gravitational constant** is a physical constant used in calculating the gravitational
+attraction between two objects.
+
+**Links:**
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Gravitational_constant>`__.
+"""
+
 
 __all__ = [
     "standard_conditions_temperature",
diff --git a/symplyphysics/symbols/__init__.py b/symplyphysics/symbols/__init__.py
index 50e7b6e2c..9901f0677 100644
--- a/symplyphysics/symbols/__init__.py
+++ b/symplyphysics/symbols/__init__.py
@@ -38,6 +38,7 @@
     "any_quantity",
     "fractional_change",
     "exponential_decay_constant",
+    "specific_energy",
     # chemistry,
     "mass_fraction",
     "amount_of_substance",
@@ -103,6 +104,12 @@
     "degrees_of_freedom",
     "angular_momentum",
     "arc_length",
+    "semimajor_axis",
+    "semiminor_axis",
+    "latitude",
+    "longitude",
+    "sector_speed",
+    "kepler_constant",
     # electrodynamics
     "admittance",
     "electrical_conductance",
diff --git a/symplyphysics/symbols/basic.py b/symplyphysics/symbols/basic.py
index 0cab1ab64..b5fb60a0e 100644
--- a/symplyphysics/symbols/basic.py
+++ b/symplyphysics/symbols/basic.py
@@ -61,6 +61,11 @@
 the performance of work and in the form of heat and light.
 """
 
+specific_energy = SymbolNew("epsilon", units.energy / units.mass, display_latex="\\varepsilon")
+"""
+**Specific energy** is defined as :symbols:`energy` per unit :symbols:`mass`.
+"""
+
 power = SymbolNew("P", units.power)
 """
 **Power** is the amount of energy transferred or converted per unit time.
diff --git a/symplyphysics/symbols/classical_mechanics.py b/symplyphysics/symbols/classical_mechanics.py
index d69df26c2..191822c93 100644
--- a/symplyphysics/symbols/classical_mechanics.py
+++ b/symplyphysics/symbols/classical_mechanics.py
@@ -93,6 +93,16 @@
 the sphere.
 """
 
+semimajor_axis = SymbolNew("a", units.length)
+"""
+**Semi-major axis**, or **major semiaxis**, is the longest semidiameter of an ellipse.
+"""
+
+semiminor_axis = SymbolNew("b", units.length)
+"""
+**Semi-minor axis**, or **minor semiaxis**, is the smallest semidiameter of an ellipse.
+"""
+
 thickness = SymbolNew("h", units.length)
 """
 **Thickness** is a measure of a size of an object, usually the separation between two layers, or
@@ -362,3 +372,35 @@
 **Angular momentum**, sometimes called **rotational momentum**, is the rotational analog of linear
 :symbols:`momentum`.
 """
+
+latitude = SymbolNew("phi", angle_type, display_latex="\\phi")
+"""
+**Latitude** is a coordinate that specifies the north-south position of a point on the surface of
+the Earth or another celestial body. Its value ranges from :math:`-90^\circ` at the south pole
+to :math:`90^\circ` at the north pole, with :math:`0^\circ` at the Equator.
+"""
+
+longitude = SymbolNew("lambda", angle_type, display_latex="\\lambda")
+"""
+**Longitude** is a coordinate that specifies the east-west position of a point on the
+surface of the Earth, or another celestial body. The prime meridian defines :math:`0^\circ`
+longitude; positive longitudes are east of the prime meridian, and the negative ones are west.
+"""
+
+sector_speed = SymbolNew("sigma", units.area / units.time, display_latex="\\sigma")
+"""
+**Areal speed**, also called **sector speed** or **sectorial speed**, is a quantity that indicates
+the rate of change at which :symbols:`area` is swept out by a particle as it moves along a curve.
+The calculation of the area swept can be found in the link below.
+
+**Links:**
+
+#. `Wikipedia <https://en.wikipedia.org/wiki/Areal_velocity#/media/File:ArealVelocity_with_curved_area.svg>`__.
+"""
+
+kepler_constant = SymbolNew("K", units.length**3 / units.time**2, display_latex="\\mathfrak{K}")
+"""
+The **Kepler's constant** is the constant of proportionality in Kepler's third law of planetary
+motion, namely the ratio between the square of the period of the planet to the semi-major axis
+of the planet's orbit. It is constant for all objects orbiting around the same object.
+"""
diff --git a/test/gravity/relative_acceleration_from_force_gravity_test.py b/test/gravity/vector/relative_acceleration_from_force_gravity_test.py
similarity index 98%
rename from test/gravity/relative_acceleration_from_force_gravity_test.py
rename to test/gravity/vector/relative_acceleration_from_force_gravity_test.py
index 9567866db..f36168e70 100644
--- a/test/gravity/relative_acceleration_from_force_gravity_test.py
+++ b/test/gravity/vector/relative_acceleration_from_force_gravity_test.py
@@ -7,7 +7,7 @@
     Quantity,
     QuantityVector,
 )
-from symplyphysics.laws.gravity import relative_acceleration_from_force as law
+from symplyphysics.laws.gravity.vector import relative_acceleration_from_force as law
 from symplyphysics.core.approx import assert_equal_vectors
 
 Args = namedtuple("Args", "g ac f m ar")