diff --git a/.gitignore b/.gitignore
index 195140c..6fe18fd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -49,6 +49,8 @@ cover/
 
 # Sphinx documentation
 docs/_build/
+docs/source/api
+docs/source/modules
 
 # Jupyter Notebook
 **/*.ipynb_checkpoints
diff --git a/docs/README.rst b/docs/README.rst
new file mode 100644
index 0000000..0601145
--- /dev/null
+++ b/docs/README.rst
@@ -0,0 +1,22 @@
+SSAPy Documentation
+===================
+
+Building the docs
+-----------------
+
+First, build and _install_ `ssapy` following the `Installing SSAPy <https://LLNL.github.io/SSAPy/installation.html>`_ section of the documentation.
+
+Then, to build the HTML documentation locally, from the `docs` directory, run
+
+.. code-block:: bash
+
+    make html
+
+(run it twice to generate all files.)
+
+Then open `_build/html/index.html` to browse the docs locally.
+
+Alternatively, you can run `sphinx-autobuild . _build/html` to start a server that watches for changes in `/docs`
+and regenerates the HTML docs automatically while serving them at http://127.0.0.1:8000/.
+
+Note that if you updated docstrings, you'll need to re-build and re-install ssapy before re-building the docs.
diff --git a/docs/source/installation.rst b/docs/source/installation.rst
index 7529184..7c037f4 100644
--- a/docs/source/installation.rst
+++ b/docs/source/installation.rst
@@ -34,7 +34,7 @@ These requirements can be easily installed on most modern macOS and Linux system
        .. code-block:: console
 
           brew update
-          brew install gcc git git-lfs python3 svn graphviz
+          brew install gcc git git-lfs python3 graphviz
 
 Installation
 ------------
diff --git a/ssapy/accel.py b/ssapy/accel.py
index cceabc6..8a15282 100644
--- a/ssapy/accel.py
+++ b/ssapy/accel.py
@@ -1,3 +1,7 @@
+"""
+Classes for modeling accelerations.
+"""
+
 import numpy as np
 
 from . import _ssapy
@@ -101,6 +105,7 @@ def __eq__(self, rhs):
 
 
 class AccelProd(Accel):
+    """Acceleration defined as the product of an acceleration with a constant factor."""
     def __init__(self, accel, factor):
         super().__init__()
         self.accel = accel
diff --git a/ssapy/body.py b/ssapy/body.py
index 1c3a344..5058d10 100644
--- a/ssapy/body.py
+++ b/ssapy/body.py
@@ -1,3 +1,7 @@
+"""
+Classes representing celestial bodies.
+"""
+
 import erfa
 import numpy as np
 from .utils import _gpsToTT, iers_interp
diff --git a/ssapy/constants.py b/ssapy/constants.py
index 776e021..668dbed 100644
--- a/ssapy/constants.py
+++ b/ssapy/constants.py
@@ -1,3 +1,7 @@
+"""
+A collection of physical constants.
+"""
+
 import numpy as np
 
 # GM
diff --git a/ssapy/correlate_tracks.py b/ssapy/correlate_tracks.py
index fbefddf..d0cd27a 100644
--- a/ssapy/correlate_tracks.py
+++ b/ssapy/correlate_tracks.py
@@ -2012,7 +2012,7 @@ def summarize_tracklet(arc):
     pmra, pmdec: the proper motion in right ascension and declination
     dpmra, dpmdec: the uncertainty in the proper motion in RA and declination
     """
-    unit = ssapy.utils.lb2unit(arc['ra'].to(u.rad).value,
+    unit = ssapy.utils.lb_to_unit(arc['ra'].to(u.rad).value,
                                arc['dec'].to(u.rad).value)
     meanunit = np.mean(unit, axis=0)
     # defines tangent plane
@@ -2033,7 +2033,7 @@ def summarize_tracklet(arc):
     rasig = np.sqrt(np.diag(racovar))
     decp, deccovar = np.polyfit(dt, dectp, 1, w=1./dpos, cov='unscaled')
     decsig = np.sqrt(np.diag(deccovar))
-    meanra, meandec = ssapy.utils.unit2lb(meanunit[None, :])
+    meanra, meandec = ssapy.utils.unit_to_lb(meanunit[None, :])
     rap[1] += meanra[0]
     decp[1] += meandec[0]
     out = [(rap[1]*u.rad, decp[1]*u.rad),
diff --git a/ssapy/ellipsoid.py b/ssapy/ellipsoid.py
index e510267..3be1d43 100644
--- a/ssapy/ellipsoid.py
+++ b/ssapy/ellipsoid.py
@@ -1,3 +1,11 @@
+"""
+Class to handle transformations between ECEF x,y,z coords and geodetic
+longitude, latitude, and height.
+
+Technically, only handles a one-axis ellipsoid, defined via a flattening
+parameter f, but that's good enough for simple Earth models.
+"""
+
 import numpy as np
 
 from ._ssapy import Ellipsoid
diff --git a/ssapy/gravity.py b/ssapy/gravity.py
index ecf9443..451bbf7 100644
--- a/ssapy/gravity.py
+++ b/ssapy/gravity.py
@@ -1,8 +1,12 @@
+"""
+Classes for gravity-related accelerations.
+"""
+
 # from functools import lru_cache
 
 import numpy as np
 
-from .accel import Accel, _invnorm
+from .accel import Accel as _Accel, _invnorm
 from .utils import find_file, norm
 from . import _ssapy
 
@@ -220,7 +224,7 @@ def __eq__(self, rhs):
         )
 
 
-class AccelThirdBody(Accel):
+class AccelThirdBody(_Accel):
     """Acceleration due to a third body.
 
     Parameters
@@ -254,7 +258,7 @@ def __call__(self, r, v, t, **kwargs):
         return -self.body.mu * (s / norm(s)**3 - d / norm(d)**3)
 
 
-class AccelHarmonic(Accel):
+class AccelHarmonic(_Accel):
     """Acceleration due to a harmonic potential.
 
     Parameters
diff --git a/ssapy/io.py b/ssapy/io.py
index 0927a60..12a1399 100644
--- a/ssapy/io.py
+++ b/ssapy/io.py
@@ -1,3 +1,7 @@
+"""
+Collection of functions to read from and write to various file formats.
+"""
+
 import datetime
 import numpy as np
 from astropy.time import Time
diff --git a/ssapy/orbit.py b/ssapy/orbit.py
index ff240a0..f481c0e 100644
--- a/ssapy/orbit.py
+++ b/ssapy/orbit.py
@@ -1,18 +1,3 @@
-import numpy as np
-from astropy.time import Time
-from .utils import (
-    normSq, normed, lazy_property, unitAngle3, sunPos, _gpsToTT, _wrapToPi,
-    teme_to_gcrf, gcrf_to_teme, iers_interp
-)
-from .constants import EARTH_MU
-from .propagator import KeplerianPropagator
-
-try:
-    import erfa
-except ImportError:
-    # Let this raise
-    import astropy._erfa as erfa
-
 """
 Module to handle Earth satellite state vector manipulation and propagation.
 
@@ -27,6 +12,22 @@
 """
 
 
+import numpy as np
+import astropy
+from .utils import (
+    normSq, normed, lazy_property as _lazy_property, unitAngle3, sunPos,
+    _gpsToTT, _wrapToPi, teme_to_gcrf, gcrf_to_teme, iers_interp
+)
+from .constants import EARTH_MU
+from .propagator import KeplerianPropagator as _KeplerianPropagator
+
+try:
+    import erfa
+except ImportError:
+    # Let this raise
+    import astropy._erfa as erfa
+
+
 # Conversion routines for anomalies
 def _ellipticalEccentricToTrueAnomaly(E, e):
     """Compute true anomaly from eccentric anomaly for elliptical orbit.
@@ -667,7 +668,7 @@ class Orbit:
     # Internally use equinoctial elements, calculating cartesian and keplerian
     # on demand.  Main ctor uses cartesian though
     def __init__(self, r, v, t, mu=EARTH_MU, propkw=None):
-        if isinstance(t, Time):
+        if isinstance(t, astropy.time.Time):
             t = t.gps
         self.r, self.v = np.broadcast_arrays(r, v)
         self.propkw = dict() if propkw is None else propkw
@@ -778,7 +779,7 @@ def fromEquinoctialElements(
             The Orbit with given parameters.
         """
         import numbers
-        if isinstance(t, Time):
+        if isinstance(t, astropy.time.Time):
             t = t.gps
         if not all(isinstance(x, numbers.Real) for x in (a, hx, hy, ex, ey, lv, t)):
             a, hx, hy, ex, ey, lv, t = np.broadcast_arrays(
@@ -849,7 +850,7 @@ def fromKeplerianElements(
             The Orbit with given parameters.
         """
         import numbers
-        if isinstance(t, Time):
+        if isinstance(t, astropy.time.Time):
             t = t.gps
         if not all(isinstance(x, numbers.Real) for x in (a, e, i, pa, raan, trueAnomaly, t)):
             a, e, i, pa, raan, trueAnomaly, t = np.broadcast_arrays(
@@ -945,10 +946,10 @@ def fromKozaiMeanKeplerianElements(
             )
         meanMotion = np.sqrt(mu / np.abs(a**3)) * 60.0  # rad/min
 
-        if not isinstance(t, Time):
-            t = Time(t, format='gps')
+        if not isinstance(t, astropy.time.Time):
+            t = astropy.time.Time(t, format='gps')
         mjd = t.mjd
-        # epoch = mjd - Time("1949-12-31T00:00:00").mjd
+        # epoch = mjd - astropy.time.Time("1949-12-31T00:00:00").mjd
         epoch = mjd - 33281.0
 
         sat = Satrec()
@@ -978,7 +979,7 @@ def fromKozaiMeanKeplerianElements(
         v *= 1e3  # km/s -> m/s
         return Orbit(r, v, t, mu=mu, propkw=propkw)
 
-    @lazy_property
+    @_lazy_property
     def kozaiMeanKeplerianElements(self):
         """Kozai mean Keplerian elements in TEME frame
         (a, e, i, pa, raan, trueAnomaly)
@@ -1077,7 +1078,7 @@ def fromTLETuple(cls, tle, propkw=None):
         obj.propkw = propkw if propkw is not None else dict()
         return obj
 
-    def at(self, t, propagator=KeplerianPropagator()):
+    def at(self, t, propagator=_KeplerianPropagator()):
         """Propagate this orbit to time t.
 
         Parameters
@@ -1451,69 +1452,69 @@ def _setPQKeplerian(self):
             -crsp - cosI * srcp, -srsp + cosI * crcp, sinI * cosPa
         ]).T
 
-    @lazy_property
+    @_lazy_property
     def _pEq(self):
         self._setPQEquinoctial()
         return self._pEq
 
-    @lazy_property
+    @_lazy_property
     def _qEq(self):
         self._setPQEquinoctial()
         return self._qEq
 
-    @lazy_property
+    @_lazy_property
     def _pK(self):
         self._setPQKeplerian()
         return self._pK
 
-    @lazy_property
+    @_lazy_property
     def _qK(self):
         self._setPQKeplerian()
         return self._qK
 
-    @lazy_property
+    @_lazy_property
     def a(self):
         """Semimajor axis in meters.
         """
         self._setEquinoctial()
         return self.a
 
-    @lazy_property
+    @_lazy_property
     def hx(self):
         """First component of equinoctial inclination vector.
         """
         self._setEquinoctial()
         return self.hx
 
-    @lazy_property
+    @_lazy_property
     def hy(self):
         """Second component of equinoctial inclination vector.
         """
         self._setEquinoctial()
         return self.hy
 
-    @lazy_property
+    @_lazy_property
     def ex(self):
         """First component of equinoctial eccentricity vector.
         """
         self._setEquinoctial()
         return self.ex
 
-    @lazy_property
+    @_lazy_property
     def ey(self):
         """Second component of equinoctial eccentricity vector.
         """
         self._setEquinoctial()
         return self.ey
 
-    @lazy_property
+    @_lazy_property
     def lv(self):
         """True longitude in radians.
         """
         self._setEquinoctial()
         return self.lv
 
-    @lazy_property
+    @_lazy_property
     def lE(self):
         """Eccentric longitude in radians.
         """
@@ -1540,7 +1541,7 @@ def lE(self):
                 )
             return lE
 
-    @lazy_property
+    @_lazy_property
     def lM(self):
         """Mean longitude in radians.
         """
@@ -1567,48 +1568,48 @@ def lM(self):
                 )
             return lM
 
-    @lazy_property
+    @_lazy_property
     def e(self):
         """Eccentricity.
         """
         return np.sqrt(self.ex * self.ex + self.ey * self.ey)
 
-    @lazy_property
+    @_lazy_property
     def i(self):
         """Inclination in radians.
         """
         self._setKeplerian()
         return self.i
 
-    @lazy_property
+    @_lazy_property
     def pa(self):
         """Periapsis argument in radians.
         """
         self._setKeplerian()
         return self.pa
 
-    @lazy_property
+    @_lazy_property
     def lonPa(self):
         """Longitude of periapsis argument in radians.
         """
         self._setEquinoctial()
         return self.lonPa
 
-    @lazy_property
+    @_lazy_property
     def raan(self):
         """Right ascension of the ascending node in radians.
         """
         self._setKeplerian()
         return self.raan
 
-    @lazy_property
+    @_lazy_property
     def trueAnomaly(self):
         """True anomaly in radians.
         """
         self._setEquinoctial()
         return self.trueAnomaly
 
-    @lazy_property
+    @_lazy_property
     def eccentricAnomaly(self):
         """Eccentric anomaly in radians.
         """
@@ -1635,7 +1636,7 @@ def eccentricAnomaly(self):
                 )
             return eccentricAnomaly
 
-    @lazy_property
+    @_lazy_property
     def meanAnomaly(self):
         """Mean anomaly in radians.
         """
@@ -1670,7 +1671,7 @@ def meanAnomaly(self):
                 )
             return meanAnomaly
 
-    @lazy_property
+    @_lazy_property
     def period(self):
         """Orbital period in seconds.
         """
@@ -1686,7 +1687,7 @@ def period(self):
                 period[~wBound] = np.inf
             return period
 
-    @lazy_property
+    @_lazy_property
     def meanMotion(self):
         """Mean motion in radians per second.
         """
@@ -1704,31 +1705,31 @@ def keplerianElements(self):
         """
         return self.a, self.e, self.i, self.pa, self.raan, self.trueAnomaly
 
-    @lazy_property
+    @_lazy_property
     def p(self):
         """Semi-latus rectum in meters.
         """
         return self.a * (1.0 - self.e**2)
 
-    @lazy_property
+    @_lazy_property
     def angularMomentum(self):
         """(Specific) angular momentum vector in m^2/s.
         """
         return np.cross(self.r, self.v)
 
-    @lazy_property
+    @_lazy_property
     def energy(self):
         """(Specific) orbital energy in m^2/s^2.
         """
         return -0.5 * self.mu / self.a
 
-    @lazy_property
+    @_lazy_property
     def LRL(self):
         """Laplace-Runge-Lenz vector in m^3/s^2.
         """
         return -self.mu * normed(self.r) - np.cross(self.angularMomentum, self.v)
 
-    @lazy_property
+    @_lazy_property
     def periapsis(self):
         """Periapsis coordinate of orbit in meters.
         """
@@ -1739,7 +1740,7 @@ def periapsis(self):
         else:
             return normed(self.LRL) * pdist[:, None]
 
-    @lazy_property
+    @_lazy_property
     def apoapsis(self):
         """Apoapsis coordinate of orbit in meters.
         """
@@ -1759,7 +1760,7 @@ def apoapsis(self):
                 apoapsis[~wBound] = np.array([np.inf, np.inf, np.inf])
             return apoapsis
 
-    @lazy_property
+    @_lazy_property
     def tle(self):
         """Two line element for this Orbit.
         """
@@ -1844,7 +1845,7 @@ def getRV(self, time):
 
         if self.fast:
             import numbers
-            if isinstance(time, Time):
+            if isinstance(time, astropy.time.Time):
                 time = time.gps
             if isinstance(time, numbers.Real):
                 scalar = True
@@ -1879,8 +1880,8 @@ def getRV(self, time):
                 return outR[0], outV[0]
             return outR, outV
         else:
-            if not isinstance(time, Time):
-                time = Time(time, format='gps')
+            if not isinstance(time, astropy.time.Time):
+                time = astropy.time.Time(time, format='gps')
             r, v = self._location.get_gcrs_posvel(time)
             return r.xyz.to(u.m).T.value, v.xyz.to(u.m / u.s).T.value
 
@@ -1898,7 +1899,7 @@ def sunAlt(self, time):
         alt : array_like(n,)
             Altitude of sun in radians.
         """
-        if isinstance(time, Time):
+        if isinstance(time, astropy.time.Time):
             time = time.gps
         ro, _ = self.getRV(time)
         r_sun = sunPos(time, fast=False)
@@ -1926,13 +1927,13 @@ class OrbitalObserver:
     getRV(time)
         Get position and velocity at specified time(s).
     """
-    def __init__(self, orbit, propagator=KeplerianPropagator()):
+    def __init__(self, orbit, propagator=_KeplerianPropagator()):
         self.orbit = orbit
         self.propagator = propagator
 
     def __repr__(self):
         out = "OrbitalObserver({!r}".format(self.orbit)
-        if self.propagator != KeplerianPropagator():
+        if self.propagator != _KeplerianPropagator():
             out += ", propagator={!r}".format(self.propagator)
         out += ")"
         return out
diff --git a/ssapy/orbit_solver.py b/ssapy/orbit_solver.py
index 1ced042..df50c14 100644
--- a/ssapy/orbit_solver.py
+++ b/ssapy/orbit_solver.py
@@ -1,5 +1,5 @@
 """
-Classes to solve for Keplerina orbital parameters from different initial inputs
+Classes to solve for Keplerian orbital parameters from different initial inputs
 
 # TODO: Provide overview of different solvers
 
diff --git a/ssapy/particles.py b/ssapy/particles.py
index 0e20bb7..f753025 100644
--- a/ssapy/particles.py
+++ b/ssapy/particles.py
@@ -1,5 +1,9 @@
+"""
+Classes for sampling orbit model parameters.
+"""
+
 import numpy as np
-from .orbit import Orbit
+from .orbit import Orbit as _Orbit
 from .compute import rv
 from . import utils
 
@@ -113,7 +117,7 @@ def orbits(self):
         """ Get a list of orbits corresponding to each stored particle
         """
         if self._orbits is None:
-            self._orbits = [Orbit(p[0:3], p[3:6], t=self.epoch) for p in self.particles]
+            self._orbits = [_Orbit(p[0:3], p[3:6], t=self.epoch) for p in self.particles]
         return self._orbits
 
     @property
@@ -260,7 +264,7 @@ def draw_orbit(self):
         prob = utils.get_normed_weights(self.ln_wts)
         particle_ndx = np.random.choice(self.num_particles, size=1, p=prob)
         p = self.particles[particle_ndx, :][0]
-        return [Orbit(p[0:3], p[3:6], t=self.epoch)]
+        return [_Orbit(p[0:3], p[3:6], t=self.epoch)]
 
     def mean(self):
         """ Evaluate the weighted mean of all particles
diff --git a/ssapy/plotUtils.py b/ssapy/plotUtils.py
index 2030440..ee9a6bb 100644
--- a/ssapy/plotUtils.py
+++ b/ssapy/plotUtils.py
@@ -1,3 +1,8 @@
+"""
+Utility functions for plotting.
+"""
+
+
 from .body import get_body
 from .compute import groundTrack
 from .constants import RGEO, EARTH_RADIUS
diff --git a/ssapy/propagator.py b/ssapy/propagator.py
index 9196651..f0f6fb7 100644
--- a/ssapy/propagator.py
+++ b/ssapy/propagator.py
@@ -1,3 +1,7 @@
+"""
+Classes for propagating orbits.
+"""
+
 from abc import ABC, abstractmethod
 import numpy as np
 from astropy.time import Time
diff --git a/ssapy/rvsampler.py b/ssapy/rvsampler.py
index e09c58c..efbfc82 100644
--- a/ssapy/rvsampler.py
+++ b/ssapy/rvsampler.py
@@ -592,7 +592,7 @@ def vR(r, ub=np.inf, verbose=False, square=False):
     rRange = sol[0]
     vlos = sgn*np.sqrt(vr2(rRange))
     import ssapy.utils
-    ra, dec = ssapy.utils.unit2lb(np.array([lineOfSight]))
+    ra, dec = ssapy.utils.unit_to_lb(np.array([lineOfSight]))
     rGuess, vGuess = radecRateObsToRV(
         ra[0], dec[0], rRange, muAlpha, muDelta, vlos, rStation, vGroundSky)
     # ignores light-time correction, etc.
diff --git a/tests/test_orbit.py b/tests/test_orbit.py
index c994cf8..27983af 100644
--- a/tests/test_orbit.py
+++ b/tests/test_orbit.py
@@ -14,7 +14,7 @@
 from ssapy.orbit import _hyperbolicEccentricToTrueLongitude, _hyperbolicTrueToEccentricLongitude
 from ssapy.orbit import _ellipticalEccentricToMeanLongitude, _ellipticalMeanToEccentricLongitude
 from ssapy.orbit import _hyperbolicEccentricToMeanLongitude, _hyperbolicMeanToEccentricLongitude
-from ssapy.utils import normed, norm, teme2gcrf
+from ssapy.utils import normed, norm, teme_to_gcrf
 from ssapy_test_helpers import timer, checkAngle, checkSphere
 
 
@@ -1682,7 +1682,7 @@ def test_sgp4():
             times.append(Time(timestr))
         times = Time(times)
 
-        rot = teme2gcrf(times[0])
+        rot = teme_to_gcrf(times[0])
         x = np.cos(np.deg2rad(b3['azimuthAngle'])) * np.cos(np.deg2rad(b3['polarAngle']))
         y = np.sin(np.deg2rad(b3['azimuthAngle'])) * np.cos(np.deg2rad(b3['polarAngle']))
         z = np.sin(np.deg2rad(b3['polarAngle']))
@@ -1699,7 +1699,7 @@ def test_sgp4():
             vsgp4.append(v)
         rsgp4 = np.array(rsgp4)*1e3
         vsgp4 = np.array(vsgp4)*1e3
-        rot = teme2gcrf(orbit.t)
+        rot = teme_to_gcrf(orbit.t)
         rsgp4 = np.dot(rot, rsgp4.T).T
         vsgp4 = np.dot(rot, vsgp4.T).T
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 2b2141e..ecd3042 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -82,27 +82,27 @@ def test_angular_conversions():
     np.random.seed(seed)
     npts = 10000
     uv = normed(np.random.randn(npts, 3))
-    lb = utils.unit2lb(uv)
-    tp = utils.unit2tp(uv)
+    lb = utils.unit_to_lb(uv)
+    tp = utils.unit_to_tp(uv)
     # round trips
     # 3 systems, back and forth to all other systems -> 6 tests.
     np.testing.assert_allclose(uv,
-                               utils.lb2unit(*utils.unit2lb(uv)),
+                               utils.lb_to_unit(*utils.unit_to_lb(uv)),
                                rtol=0, atol=1e-10)
     np.testing.assert_allclose(uv,
-                               utils.tp2unit(*utils.unit2tp(uv)),
+                               utils.tp_to_unit(*utils.unit_to_tp(uv)),
                                rtol=0, atol=1e-10)
     np.testing.assert_allclose(np.concatenate(tp),
-                               np.concatenate(utils.unit2tp(utils.tp2unit(*tp))),
+                               np.concatenate(utils.unit_to_tp(utils.tp_to_unit(*tp))),
                                rtol=0, atol=1e-10)
     np.testing.assert_allclose(np.concatenate(tp),
-                               np.concatenate(utils.lb2tp(*utils.tp2lb(*tp))),
+                               np.concatenate(utils.lb_to_tp(*utils.tp_to_lb(*tp))),
                                rtol=0, atol=1e-10)
     np.testing.assert_allclose(np.concatenate(lb),
-                               np.concatenate(utils.unit2lb(utils.lb2unit(*lb))),
+                               np.concatenate(utils.unit_to_lb(utils.lb_to_unit(*lb))),
                                rtol=0, atol=1e-10)
     np.testing.assert_allclose(np.concatenate(lb),
-                               np.concatenate(utils.tp2lb(*utils.lb2tp(*lb))),
+                               np.concatenate(utils.tp_to_lb(*utils.lb_to_tp(*lb))),
                                rtol=0, atol=1e-10)
 
     # check tangent plane round tripping.
@@ -111,33 +111,33 @@ def test_angular_conversions():
     # so we need to make up some lcen, bcen to project from.
     noise = np.random.randn(npts, 3)*0.01
     uv2 = normed(uv + noise)
-    lcen, bcen = utils.unit2lb(uv2)
-    xy = utils.lb2tan(*lb, lcen=lcen, bcen=bcen)
+    lcen, bcen = utils.unit_to_lb(uv2)
+    xy = utils.lb_to_tan(*lb, lcen=lcen, bcen=bcen)
 
     # vector lcen, bcen
     np.testing.assert_allclose(
         np.concatenate(lb),
-        np.concatenate(utils.tan2lb(*utils.lb2tan(*lb, lcen=lcen, bcen=bcen),
+        np.concatenate(utils.tan_to_lb(*utils.lb_to_tan(*lb, lcen=lcen, bcen=bcen),
                                      lcen=lcen, bcen=bcen)))
     np.testing.assert_allclose(
         np.concatenate(xy),
-        np.concatenate(utils.lb2tan(*utils.tan2lb(*xy, lcen=lcen, bcen=bcen),
+        np.concatenate(utils.lb_to_tan(*utils.tan_to_lb(*xy, lcen=lcen, bcen=bcen),
                                      lcen=lcen, bcen=bcen)))
 
     # single lcen, bcen; careful to choose all points to be on same hemisphere
     uv2 = uv.copy()
     uv2[:, 0] = np.abs(uv2[:, 0])
     lcen, bcen = (0, 0)
-    lb2 = utils.unit2lb(uv2)
-    xy2 = utils.lb2tan(*lb2, lcen=lcen, bcen=bcen)
+    lb2 = utils.unit_to_lb(uv2)
+    xy2 = utils.lb_to_tan(*lb2, lcen=lcen, bcen=bcen)
 
     np.testing.assert_allclose(
         np.concatenate(lb2),
-        np.concatenate(utils.tan2lb(*utils.lb2tan(*lb2, lcen=lcen, bcen=bcen),
+        np.concatenate(utils.tan_to_lb(*utils.lb_to_tan(*lb2, lcen=lcen, bcen=bcen),
                                      lcen=lcen, bcen=bcen)))
     np.testing.assert_allclose(
         np.concatenate(xy2),
-        np.concatenate(utils.lb2tan(*utils.tan2lb(*xy2, lcen=lcen, bcen=bcen),
+        np.concatenate(utils.lb_to_tan(*utils.tan_to_lb(*xy2, lcen=lcen, bcen=bcen),
                                      lcen=lcen, bcen=bcen)))