diff --git a/boule/__init__.py b/boule/__init__.py
index 00640f11..e83adcfd 100644
--- a/boule/__init__.py
+++ b/boule/__init__.py
@@ -2,6 +2,7 @@
 # Import functions/classes to make the public API
 from . import version
 from .ellipsoid import Ellipsoid
+from .sphere import Sphere
 from .realizations import WGS84, GRS80, MARS
 
 
diff --git a/boule/ellipsoid.py b/boule/ellipsoid.py
index 441d9d6d..8ed4217d 100644
--- a/boule/ellipsoid.py
+++ b/boule/ellipsoid.py
@@ -1,6 +1,7 @@
 """
 Module for defining and setting the reference ellipsoid.
 """
+from warnings import warn
 import attr
 import numpy as np
 
@@ -41,30 +42,30 @@ class Ellipsoid:
     Examples
     --------
 
-    We can define a reference unit sphere by using 0 as the flattening:
+    We can define an ellipsoid with flattening equal to 0.5 and unit semimajor axis:
 
-    >>> sphere = Ellipsoid(
-    ...     name="sphere",
-    ...     long_name="Unit sphere",
+    >>> ellipsoid = Ellipsoid(
+    ...     name="oblate-ellipsoid",
+    ...     long_name="Oblate Ellipsoid",
     ...     semimajor_axis=1,
-    ...     flattening=0,
+    ...     flattening=0.5,
     ...     geocentric_grav_const=1,
     ...     angular_velocity=0
     ... )
-    >>> print(sphere) # doctest: +ELLIPSIS
-    Ellipsoid(name='sphere', ...)
-    >>> print(sphere.long_name)
-    Unit sphere
-    >>> print("{:.2f}".format(sphere.semiminor_axis))
-    1.00
-    >>> print("{:.2f}".format(sphere.mean_radius))
-    1.00
-    >>> print("{:.2f}".format(sphere.linear_eccentricity))
-    0.00
-    >>> print("{:.2f}".format(sphere.first_eccentricity))
-    0.00
-    >>> print("{:.2f}".format(sphere.second_eccentricity))
-    0.00
+    >>> print(ellipsoid) # doctest: +ELLIPSIS
+    Ellipsoid(name='oblate-ellipsoid', ...)
+    >>> print(ellipsoid.long_name)
+    Oblate Ellipsoid
+    >>> print("{:.2f}".format(ellipsoid.semiminor_axis))
+    0.50
+    >>> print("{:.2f}".format(ellipsoid.mean_radius))
+    0.83
+    >>> print("{:.2f}".format(ellipsoid.linear_eccentricity))
+    0.87
+    >>> print("{:.2f}".format(ellipsoid.first_eccentricity))
+    0.87
+    >>> print("{:.2f}".format(ellipsoid.second_eccentricity))
+    1.73
 
     """
 
@@ -76,6 +77,19 @@ class Ellipsoid:
     long_name = attr.ib(default=None)
     reference = attr.ib(default=None)
 
+    @flattening.validator
+    def check_flattening(
+        self, flattening, value
+    ):  # pylint: disable=no-self-use,unused-argument
+        """
+        Check if flattening is not equal (or almost) zero
+        """
+        warn_msg = "Use boule.Sphere for representing ellipsoids with zero flattening."
+        if value == 0:
+            warn("Flattening equal to zero. " + warn_msg)
+        if value < 1e-7:
+            warn("Flattening '{}' too close to zero. ".format(value) + warn_msg)
+
     @property
     def semiminor_axis(self):
         "The small (polar) axis of the ellipsoid [meters]"
diff --git a/boule/sphere.py b/boule/sphere.py
new file mode 100644
index 00000000..498ec3ba
--- /dev/null
+++ b/boule/sphere.py
@@ -0,0 +1,161 @@
+"""
+Module for defining and setting the reference sphere.
+"""
+import attr
+import numpy as np
+
+from . import Ellipsoid
+
+
+# Don't let ellipsoid parameters be changed to avoid messing up calculations
+# accidentally.
+@attr.s(frozen=True)
+class Sphere(Ellipsoid):
+    """
+    Reference spherical ellipsoid
+
+    Represents ellipsoids with zero flattening (spheres). Inherits methods and
+    properties of the :class:`boule.Ellipsoid`, guaranteeing no singularities
+    due to zero flattening (and thus zero eccentricity).
+
+    All parameters are in SI units.
+
+    Parameters
+    ----------
+    name : str
+        A short name for the ellipsoid, for example ``'MOON'``.
+    radius : float
+        The radius of the sphere [meters].
+    geocentric_grav_const : float
+        The geocentric gravitational constant (GM) [m^3 s^-2].
+    angular_velocity : float
+        The angular velocity of the rotating ellipsoid (omega) [rad s^-1].
+    long_name : str or None
+        A long name for the ellipsoid, for example ``"Moon Reference System"``
+        (optional).
+    reference : str or None
+        Citation for the ellipsoid parameter values (optional).
+
+    Examples
+    --------
+
+    We can define a unit sphere:
+
+    >>> sphere = Sphere(
+    ...     name="sphere",
+    ...     radius=1,
+    ...     geocentric_grav_const=1,
+    ...     angular_velocity=0,
+    ...     long_name="Spherical Ellipsoid",
+    ... )
+    >>> print(sphere) # doctest: +ELLIPSIS
+    Sphere(name='sphere', ...)
+    >>> print(sphere.long_name)
+    Spherical Ellipsoid
+    >>> print("{:.2f}".format(sphere.semiminor_axis))
+    1.00
+    >>> print("{:.2f}".format(sphere.mean_radius))
+    1.00
+    >>> print("{:.2f}".format(sphere.flattening))
+    0.00
+    >>> print("{:.2f}".format(sphere.linear_eccentricity))
+    0.00
+    >>> print("{:.2f}".format(sphere.first_eccentricity))
+    0.00
+    >>> print("{:.2f}".format(sphere.second_eccentricity))
+    0.00
+    >>> print(sphere.normal_gravity(latitude=0, height=1))
+    25000.0
+    >>> print(sphere.normal_gravity(latitude=90, height=1))
+    25000.0
+    """
+
+    name = attr.ib()
+    radius = attr.ib()
+    geocentric_grav_const = attr.ib()
+    angular_velocity = attr.ib()
+    long_name = attr.ib(default=None)
+    reference = attr.ib(default=None)
+    # semimajor_axis and flattening shouldn't be defined on initialization:
+    #   - semimajor_axis will be equal to radius
+    #   - flattening will be equal to zero
+    semimajor_axis = attr.ib(init=False)
+    flattening = attr.ib(init=False, default=0)
+
+    @semimajor_axis.default
+    def _set_semimajor_axis(self):
+        "The semimajor axis should be the radius"
+        return self.radius
+
+    def normal_gravity(self, latitude, height):
+        """
+        Calculate normal gravity at any latitude and height
+
+        Computes the magnitude of the gradient of the gravity potential
+        (gravitational + centrifugal) generated by the sphere at the given
+        latitude and height.
+
+        Parameters
+        ----------
+        latitude : float or array
+            The latitude where the normal gravity will be computed (in degrees).
+        height : float or array
+            The height (above the sphere) of computation point (in meters).
+
+        Returns
+        -------
+        gamma : float or array
+            The normal gravity in mGal.
+
+        References
+        ----------
+        [Heiskanen-Moritz]_
+        """
+        return self._gravity_sphere(height) + self._centrifugal_force(latitude, height)
+
+    def _gravity_sphere(self, height):
+        """
+        Calculate the gravity generated by a solid sphere (mGal)
+        """
+        return 1e5 * self.geocentric_grav_const / (self.radius + height) ** 2
+
+    def _centrifugal_force(self, latitude, height):
+        """
+        Calculate the centrifugal force due to the rotation of the sphere (mGal)
+        """
+        return 1e5 * (
+            (-1)
+            * self.angular_velocity ** 2
+            * (self.radius + height)
+            * np.cos(np.radians(latitude))
+        )
+
+    @property
+    def gravity_equator(self):
+        """
+        The norm of the gravity vector at the equator on the sphere [m/s^2]
+
+        Overrides the inherited method from :class:`boule.Ellipsoid` to avoid
+        singularities due to zero flattening.
+        """
+        return self._gravity_on_surface
+
+    @property
+    def gravity_pole(self):
+        """
+        The norm of the gravity vector at the poles on the sphere [m/s^2]
+
+        Overrides the inherited method from :class:`boule.Ellipsoid` to avoid
+        singularities due to zero flattening.
+        """
+        return self._gravity_on_surface
+
+    @property
+    def _gravity_on_surface(self):
+        """
+        Compute norm of the gravity vector on the surface of the sphere [m/s^2]
+
+        Due to rotational symmetry, the norm of the gravity vector is the same
+        on every point of the surface.
+        """
+        return self.geocentric_grav_const / self.radius ** 2
diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py
index fe07cd4b..9931889c 100644
--- a/boule/tests/test_ellipsoid.py
+++ b/boule/tests/test_ellipsoid.py
@@ -2,6 +2,7 @@
 """
 Test the base Ellipsoid class.
 """
+import warnings
 import pytest
 import numpy as np
 import numpy.testing as npt
@@ -12,47 +13,30 @@
 ELLIPSOID_NAMES = [e.name for e in ELLIPSOIDS]
 
 
-@pytest.fixture
-def sphere():
-    "A spherical ellipsoid"
-    ellipsoid = Ellipsoid(
-        name="unit_sphere",
-        semimajor_axis=1.0,
-        flattening=0,
-        geocentric_grav_const=0,
-        angular_velocity=0,
-    )
-    return ellipsoid
-
-
-def test_geodetic_to_spherical_with_spherical_ellipsoid(sphere):
-    "Test geodetic to geocentric conversion if ellipsoid is a sphere."
-    rtol = 1e-10
-    size = 5
-    longitude = np.linspace(0, 180, size)
-    latitude = np.linspace(-90, 90, size)
-    height = np.linspace(-0.2, 0.2, size)
-    sph_longitude, sph_latitude, radius = sphere.geodetic_to_spherical(
-        longitude, latitude, height
-    )
-    npt.assert_allclose(sph_longitude, longitude, rtol=rtol)
-    npt.assert_allclose(sph_latitude, latitude, rtol=rtol)
-    npt.assert_allclose(radius, sphere.mean_radius + height, rtol=rtol)
-
-
-def test_spherical_to_geodetic_with_spherical_ellipsoid(sphere):
-    "Test spherical to geodetic conversion if ellipsoid is a sphere."
-    rtol = 1e-10
-    size = 5
-    spherical_longitude = np.linspace(0, 180, size)
-    spherical_latitude = np.linspace(-90, 90, size)
-    radius = np.linspace(0.8, 1.2, size)
-    longitude, latitude, height = sphere.spherical_to_geodetic(
-        spherical_longitude, spherical_latitude, radius
-    )
-    npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
-    npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
-    npt.assert_allclose(radius, sphere.mean_radius + height, rtol=rtol)
+def test_ellipsoid_zero_flattening():
+    """
+    Check if error is raised after passing zero flattening
+    """
+    # Test with zero flattening
+    with warnings.catch_warnings(record=True) as warn:
+        Ellipsoid(
+            name="zero-flattening",
+            semimajor_axis=1,
+            flattening=0,
+            geocentric_grav_const=1,
+            angular_velocity=0,
+        )
+        assert len(warn) >= 1
+    # Test with almost zero flattening
+    with warnings.catch_warnings(record=True) as warn:
+        Ellipsoid(
+            name="almost-zero-flattening",
+            semimajor_axis=1,
+            flattening=1e-8,
+            geocentric_grav_const=1,
+            angular_velocity=0,
+        )
+        assert len(warn) >= 1
 
 
 @pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
diff --git a/boule/tests/test_sphere.py b/boule/tests/test_sphere.py
new file mode 100644
index 00000000..f4daedd5
--- /dev/null
+++ b/boule/tests/test_sphere.py
@@ -0,0 +1,157 @@
+# pylint: disable=redefined-outer-name
+"""
+Test the base Sphere class.
+"""
+import pytest
+import numpy as np
+import numpy.testing as npt
+
+from .. import Sphere
+
+
+@pytest.fixture
+def sphere():
+    "A spherical ellipsoid"
+    ellipsoid = Sphere(
+        name="unit_sphere", radius=1.0, geocentric_grav_const=2.0, angular_velocity=1.3,
+    )
+    return ellipsoid
+
+
+def test_sphere_flattening(sphere):
+    """
+    Check if flattening property is equal to zero
+    """
+    assert sphere.flattening == 0
+
+
+def test_sphere_semimajor_axis(sphere):
+    """
+    Check if semimajor_axis is equal to the radius
+    """
+    npt.assert_allclose(sphere.semimajor_axis, sphere.radius)
+
+
+def test_geodetic_to_spherical(sphere):
+    "Test geodetic to geocentric conversion on spherical ellipsoid."
+    rtol = 1e-10
+    size = 5
+    longitude = np.linspace(0, 180, size)
+    latitude = np.linspace(-90, 90, size)
+    height = np.linspace(-0.2, 0.2, size)
+    sph_longitude, sph_latitude, radius = sphere.geodetic_to_spherical(
+        longitude, latitude, height
+    )
+    npt.assert_allclose(sph_longitude, longitude, rtol=rtol)
+    npt.assert_allclose(sph_latitude, latitude, rtol=rtol)
+    npt.assert_allclose(radius, sphere.radius + height, rtol=rtol)
+
+
+def test_spherical_to_geodetic(sphere):
+    "Test spherical to geodetic conversion on spherical ellipsoid."
+    rtol = 1e-10
+    size = 5
+    spherical_longitude = np.linspace(0, 180, size)
+    spherical_latitude = np.linspace(-90, 90, size)
+    radius = np.linspace(0.8, 1.2, size)
+    longitude, latitude, height = sphere.spherical_to_geodetic(
+        spherical_longitude, spherical_latitude, radius
+    )
+    npt.assert_allclose(spherical_longitude, longitude, rtol=rtol)
+    npt.assert_allclose(spherical_latitude, latitude, rtol=rtol)
+    npt.assert_allclose(radius, sphere.radius + height, rtol=rtol)
+
+
+def test_ellipsoidal_properties(sphere):
+    """
+    Check inherited properties from Ellipsoid
+    """
+    npt.assert_allclose(sphere.semiminor_axis, sphere.radius)
+    npt.assert_allclose(sphere.linear_eccentricity, 0)
+    npt.assert_allclose(sphere.first_eccentricity, 0)
+    npt.assert_allclose(sphere.second_eccentricity, 0)
+    npt.assert_allclose(sphere.mean_radius, sphere.radius)
+    npt.assert_allclose(
+        sphere.emm,
+        sphere.angular_velocity ** 2
+        * sphere.radius ** 3
+        / sphere.geocentric_grav_const,
+    )
+
+
+def test_prime_vertical_radius(sphere):
+    """
+    Check prime vertical radius
+    """
+    latitudes = np.linspace(-90, 90, 18)
+    npt.assert_allclose(
+        sphere.prime_vertical_radius(np.sin(np.radians(latitudes))), sphere.radius
+    )
+
+
+def test_geocentric_radius(sphere):
+    """
+    Check geocentric radius
+    """
+    latitudes = np.linspace(-90, 90, 18)
+    for geodetic in (True, False):
+        npt.assert_allclose(
+            sphere.geocentric_radius(latitudes, geodetic=geodetic), sphere.radius
+        )
+
+
+def test_gravity_on_equator_and_poles(sphere):
+    """
+    Check gravity on equator and poles
+    """
+    expected_value = sphere.geocentric_grav_const / sphere.radius ** 2
+    npt.assert_allclose(sphere.gravity_pole, expected_value)
+    npt.assert_allclose(sphere.gravity_equator, expected_value)
+
+
+def test_normal_gravity_no_rotation():
+    """
+    Check normal gravity without rotation
+    """
+    gm_constant = 3
+    radius = 1
+    sphere = Sphere(
+        name="sphere",
+        radius=radius,
+        geocentric_grav_const=gm_constant,
+        angular_velocity=0,
+    )
+    # Create a set of points a different latitudes and same height
+    for height in [1, 2, 3, 4]:
+        latitudes = np.linspace(-90, 90, 19)
+        heights = height * np.ones_like(latitudes)
+        # Check if normal gravity is equal on every point (rotational symmetry)
+        expected_gravity = 1e5 * gm_constant / (radius + height) ** 2
+        npt.assert_allclose(expected_gravity, sphere.normal_gravity(latitudes, heights))
+
+
+def test_normal_gravity_only_rotation():
+    """
+    Check normal gravity only with rotation (no gravitational attraction)
+    """
+    radius = 1
+    omega = 2
+    sphere = Sphere(
+        name="sphere", radius=radius, geocentric_grav_const=0, angular_velocity=omega
+    )
+    # Check normal gravity on the equator
+    for height in [1, 2, 3, 4]:
+        expected_value = -1e5 * (omega ** 2) * (radius + height)
+        npt.assert_allclose(
+            expected_value, sphere.normal_gravity(latitude=0, height=height),
+        )
+    # Check normal gravity on the poles (must be equal to zero)
+    for height in [1, 2, 3, 4]:
+        assert sphere.normal_gravity(latitude=90, height=height) < 1e-15
+        assert sphere.normal_gravity(latitude=-90, height=height) < 1e-15
+    # Check normal gravity at 45 degrees latitude
+    for height in [1, 2, 3, 4]:
+        expected_value = -1e5 * (omega ** 2) * (radius + height) * np.sqrt(2) / 2
+        npt.assert_allclose(
+            expected_value, sphere.normal_gravity(latitude=45, height=height),
+        )
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 74052c51..225d36a9 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -14,10 +14,12 @@ Reference Ellipsoid
    :toctree: generated/
 
     Ellipsoid
+    Sphere
 
-All ellipsoids are instances of the :class:`~boule.Ellipsoid` class. See the
-class reference documentation for a list its derived physical properties
-(attributes) and computations/transformations that it can perform (methods).
+All ellipsoids are instances of the :class:`~boule.Ellipsoid` or the
+:class:`~boule.Sphere` classes. See the class reference documentation for
+a list their derived physical properties (attributes) and
+computations/transformations that they can perform (methods).
 
 Utilities
 ---------
diff --git a/doc/references.rst b/doc/references.rst
index faeac11e..924b508e 100644
--- a/doc/references.rst
+++ b/doc/references.rst
@@ -2,6 +2,7 @@ References
 ==========
 
 .. [Ardalan2009] Ardalan, A. A., Karimi, R., & Grafarend, E. W. (2009). A New Reference Equipotential Surface, and Reference Ellipsoid for the Planet Mars. Earth, Moon, and Planets, 106(1), 1. https://doi.org/10.1007/s11038-009-9342-7
+.. [Heiskanen-Moritz] Heiskanen, W.A. & Moritz, H. (1967) Physical Geodesy. W.H. Freeman and Company.
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [LiGotze2001] Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi:`10.1190/1.1487109 <https://doi.org/10.1190/1.1487109>`__
 .. [Vermeille2002] Vermeille, H., 2002. Direct transformation from geocentric coordinates to geodetic coordinates. Journal of Geodesy. 76. 451-454. doi:`10.1007/s00190-002-0273-6 <https://doi.org/10.1007/s00190-002-0273-6>`__
diff --git a/tutorials/overview.py b/tutorials/overview.py
index cc2ad885..a540ffbe 100644
--- a/tutorials/overview.py
+++ b/tutorials/overview.py
@@ -62,28 +62,36 @@
 print(bl.MARS.gravity_equator)
 
 ###############################################################################
-# You can also define your own ellipsoid. For example, this would be the
-# definition of a sphere with 1000 m radius and dummy values for :math:`GM` and
-# :math:`\omega`:
-
-sphere = bl.Ellipsoid(
-    name="Sphere",
-    long_name="Ellipsoid with 0 flattening",
-    flattening=0,
+# You can also define your own ellipsoid. For example, this would be a definition of an
+# ellipsoid with 1000 m semimajor axis, flattening equal to 0.5 and dummy values for
+# :math:`GM` and :math:`\omega`:
+
+ellipsoid = bl.Ellipsoid(
+    name="Ellipsoid",
+    long_name="Ellipsoid with 0.5 flattening",
+    flattening=0.5,
     semimajor_axis=1000,
     geocentric_grav_const=1,
     angular_velocity=1,
 )
-print(sphere)
-print(sphere.semiminor_axis)
-print(sphere.first_eccentricity)
+print(ellipsoid)
+print(ellipsoid.semiminor_axis)
+print(ellipsoid.first_eccentricity)
+
 
 ###############################################################################
-# However, the equations for calculating gravity are not suited for the 0
-# flattening case. **So don't define reference spheres like this.** This is due
-# to the first eccentricity being 0 (it appears in divisions in the equations).
+# If the ellipsoid has zero flattening (a sphere), you must use the
+# :class:`boule.Sphere` class instead. For example, this would be the definition of
+# a sphere with 1000 m radius and dummy values for :math:`GM` and :math:`\omega`:
 
-print(sphere.gravity_pole)
+sphere = bl.Sphere(
+    name="Sphere",
+    long_name="Ellipsoid with 0 flattening",
+    radius=1000,
+    geocentric_grav_const=1,
+    angular_velocity=1,
+)
+print(sphere)
 
 ###############################################################################
 # Computations