diff --git a/bilby/gw/eos/eos.py b/bilby/gw/eos/eos.py
index 3faa866697cf3d45f3f9cb2535ee27d2a8db356c..17d99b9e0c350bba3a366372622afa7d6c191812 100644
--- a/bilby/gw/eos/eos.py
+++ b/bilby/gw/eos/eos.py
@@ -34,28 +34,28 @@ valid_eos_dict = dict(zip(valid_eos_names, valid_eos_file_paths))
 
 
 class TabularEOS(object):
-    """Read in EOS tables and interpolate
     """
+    Given a valid eos input format, such as 2-D array, an ascii file, or a string, parse, and interpolate
 
-    def __init__(self, eos, sampling_flag=False, warning_flag=False):
-        """Given a valid eos input format, such as 2-D array, an ascii file, or a string, parse, and interpolate
+    Parameters
+    ----------
+    eos (`numpy.ndarray`, `str`, ASCII TABLE):
+        if `numpy.ndarray` then user supplied pressure-density 2D numpy array.
+        if `str` then given a valid eos name, relevant preshipped ASCII table will be loaded
+            if ASCII TABLE then given viable file extensions, which include .txt,.dat, etc (np.loadtxt used),
+            read in pressure density from file.
 
-        Parameters
-        ----------
-            eos (`numpy.ndarray`, `str`, ASCII TABLE):
-                if `numpy.ndarray` then user supplied pressure-density 2D numpy array.
-                if `str` then given a valid eos name, relevant preshipped ASCII table will be loaded
-                if ASCII TABLE then given viable file extensions, which include .txt,.dat, etc (np.loadtxt used),
-                    read in pressure density from file.
+    sampling_flag (`bool`): Do you plan on sampling the parameterized EOS? Highly recommended. Defaults to False.
 
-            sampling_flag (`bool`): Do you plan on sampling the parameterized EOS?
+    warning_flag (`bool`): Keeps track of status of various physical checks on EoS.
 
-            warning_flag (`bool`):
+    Attributes:
+        msg (str): Human readable string describing the exception.
+        code (int): Exception error code.
+    """
+
+    def __init__(self, eos, sampling_flag=False, warning_flag=False):
 
-        Attributes:
-            msg (str): Human readable string describing the exception.
-            code (int): Exception error code.
-        """
         self.sampling_flag = sampling_flag
         self.warning_flag = warning_flag
 
@@ -137,6 +137,12 @@ class TabularEOS(object):
         """
         Find value of energy_from_pressure
         as in lalsimulation, return e = K * p**(3./5.) below min pressure
+
+        :param pressure (`float`): pressure in geometerized units.
+        :param interp_type (`str`): String specifying which interpolation type to use.
+                                    Currently implemented: 'CubicSpline', 'linear'.
+
+        :param energy_density (`float`): energy-density in geometerized units.
         """
         pressure = np.atleast_1d(pressure)
         energy_returned = np.zeros(pressure.size)
@@ -167,6 +173,12 @@ class TabularEOS(object):
         """
         Find p(h)
         as in lalsimulation, return p = K * h**(5./2.) below min enthalpy
+
+        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
+        :interp_type (`str`): String specifying interpolation type.
+                              Current implementations are 'CubicSpline', 'linear'.
+
+        :return pressure (`float`): pressure in geometerized units.
         """
         pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy)
         pressure_returned = np.zeros(pseudo_enthalpy.size)
@@ -196,6 +208,12 @@ class TabularEOS(object):
         """
         Find energy_density_from_pseudo_enthalpy(pseudo_enthalpy)
         as in lalsimulation, return e = K * h**(3./2.) below min enthalpy
+
+        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
+        :param interp_type (`str`): String specifying interpolation type.
+                                    Current implementations are 'CubicSpline', 'linear'.
+
+        :return energy_density (`float`): energy-density in geometerized units.
         """
         pseudo_enthalpy = np.atleast_1d(pseudo_enthalpy)
         energy_returned = np.zeros(pseudo_enthalpy.size)
@@ -224,6 +242,12 @@ class TabularEOS(object):
         """
         Find h(epsilon)
         as in lalsimulation, return h = K * e**(2./3.) below min enthalpy
+
+        :param energy_density (`float`): energy-density in geometerized units.
+        :param interp_type (`str`): String specifying interpolation type.
+                                    Current implementations are 'CubicSpline', 'linear'.
+
+        :return pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
         """
         energy_density = np.atleast_1d(energy_density)
         pseudo_enthalpy_returned = np.zeros(energy_density.size)
@@ -252,6 +276,14 @@ class TabularEOS(object):
     def dedh(self, pseudo_enthalpy, rel_dh=1e-5, interp_type='CubicSpline'):
         """
         Value of [depsilon/dh](p)
+
+        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
+        :param interp_type (`str`): String specifying interpolation type.
+                                    Current implementations are 'CubicSpline', 'linear'.
+        :param rel_dh (`float`): Relative step size in pseudo-enthalpy space.
+
+        :return dedh (`float`): Derivative of energy-density with respect to pseudo-enthalpy
+                                evaluated at `pseudo_enthalpy` in geometerized units.
         """
 
         # step size=fraction of value
@@ -265,6 +297,14 @@ class TabularEOS(object):
     def dedp(self, pressure, rel_dp=1e-5, interp_type='CubicSpline'):
         """
         Find value of [depsilon/dp](p)
+
+        :param pressure (`float`): pressure in geometerized units.
+        :param rel_dp (`float`): Relative step size in pressure space.
+        :param interp_type (`float`): String specifying interpolation type.
+                                      Current implementations are 'CubicSpline', 'linear'.
+
+        :return dedp (`float`): Derivative of energy-density with respect to pressure
+                                evaluated at `pressure`.
         """
 
         # step size=fraction of value
@@ -282,6 +322,12 @@ class TabularEOS(object):
 
         Assumes the equation
         vs = c (de/dp)^{-1/2}
+
+        :param pseudo_enthalpy (`float`): Dimensionless pseudo-enthalpy.
+        :param interp_type (`str`): String specifying interpolation type.
+                                    Current implementations are 'CubicSpline', 'linear'.
+
+        :return v_s (`float`): Speed of sound at `pseudo-enthalpy` in geometerized units.
         """
         pressure = self.pressure_from_pseudo_enthalpy(pseudo_enthalpy, interp_type=interp_type)
         return self.dedp(pressure, interp_type=interp_type) ** -0.5
@@ -341,6 +387,28 @@ class TabularEOS(object):
         self.pseudo_enthalpy = np.array(hdat)
 
     def plot(self, rep, xlim=None, ylim=None, units=None):
+        """
+        Given a representation in the form 'energy_density-pressure', plot the EoS in that space.
+
+        Parameters
+        ----------
+        rep: str
+            Representation to plot. For example, plotting in energy_density-pressure space,
+            specify 'energy_density-pressure'
+        xlim: list
+            Plotting bounds for x-axis in the form [low, high].
+            Defaults to 'None' which will plot from 10% below min x value to 10% above max x value
+        ylim: list
+            Plotting bounds for y-axis in the form [low, high].
+            Defaults to 'None' which will plot from 10% below min y value to 10% above max y value
+        units: str
+            Specifies unit system to plot. Currently can plot in CGS:'cgs', SI:'si', or geometerized:'geom'
+
+        Returns
+        -------
+        fig: matplotlib.figure.Figure
+            EOS plot.
+        """
 
         # Set data based on specified representation
         varnames = rep.split('-')
@@ -418,34 +486,34 @@ def spectral_adiabatic_index(gammas, x):
 
 
 class SpectralDecompositionEOS(TabularEOS):
+    """
+    Parameterized EOS using a spectral
+    decomposition per Lindblom
+    arXiv: 1009.0738v2. Inherits from TabularEOS.
+
+    Parameters
+    ----------
+    gammas: list
+        List of adiabatic expansion parameters used
+        to construct the equation of state in various
+        spaces.
+    p0: float
+        The starting point in pressure of the high-density EoS. This is stitched to
+        the low-density portion of the SLY EoS model. The default value chosen is set to
+        a sufficiently low pressure so that the high-density EoS will never be
+        overconstrained.
+    e0/c**2: float
+        The starting point in energy-density of the high-density EoS. This is stitched to
+        the low-density portion of the SLY EoS model. The default value chosen is set to
+        a sufficiently low energy density so that the high-density EoS will never be
+        overconstrained.
+    xmax: float
+        highest dimensionless pressure value in EoS
+    npts: float (optional)
+        number of points in pressure-energy density data.
+    """
 
     def __init__(self, gammas, p0=3.01e33, e0=2.03e14, xmax=None, npts=100, sampling_flag=False, warning_flag=False):
-        """
-        Parameterized EOS using a spectral
-        decomposition per Lindblom
-        arXiv: 1009.0738v2
-
-        Parameters
-        ----------
-        gammas: list
-            List of adiabatic expansion parameters used
-            to construct the equation of state in various
-            spaces.
-        p0: float
-            The starting point in pressure of the high-density EoS. This is stitched to
-            the low-density portion of the SLY EoS model. The default value chosen is set to
-            a sufficiently low pressure so that the high-density EoS will never be
-            overconstrained.
-        e0/c**2: float
-            The starting point in energy-density of the high-density EoS. This is stitched to
-            the low-density portion of the SLY EoS model. The default value chosen is set to
-            a sufficiently low energy density so that the high-density EoS will never be
-            overconstrained.
-        xmax: float
-            highest dimensionless pressure value in EoS
-        npts: float (optional)
-            number of points in pressure-energy density data.
-        """
         self.warning_flag = warning_flag
         self.gammas = gammas
         self.p0 = p0
@@ -566,16 +634,18 @@ class SpectralDecompositionEOS(TabularEOS):
 
 class EOSFamily(object):
 
-    def __init__(self, eos, npts=500):
-        """Create a EOS family and get mass-radius information
+    """Create a EOS family and get mass-radius information
 
-        Parameters:
-            eos (`object`): Supply a `TabularEOS` class (or subclass)
+    Parameters:
+        eos (`object`): Supply a `TabularEOS` class (or subclass)
+        npts (`float`): Number of points to calculate for mass-radius relation.
+                        Default is 500.
 
-        Note:
-            The mass-radius and mass-k2 data should be
-            populated here via the TOV solver upon object construction.
+    Note:
+        The mass-radius and mass-k2 data should be
+        populated here via the TOV solver upon object construction.
         """
+    def __init__(self, eos, npts=500):
         self.eos = eos
 
         # FIXME: starting_energy_density is set somewhat arbitrarily