diff --git a/lenstronomy/LightModel/Profiles/hernquist.py b/lenstronomy/LightModel/Profiles/hernquist.py
index 7a33bb4a6..ce7235339 100644
--- a/lenstronomy/LightModel/Profiles/hernquist.py
+++ b/lenstronomy/LightModel/Profiles/hernquist.py
@@ -1,10 +1,11 @@
 import lenstronomy.Util.param_util as param_util
+import numpy as np
 
 __all__ = ["Hernquist", "HernquistEllipse"]
 
 
 class Hernquist(object):
-    """Class for pseudo Jaffe lens light (2d projected light/mass distribution."""
+    """Class for Hernquist lens light (2d projected light/mass distribution)."""
 
     def __init__(self):
         from lenstronomy.LensModel.Profiles.hernquist import Hernquist as Hernquist_lens
@@ -27,13 +28,13 @@ def __init__(self):
     def function(self, x, y, amp, Rs, center_x=0, center_y=0):
         """
 
-        :param x:
-        :param y:
-        :param amp:
+        :param x: x-position
+        :param y: y-position
+        :param amp: surface brightness amplitude
         :param Rs: scale radius: half-light radius = Rs / 0.551
-        :param center_x:
-        :param center_y:
-        :return:
+        :param center_x: centroid in x-direction
+        :param center_y: centroid in y-direction
+        :return: surface brightness
         """
         rho0 = self.lens.sigma2rho(amp, Rs)
         return self.lens.density_2d(x, y, rho0, Rs, center_x, center_y)
@@ -41,14 +42,26 @@ def function(self, x, y, amp, Rs, center_x=0, center_y=0):
     def light_3d(self, r, amp, Rs):
         """
 
-        :param r:
-        :param amp:
-        :param Rs:
+        :param r: 3d radius (in angular units)
+        :param amp: surface brightness amplitude
+        :param Rs: scale radius: half-light radius = Rs / 0.551
         :return:
         """
         rho0 = self.lens.sigma2rho(amp, Rs)
         return self.lens.density(r, rho0, Rs)
 
+    @staticmethod
+    def total_flux(amp, Rs):
+        """
+
+        :param amp: surface brightness amplitude
+        :param Rs: scale radius: half-light radius = Rs / 0.551
+        :return: total integrated flux of profile
+        """
+        rhos = amp / Rs
+        m_tot = 2 * np.pi * rhos * Rs**3
+        return m_tot
+
 
 class HernquistEllipse(object):
     """Class for elliptical pseudo Jaffe lens light (2d projected light/mass
@@ -81,15 +94,15 @@ def __init__(self):
     def function(self, x, y, amp, Rs, e1, e2, center_x=0, center_y=0):
         """
 
-        :param x:
-        :param y:
-        :param amp:
-        :param Rs:
-        :param e1:
-        :param e2:
-        :param center_x:
-        :param center_y:
-        :return:
+        :param x: x-position
+        :param y: y-position
+        :param amp: surface brightness amplitude
+        :param Rs: scale radius: half-light radius = Rs / 0.551
+        :param e1: eccentricity component
+        :param e2: eccentricity component
+        :param center_x: centroid in x-direction
+        :param center_y: centroid in y-direction
+        :return: surface brightness
         """
         x_, y_ = param_util.transform_e1e2_product_average(
             x, y, e1, e2, center_x, center_y
@@ -99,12 +112,23 @@ def function(self, x, y, amp, Rs, e1, e2, center_x=0, center_y=0):
     def light_3d(self, r, amp, Rs, e1=0, e2=0):
         """
 
-        :param r:
-        :param amp:
-        :param Rs:
-        :param e1:
-        :param e2:
-        :return:
+        :param r: 3d radius (in angular units)
+        :param amp: surface brightness amplitude
+        :param Rs: scale radius: half-light radius = Rs / 0.551
+        :param e1: eccentricity component
+        :param e2: eccentricity component
+        :return: flux density in 3d
         """
         rho0 = self.lens.sigma2rho(amp, Rs)
         return self.lens.density(r, rho0, Rs)
+
+    def total_flux(self, amp, Rs, e1=0, e2=0):
+        """
+
+        :param amp: surface brightness amplitude
+        :param Rs: scale radius: half-light radius = Rs / 0.551
+        :param e1: eccentricity component
+        :param e2: eccentricity component
+        :return: total integrated flux
+        """
+        return self.spherical.total_flux(amp=amp, Rs=Rs)
diff --git a/test/test_LightModel/test_Profiles/test_hernquist.py b/test/test_LightModel/test_Profiles/test_hernquist.py
new file mode 100644
index 000000000..a96c45d19
--- /dev/null
+++ b/test/test_LightModel/test_Profiles/test_hernquist.py
@@ -0,0 +1,24 @@
+from lenstronomy.LightModel.Profiles.hernquist import Hernquist, HernquistEllipse
+from lenstronomy.Util.util import make_grid
+import numpy as np
+import numpy.testing as npt
+
+
+class TestHernquist(object):
+
+    def setup_method(self):
+        self.hernquist = Hernquist()
+        self.hernquist_ellipse = HernquistEllipse()
+
+    def test_total_flux(self):
+        delta_pix = 0.2
+        x, y = make_grid(numPix=1000, deltapix=delta_pix)
+
+        rs, amp = 1, 1
+        total_flux = self.hernquist.total_flux(amp=amp, Rs=rs)
+        flux = self.hernquist.function(x, y, amp=amp, Rs=rs)
+        total_flux_numerics = np.sum(flux) * delta_pix**2
+        npt.assert_almost_equal(total_flux_numerics / total_flux, 1, decimal=1)
+
+        total_flux_ellipse = self.hernquist_ellipse.total_flux(amp=amp, Rs=rs)
+        npt.assert_almost_equal(total_flux_ellipse, total_flux)