diff --git a/pymap3d/__init__.py b/pymap3d/__init__.py
index 2b733a1..3767edd 100644
--- a/pymap3d/__init__.py
+++ b/pymap3d/__init__.py
@@ -52,6 +52,10 @@
     from .enu import enu2geodetic, geodetic2enu, aer2enu, enu2aer  # noqa: F401
     from .aer import ecef2aer, aer2ecef, geodetic2aer, aer2geodetic, eci2aer, aer2eci  # noqa: F401
     from .los import lookAtSpheroid  # noqa: F401
-    from .lox import isometric, meridian_dist, loxodrome_inverse  # noqa: F401
+    from .lox import loxodrome_direct, loxodrome_inverse, meridianarc, departure, meanm
+    from .latitude import geodetic2isometric, isometric2geodetic, geodetic2rectifying, rectifying2geodetic, geodetic2conformal, conformal2geodetic, geodetic2parametric, parametric2geodetic, geodetic2geocentric, geocentric2geodetic, geodetic2authalic, authalic2geodetic
+    from .rcurve import rcurve_meridian, rcurve_parallel, rcurve_transverse
+    from .rsphere import rsphere_authalic, rsphere_biaxial, rsphere_curve, rsphere_eqavol, rsphere_euler, rsphere_rectifying, rsphere_triaxial
+    from .util import wrapToPi, wrapTo2Pi, sph2cart, cart2sph, pol2cart, cart2pol
 else:  # pure Python only
     from .math import *  # type: ignore  # noqa: F401, F403
diff --git a/pymap3d/ecef.py b/pymap3d/ecef.py
index 7814d7f..f7961b7 100644
--- a/pymap3d/ecef.py
+++ b/pymap3d/ecef.py
@@ -126,7 +126,7 @@ def ecef2geodetic(x: float, y: float, z: float,
     eps = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E**2) * sin(Beta)) / (ell.semimajor_axis * huE * 1 / cos(Beta) - E**2 * cos(Beta))
 
     Beta += eps
-# %% final output
+#  final output
     lat = arctan(ell.semimajor_axis / ell.semiminor_axis * tan(Beta))
 
     lon = arctan2(y, x)
diff --git a/pymap3d/eci.py b/pymap3d/eci.py
index 38026fc..f8ae26b 100644
--- a/pymap3d/eci.py
+++ b/pymap3d/eci.py
@@ -41,14 +41,14 @@ def eci2ecef(x: float, y: float, z: float = None,
     z : float
         target z ECEF coordinate
     """
-# %%
+# 
     # FIXME: temporary old API, which was a single N x 3 numpy.ndarray
     if z is None and isinstance(y, (str, datetime)):
         time = y
         z = x[:, 2]
         y = x[:, 1]
         x = x[:, 0]
-# %%
+# 
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     z = np.atleast_1d(z)
@@ -112,14 +112,14 @@ def ecef2eci(x: float, y: float, z: float = None,
     z : float
         target z ECI coordinate
     """
-# %%
+# 
     # FIXME: temporary old API, which was a single N x 3 numpy.ndarray
     if z is None and isinstance(y, (str, datetime)):
         time = y
         z = x[:, 2]
         y = x[:, 1]
         x = x[:, 0]
-# %%
+# 
     x = np.atleast_1d(x)
     y = np.atleast_1d(y)
     z = np.atleast_1d(z)
diff --git a/pymap3d/ellipsoid.py b/pymap3d/ellipsoid.py
index 1ed4939..8fe713f 100644
--- a/pymap3d/ellipsoid.py
+++ b/pymap3d/ellipsoid.py
@@ -1,4 +1,4 @@
-
+from math import sqrt
 
 class Ellipsoid:
     """
@@ -26,19 +26,27 @@ def __init__(self, model='wgs84'):
             self.semimajor_axis = 6378137.
             self.flattening = 1 / 298.2572235630
             self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'wgs72':
             self.semimajor_axis = 6378135.
             self.flattening = 1 / 298.26
             self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'grs80':
             """https://en.wikipedia.org/wiki/GRS_80"""
             self.semimajor_axis = 6378137.
             self.flattening = 1 / 298.257222100882711243
             self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'clarke1866':
             self.semimajor_axis = 6378206.4
             self.semiminor_axis = 6356583.8
             self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'mars':
             """
             https://tharsis.gsfc.nasa.gov/geodesy.html
@@ -46,18 +54,26 @@ def __init__(self, model='wgs84'):
             self.semimajor_axis = 3396900.
             self.semiminor_axis = 3376097.80585952
             self.flattening = 1 / 163.295274386012
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'moon':
             self.semimajor_axis = 1738000.
             self.semiminor_axis = self.semimajor_axis
             self.flattening = 0.
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'venus':
             self.semimajor_axis = 6051000.
             self.semiminor_axis = self.semimajor_axis
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
             self.flattening = 0.
         elif model == 'jupiter':
             self.semimajor_axis = 71492000.
             self.flattening = 1 / 15.415446277169725
             self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'io':
             """
             https://doi.org/10.1006/icar.1998.5987
@@ -65,9 +81,13 @@ def __init__(self, model='wgs84'):
             self.semimajor_axis = 1829.7
             self.flattening = 1 / 131.633093525179
             self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         elif model == 'pluto':
             self.semimajor_axis = 1187000.
             self.semiminor_axis = self.semimajor_axis
             self.flattening = 0.
+            self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) 
+            self.eccentricity = sqrt(2*self.flattening * self.flattening**2) 
         else:
             raise NotImplementedError('{} model not implemented, let us know and we will add it (or make a pull request)'.format(model))
diff --git a/pymap3d/latitude.py b/pymap3d/latitude.py
new file mode 100644
index 0000000..a609367
--- /dev/null
+++ b/pymap3d/latitude.py
@@ -0,0 +1,567 @@
+""" compute auxiliary latitudes and their inverse"""
+import numpy as np
+from .ellipsoid import Ellipsoid
+from numpy import radians, sin, tan, exp, arctan, log, degrees, sqrt, pi
+
+__all__ = ['geodetic2isometric', 'isometric2geodetic', 'geodetic2rectifying', 'rectifying2geodetic', 'geodetic2conformal', 'conformal2geodetic', 'geodetic2parametric', 'parametric2geodetic', 'geodetic2geocentric', 'geocentric2geodetic', 'geodetic2authalic', 'authalic2geodetic']
+
+def geodetic2isometric(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    """
+     computes isometric latitude of a point on an ellipsoid
+
+     Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     isolat : float or numpy.ndarray of float
+         isometric latiude
+
+     Notes
+     -----
+
+     Isometric latitude is an auxiliary latitude proportional to the spacing
+     of parallels of latitude on an ellipsoidal mercator projection.
+
+     Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
+     School of Mathematical and Geospatial Sciences, RMIT University,
+     January 2010
+    """
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = np.deg2rad(lat)
+
+    x = ell.eccentricity * np.sin(lat)
+    y = (1 - x)/(1 + x)
+    z = np.pi/4 + lat/2
+
+#   calculate the isometric latitude
+    isolat = np.log(np.tan(z) * (y**(ell.eccentricity /2)))
+
+    if deg is True:
+        isolat = np.degrees(isolat)
+
+    return isolat
+
+def isometric2geodetic(isolat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from isometric latitude to geodetic latitude
+    
+    Parameters
+     ----------
+
+     isolat : float or numpy.ndarray of float
+         isometric latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     lat : float or numpy.ndarray of float
+         geodetic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        isolat = radians(isolat)
+
+    cnflat = 2 * arctan(exp(isolat)) - (pi/2)
+    lat = conformal2geodetic(cnflat, ell, deg=False)
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
+
+def geodetic2rectifying(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geodetic latitude to rectifying latitude
+    
+    Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     reclat : float or numpy.ndarray of float
+         rectifying latiude
+
+     Notes
+     -----
+
+    Rectifying latitude is an auxiliary latitude used to map an ellipsoid to a 
+    sphere such that correct distances along meridians are preserved.
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    n = ell.thirdflattening
+    f1 = 3*n/2 - 9*n**3/16
+    f2 = 15*n**2/16 - 15*n**4/32
+    f3 = 35*n**3/48
+    f4 = 315*n**4/512
+
+    if deg is True:
+        lat = radians(lat)
+
+    reclat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) + f4*sin(8*lat)
+
+    if deg is True:
+        reclat = degrees(reclat)
+    
+    return reclat
+
+
+def rectifying2geodetic(reclat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from rectifying latitude to geodetic latitude
+    
+    Parameters
+     ----------
+
+     reclat : float or numpy.ndarray of float
+         rectifying latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     lat : float or numpy.ndarray of float
+         geodetic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    n = ell.thirdflattening
+    f1 = 3*n/2 - 27*n**3/32
+    f2 = 21*n**2/16 - 55*n**4/32
+    f3 = 151*n**3/96
+    f4 = 1097*n**4/512
+
+    if deg is True:
+        reclat = radians(reclat)
+
+    lat = reclat + f1*sin(2*reclat) + f2*sin(4*reclat) + f3*sin(6*reclat) + f4*sin(8*reclat)
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
+
+
+def geodetic2conformal(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geodetic latitude to conformal latitude
+    
+    Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     cnflat : float or numpy.ndarray of float
+         conformal latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    e = ell.eccentricity
+    f1 = 1 - e*sin(lat)
+    f2 = 1 + e*sin(lat)
+    f3 = 1 - sin(lat)
+    f4 = 1 + sin(lat)
+
+    #  compute conformal latitudes with correction for points at +90
+    cnflat = np.where(f3 == 0,  pi/2, 2 * arctan(sqrt((f4/f3) * ((f1/f2)**e))) - (pi/2))
+
+    if deg is True:
+        cnflat = degrees(cnflat)
+
+    return cnflat
+
+def conformal2geodetic(cnflat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from conformal latitude to geodetic latitude
+
+    Parameters
+    ----------
+
+    cnflat : float or numpy.ndarray of float
+        conformal latitude
+    ell : Ellipsoid, optional
+        reference ellipsoid (default WGS84)
+    deg : bool, optional
+        degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+
+    lat : float or numpy.ndarray of float
+        geodetic latiude
+
+    Notes
+    -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        cnflat = radians(cnflat)
+
+    e = ell.eccentricity
+    f1 = e**2/2 + 5*e**4/24 + e**6/12 + 13*e**8/360
+    f2 = 7*e**4/48 + 29*e**6/240 + 811*e**8/11520
+    f3 = 7*e**6/120 + 81*e**8/1120
+    f4 = 4279*e**8/161280
+
+    lat = cnflat + f1*sin(2*cnflat) + f2*sin(4*cnflat) + f3*sin(6*cnflat) + f4*sin(8*cnflat)
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
+
+def geodetic2parametric(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geodetic latitude to parametric latitude
+    
+    Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     parlat : float or numpy.ndarray of float
+         parametric latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    epsilon = 1E-10
+
+    lat = np.where(lat == pi/2, pi/2 - epsilon, lat)
+    lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat)
+
+    parlat = arctan(sqrt(1-(ell.eccentricity)**2) * tan(lat))
+
+    if deg is True:
+        parlat = degrees(parlat)
+
+    return parlat
+
+def parametric2geodetic(parlat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from parametric latitude to geodetic latitude
+    
+    Parameters
+     ----------
+
+     parlat : float or numpy.ndarray of float
+         parametric latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     lat : float or numpy.ndarray of float
+         geodetic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        parlat = radians(parlat)
+
+    epsilon = 1E-10
+
+    parlat = np.where(parlat == pi/2, pi/2 - epsilon, parlat)
+    parlat = np.where(parlat == -pi/2, -pi/2 + epsilon, parlat)
+
+    lat = arctan(tan(parlat)/ sqrt(1-(ell.eccentricity)**2))
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
+
+def geodetic2geocentric(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geodetic latitude to geocentric latitude
+    
+    Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     cenlat : float or numpy.ndarray of float
+         geocentric latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    epsilon = 1E-10
+
+    lat = np.where(lat == pi/2, pi/2 - epsilon, lat)
+    lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat)
+
+    cenlat = arctan((1-(ell.eccentricity)**2) * tan(lat))
+
+    if deg is True:
+        cenlat = degrees(cenlat)
+
+    return cenlat
+
+def geocentric2geodetic(cenlat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geocentric latitude to geodetic latitude
+    
+    Parameters
+     ----------
+
+     cenlat : float or numpy.ndarray of float
+         geocentric latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     lat : float or numpy.ndarray of float
+         geodetic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        cenlat = radians(cenlat)
+
+    epsilon = 1E-10
+
+    cenlat = np.where(cenlat == pi/2, pi/2 - epsilon, cenlat)
+    cenlat = np.where(cenlat == -pi/2, -pi/2 + epsilon, cenlat)
+
+    lat = arctan(tan(cenlat)/ (1-(ell.eccentricity)**2))
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
+
+def geodetic2authalic(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from geodetic latitude to authalic latitude
+    
+    Parameters
+     ----------
+
+     lat : float or numpy.ndarray of float
+         geodetic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     autlat : float or numpy.ndarray of float
+         authalic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    e = ell.eccentricity
+    f1 = e**2/3 + 31*e**4/180 + 59*e**6/560
+    f2 = 17*e**4/360 + 61*e**6/1260
+    f3 = 383*e**6/45360
+
+    autlat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) 
+
+    if deg is True:
+        autlat = degrees(autlat)
+
+    return autlat
+
+def authalic2geodetic(autlat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    converts from authalic latitude to geodetic latitude
+    
+    Parameters
+     ----------
+
+     autlat : float or numpy.ndarray of float
+         authalic latitude
+     ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+     deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+     Returns
+     -------
+
+     lat : float or numpy.ndarray of float
+         geodetic latiude
+
+     Notes
+     -----
+
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        autlat = radians(autlat)
+
+    e = ell.eccentricity
+    f1 = e**2/3 + 31*e**4/180 + 517*e**6/5040
+    f2 = 23*e**4/360 + 251*e**6/3780
+    f3 = 761*e**6/45360
+
+    lat = autlat + f1*sin(2*autlat) + f2*sin(4*autlat) + f3*sin(6*autlat)   
+
+    if deg is True:
+        lat = degrees(lat)
+
+    return lat
diff --git a/pymap3d/los.py b/pymap3d/los.py
index ef175fb..bfe47a0 100644
--- a/pymap3d/los.py
+++ b/pymap3d/los.py
@@ -72,13 +72,13 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
 
     magnitude = a**2 * b**2 * w**2 + a**2 * c**2 * v**2 + b**2 * c**2 * u**2
 
-# %%   Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth
+#    Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth
     with np.errstate(invalid='ignore'):
         d = np.where(radical > 0,
                      (value - a * b * c * np.sqrt(radical)) / magnitude,
                      np.nan)
         d[d < 0] = np.nan
-# %% cartesian to ellipsodal
+#  cartesian to ellipsodal
     lat, lon, _ = ecef2geodetic(x + d * u, y + d * v, z + d * w, deg=deg)
 
     return lat, lon, d
diff --git a/pymap3d/lox.py b/pymap3d/lox.py
index 05671af..08a7f04 100644
--- a/pymap3d/lox.py
+++ b/pymap3d/lox.py
@@ -1,149 +1,192 @@
-""" isometric latitude, meridian distance """
 import numpy as np
 from .ellipsoid import Ellipsoid
+from typing import Tuple
+from .latitude import geodetic2isometric, geodetic2rectifying, rectifying2geodetic, geodetic2authalic, authalic2geodetic
+from .rsphere import rsphere_rectifying
+from .rcurve import rcurve_parallel
+from .util import wrapToPi, sph2cart, cart2sph
+from numpy import tan, sign, sin, cos, pi, radians, degrees
 
+__all__ = ['loxodrome_inverse', 'loxodrome_direct', 'meridianarc', 'departure', 'meanm']
 
-def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True):
+def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
+                      ell: Ellipsoid = None, deg: bool = True)  -> Tuple[float, float]:
     """
-     computes isometric latitude of a point on an ellipsoid
+    computes the arc length and azimuth of the loxodrome
+    between two points on the surface of the reference ellipsoid
 
-     Parameters
-     ----------
+    Parameters
+    ----------
 
-     lat : float or numpy.ndarray of float
-         geodetic latitude
-     ell : Ellipsoid, optional
+    lat1 : float or numpy.ndarray of float
+        geodetic latitude of first point
+    lon1 : float or numpy.ndarray of float
+        geodetic longitude of first point
+    lat2 : float or numpy.ndarray of float
+        geodetic latitude of second point
+    lon2 : float or numpy.ndarray of float
+        geodetic longitude of second point
+    ell : Ellipsoid, optional
          reference ellipsoid (default WGS84)
-     deg : bool, optional
+    deg : bool, optional
          degrees input/output  (False: radians in/out)
 
-     Returns
-     -------
-
-     isolat : float or numpy.ndarray of float
-         isometric latiude
+    Results
+    -------
 
-     Notes
-     -----
+    dist : float or numpy.ndarray of float
+        distance along loxodrome
+    a12 : float or numpy.ndarray of float
+        azimuth of loxodrome (degrees/radians)
 
-     Isometric latitude is an auxiliary latitude proportional to the spacing
-     of parallels of latitude on an ellipsoidal mercator projection.
+    Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
+    School of Mathematical and Geospatial Sciences, RMIT University, January 2010
 
-     Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
-     School of Mathematical and Geospatial Sciences, RMIT University,
-     January 2010
+    [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the
+    ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985,
+    pp.223-230.
+    [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S.
+    Geological Survey Professional Paper 1395. Washington, DC: U.S.
+    Government Printing Office, pp.15-16 and pp. 44-45.
+    [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and
+    Cartography, Special Publication No. 251, Coast and Geodetic
+    Survey, U.S. Department of Commerce, Washington, DC: U.S.
+    Government Printing Office, p. 66.
     """
 
+    #   set ellipsoid parameters
     if ell is None:
         ell = Ellipsoid()
 
-    f = ell.flattening  # flattening of ellipsoid
-
     if deg is True:
-        lat = np.deg2rad(lat)
+        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
 
-    e2 = f * (2 - f)  # eccentricity-squared
-    e = np.sqrt(e2)  # eccentricity of ellipsoid
+    # compute isometric latitude of P1 and P2
+    isolat1 = geodetic2isometric(lat1, deg=False, ell=ell)
+    isolat2 = geodetic2isometric(lat2, deg=False, ell=ell)
 
-    x = e * np.sin(lat)
-    y = (1 - x) / (1 + x)
-    z = np.pi / 4 + lat / 2
+    # compute changes in isometric latitude and longitude between points
+    disolat = isolat2 - isolat1
 
-#   calculate the isometric latitude
-    isolat = np.log(np.tan(z) * (y**(e / 2)))
+    dlon = abs(wrapToPi(lon2-lon1))
+
+    # compute azimuth
+    a12 = np.arctan2(dlon, disolat)
+    cosaz = abs(cos(a12))
+
+    # compute distance along loxodromic curve
+    dist = meridianarc(lat2, lat1, deg=False, ell=ell) / cosaz
+
+    # consider degenerate case (directly east/west)
+    epsilon = 1e-10 
+    dist =  np.where(cosaz < epsilon,departure(lon1, lon2, (lat1+lat2)/2, ell, deg=False), dist)
 
     if deg is True:
-        isolat = np.degrees(isolat)
+        a12 = np.degrees(a12) % 360.
 
-    return isolat
+    return dist, a12
 
 
-def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True):
+def loxodrome_direct(lat1: float, lon1: float, rng: float, a12: float,
+            ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
-    computes the meridian distance on an ellipsoid *from the equator* to a latitude.
 
     Parameters
     ----------
 
-    lat : float or numpy.ndarray of float
-        geodetic latitude
+    lat1 : float
+        inital geodetic latitude (degrees)
+    lon1 : float
+        initial geodetic longitude (degrees)
+    rng : float
+        ground distance (meters)
+    a12 : float
+        azimuth (degrees) clockwide from north.
     ell : Ellipsoid, optional
-         reference ellipsoid (default WGS84)
-    deg : bool, optional
-         degrees input/output  (False: radians in/out)
+          reference ellipsoid
 
     Results
     -------
 
-    mdist : float or numpy.ndarray of float
-         meridian distance (degrees/radians)
+    lat2 : float
+        final geodetic latitude (degrees)
+    lon2 : float
+        final geodetic longitude (degrees)
+    a21 : float
+        reverse azimuth (degrees), at final point facing back toward the intial point
 
+"""
+    if ell is None:
+        ell = Ellipsoid()
 
-    Notes
-    -----
+    lat1 = np.atleast_1d(lat1)
+    lon1 = np.atleast_1d(lon1)
+    rng = np.atleast_1d(rng)
+    a12 = np.atleast_1d(a12)
 
-    Formula given Baeschlin, C.F., 1948,
-    "Lehrbuch Der Geodasie", Orell Fussli Verlag, Zurich, pp.47-50.
+    if rng.ndim != 1 or a12.ndim != 1:
+        raise ValueError('Range and azimuth must be scalar or vector')
 
-    Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
-    School of Mathematical and Geospatial Sciences, RMIT University, January 2010
-    """
+    if lat1.size > 1 and rng.size > 1:
+        raise ValueError('LOXODROME_DIRECT: Variable ranges are only allowed for a single point.')
+
+    if rng.size != a12.size and rng.size == 1:
+        rng = np.broadcast_to(rng, a12.size)
 
     if deg is True:
-        lat = np.radians(lat)
+        lat1, lon1, a12 = radians([lat1, lon1, a12])
+    
+    if abs(lat1) > (pi/2):
+        raise ValueError('LOXODROME_DIRECT: Input lat. must be between -90 and 90 deg., inclusive.')
 
-    #   set ellipsoid parameters
-    if ell is None:
-        ell = Ellipsoid()
+    # if range is negative, use inverse bearing
+    a12 = np.where(sign(rng) == -1, (a12 + pi) % (2*pi), a12)
+    rng = abs(rng)
 
-    a = ell.semimajor_axis
-    f = ell.flattening  # flattening of ellipsoid
+    #   compute rectifying sphere latitude and radius
+    reclat = geodetic2rectifying(lat1, ell, deg=False)
 
-    e2 = f * (2 - f)  # eccentricity-squared
+    lat2  = np.zeros_like(lat1)    
+    lon2  = np.zeros_like(lon1)
+    epsilon = 1e-10     # Set tolerance (should set this for specific machine)
 
-    # powers of eccentricity
-    e4 = e2 * e2
-    e6 = e4 * e2
-    e8 = e6 * e2
-    e10 = e8 * e2
+    # correct azimuths at either pole.
+    a12 = np.where(lat1 >= (pi/2) - epsilon, pi, a12)
+    a12 = np.where(lat1 <= (-pi/2) + epsilon, 0,  a12)
 
-    # coefficients of series expansion for meridian distance
-    A = 1 + (3 / 4) * e2 + (45 / 64) * e4 + (175 / 256) * e6 + (11025 / 16384) * e8 + (43659 / 65536) * e10
-    B = (3 / 4) * e2 + (15 / 16) * e4 + (525 / 512) * e6 + (2205 / 2048) * e8 + (72765 / 65536) * e10
-    C = (15 / 64) * e4 + (105 / 256) * e6 + (2205 / 4096) * e8 + (10395 / 16384) * e10
-    D = (35 / 512) * e6 + (315 / 2048) * e8 + (31185 / 131072) * e10
-    E = (315 / 16384) * e8 + (3465 / 65536) * e10
-    F = (693 / 131072) * e10
+    # compute the new points
+    cosaz  = cos(a12)
+    lat2 = reclat + (rng/rsphere_rectifying(ell))*cosaz # compute rectifying latitude
+    lat2 = rectifying2geodetic(lat2,ell,deg=False) #  transform to geodetic latitude
 
-    term1 = A * lat
-    term2 = (B / 2) * np.sin(2 * lat)
-    term3 = (C / 4) * np.sin(4 * lat)
-    term4 = (D / 6) * np.sin(6 * lat)
-    term5 = (E / 8) * np.sin(8 * lat)
-    term6 = (F / 10) * np.sin(10 * lat)
+    # nudge latitudes near either pole
+    lat2 = np.where(lat2 >=  (pi/2) - epsilon, (pi/2) - epsilon, lat2)
+    lat2 = np.where(lat2 <=  (-pi/2) + epsilon, (-pi/2) + epsilon , lat2)
+    lat2 = np.where(abs(cosaz)<=epsilon, lat1, lat2)
 
-    mdist = a * (1 - e2) * (term1 - term2 + term3 - term4 + term5 - term6)
+    newiso = geodetic2isometric(lat2,ell,deg=False)
+    iso    = geodetic2isometric(lat1,ell,deg=False)
+    
+    dlon = np.where(abs(cosaz)<=epsilon, sign(sin(a12)) * rng/rcurve_parallel(lat1, ell, deg=False), tan(a12) * (newiso - iso))
+    lon2 = lon1 + dlon
 
-    return mdist
+    lon2 = wrapToPi(lon2)
+    a21 = (a12 + pi)%(2*pi)
 
+    if deg is True:
+        lat2, lon2, a21 = degrees([lat2, lon2, a21])
 
-def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
-                      ell: Ellipsoid = None, deg: bool = True):
-    """
-    computes the arc length and azimuth of the loxodrome
-    between two points on the surface of the reference ellipsoid
+    return lat2, lon2, a21
+
+def meridianarc(lat1: float, lat2: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    Computes the meridian distance on an ellipsoid between two latitudes.
 
     Parameters
     ----------
 
-    lat1 : float or numpy.ndarray of float
-        geodetic latitude of first point
-    lon1 : float or numpy.ndarray of float
-        geodetic longitude of first point
-    lat2 : float or numpy.ndarray of float
-        geodetic latitude of second point
-    lon2 : float or numpy.ndarray of float
-        geodetic longitude of second point
+    lat1, lat2 : float or numpy.ndarray of float
+        geodetic latitudes
     ell : Ellipsoid, optional
          reference ellipsoid (default WGS84)
     deg : bool, optional
@@ -152,51 +195,49 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
     Results
     -------
 
-    lox_s : float or numpy.ndarray of float
-        distance along loxodrome
-    az12 : float or numpy.ndarray of float
-        azimuth of loxodrome (degrees/radians)
-
-    Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes,
-    School of Mathematical and Geospatial Sciences, RMIT University, January 2010
+    dist : float or numpy.ndarray of float
+         distance (units same as ellipsoid)
+    '''
 
-    [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the
-    ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985,
-    pp.223-230.
-    [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S.
-    Geological Survey Professional Paper 1395. Washington, DC: U.S.
-    Government Printing Office, pp.15-16 and pp. 44-45.
-    [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and
-    Cartography, Special Publication No. 251, Coast and Geodetic
-    Survey, U.S. Department of Commerce, Washington, DC: U.S.
-    Government Printing Office, p. 66.
-    """
+    if deg is True:
+        lat1, lat2 = np.radians([lat1, lat2])
 
     #   set ellipsoid parameters
     if ell is None:
         ell = Ellipsoid()
+    
+    rlat1 = geodetic2rectifying(lat1, ell, deg=False)
+    rlat2 = geodetic2rectifying(lat2, ell, deg=False)
 
-    if deg is True:
-        lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2])
+    return rsphere_rectifying(ell) * abs(rlat2 - rlat1)
 
-    # compute isometric latitude of P1 and P2
-    isolat1 = isometric(lat1, deg=False, ell=ell)
-    isolat2 = isometric(lat2, deg=False, ell=ell)
+def departure(lon1: float, lon2: float, lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    Computes is the distance along a specific parallel between two meridians.
+    '''
 
-    # compute changes in isometric latitude and longitude between points
-    disolat = isolat2 - isolat1
-    dlon = lon2 - lon1
+    if ell is None:
+        ell = Ellipsoid()
 
-    # compute azimuth
-    az12 = np.arctan2(dlon, disolat)
+    if deg is True:
+        lon1, lon2, lat = np.radians([lon1, lon2, lat])
+    
+    return rcurve_parallel(lat, ell, deg=False) * abs(wrapToPi(lon2-lon1))
 
-    # compute distance along loxodromic curve
-    m1 = meridian_dist(lat1, deg=False, ell=ell)
-    m2 = meridian_dist(lat2, deg=False, ell=ell)
-    dm = m2 - m1
-    lox_s = dm / np.cos(az12)
+
+def meanm(lat: float, lon: float, ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float]:
+    '''Computes geographic mean for geographic points on an ellipsoid'''
+    
+    if ell is None:
+        ell = Ellipsoid()
 
     if deg is True:
-        az12 = np.degrees(az12) % 360.
+        lon, lat = radians([lon, lat])
+
+    lat = geodetic2authalic(lat)
+    x,y,z = sph2cart(lon,lat,np.ones_like(lat))
+    latbar, lonbar, _ =cart2sph(np.sum(x),np.sum(y),np.sum(z))
+    latbar = authalic2geodetic(latbar, ell, deg=False)
+    lonbar = wrapToPi(lonbar)
 
-    return lox_s, az12
+    return latbar, lonbar
diff --git a/pymap3d/rcurve.py b/pymap3d/rcurve.py
new file mode 100644
index 0000000..90de10e
--- /dev/null
+++ b/pymap3d/rcurve.py
@@ -0,0 +1,46 @@
+"""compute radii of curvature for an ellipsoid"""
+from .ellipsoid import Ellipsoid
+from numpy import radians, sin, cos, sqrt
+
+__all__ = ['rcurve_parallel', 'rcurve_meridian', 'rcurve_transverse']
+
+def rcurve_parallel(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''
+    computes the radius of the small circle encompassing the globe at the specified latitude
+    '''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    f1 = 1 - ell.eccentricity**2
+    f2 = 1 - (ell.eccentricity * cos(lat))**2
+    return cos(lat) * ell.semimajor_axis * sqrt(f1 / f2)
+
+def rcurve_meridian(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''computes the meridional radius of curvature for the ellipsoid'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    f1 = ell.semimajor_axis * (1 - ell.eccentricity ** 2)
+    f2 = 1 - (ell.eccentricity * sin(lat)) ** 2
+    return f1 / sqrt(f2 ** 3)
+
+def rcurve_transverse(lat, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''computes the radius of the curve formed by a plane 
+    intersecting the ellipsoid at the latitude which is 
+    normal to the surface of the ellipsoid'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    return ell.semimajor_axis / (1 - (ell.eccentricity * sin(lat)) ** 2)
\ No newline at end of file
diff --git a/pymap3d/rsphere.py b/pymap3d/rsphere.py
new file mode 100644
index 0000000..ac59f7d
--- /dev/null
+++ b/pymap3d/rsphere.py
@@ -0,0 +1,107 @@
+""" compute radii of auxiliary spheres"""
+from .ellipsoid import Ellipsoid
+from numpy import radians, sin, cos, log, degrees, sqrt
+from .rcurve import rcurve_meridian, rcurve_transverse
+
+__all__ = ['rsphere_eqavol', 'rsphere_authalic', 'rsphere_rectifying', 'rsphere_euler', 'rsphere_curve', 'rsphere_triaxial', 'rsphere_biaxial']
+
+def rsphere_eqavol(ell: Ellipsoid = None) -> float:
+    '''computes the radius of the sphere with equal volume as the ellipsoid'''
+    if ell is None:
+        ell = Ellipsoid()
+
+    return ell.semimajor_axis * (1 - ell.flattening/3 - ell.flattening**2/9)
+
+def rsphere_authalic(ell: Ellipsoid = None) -> float:
+    '''computes the radius of the sphere with equal surface area as the ellipsoid'''
+    if ell is None:
+        ell = Ellipsoid()
+
+    if ell.eccentricity > 0:
+        f1 = ell.semimajor_axis**2 / 2
+        f2 = (1 - ell.eccentricity**2) / 2*ell.eccentricity
+        f3 = log((1+ell.eccentricity) / (1-ell.eccentricity))
+        return sqrt(f1 * (1 + f2 * f3))
+    else:
+        return ell.semimajor_axis
+
+def rsphere_rectifying(ell: Ellipsoid = None) -> float:
+    '''computes the radius of the sphere with equal meridional distances as the ellipsoid'''
+    if ell is None:
+        ell = Ellipsoid()
+
+    return ((ell.semimajor_axis**(3/2) + ell.semiminor_axis**(3/2))/2)**(2/3)
+    
+
+def rsphere_euler(lat1, lon1, lat2, lon2, ell: Ellipsoid = None, deg: bool = True) -> float:
+    '''computes the Euler radii of curvature at the midpoint of the
+     great circle arc defined by the endpoints (lat1,lon1) and (lat2,lon2)'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is False:
+        lat1, lon1, lat2, lon2 = degrees(lat1, lon1, lat2, lon2)
+
+    latmid = lat1 + (lat2 - lat1)/2  # compute the midpoint
+
+    # compute azimuth
+    _,_,az = vdist(lat1,lon1,lat2,lon2,ell)
+
+    #   compute meridional and transverse radii of curvature
+
+    rho = rcurve_meridian(latmid, ell, deg=True)
+    nu  = rcurve_transverse(latmid, ell, deg=True)
+
+    az = radians(az)
+    den = rho * sin(az)**2 + nu * cos(az)**2;
+
+    #  compute radius of the arc from point 1 to point 2
+    return rho * nu / den
+
+def rsphere_curve(lat, ell: Ellipsoid = None, deg: bool = True, method='mean') -> float:
+    '''computes the arithmetic average of the transverse and meridional
+    radii of curvature at a specified latitude point'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if deg is True:
+        lat = radians(lat)
+
+    rho = rcurve_meridian(lat, ell, deg=False)
+    nu  = rcurve_transverse(lat, ell, deg=False)
+
+    if method=='mean':
+        return (rho + nu) / 2
+    elif method=='norm':
+        return sqrt(rho * nu)
+    else:
+        raise Exception("pymap3d.rsphere.curve: method must be mean or norm")
+
+def rsphere_triaxial(ell: Ellipsoid = None, method='mean') -> float:
+    '''computes triaxial average of the semimajor and semiminor axes of the ellipsoid'''
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if method=='mean':
+        return (2*ell.semimajor_axis + ell.semiminor_axis) / 3
+    elif method=='norm':
+        return (ell.semimajor_axis**2 * ell.semiminor_axis) ** (1/3)
+    else:
+        raise Exception("pymap3d.rsphere.rsphere_triaxial: method must be mean or norm")
+
+
+def rsphere_biaxial(ell: Ellipsoid = None, method='mean') -> float:
+    '''computes biaxial average of the semimajor and semiminor axes of the ellipsoid'''
+    
+    if ell is None:
+        ell = Ellipsoid()
+
+    if method=='mean':
+        return (ell.semimajor_axis + ell.semiminor_axis) / 2
+    elif method=='norm':
+        return sqrt(ell.semimajor_axis * ell.semiminor_axis)
+    else: 
+        raise Exception("pymap3d.rsphere.rsphere_biaxial: method must be mean or norm")
diff --git a/pymap3d/sidereal.py b/pymap3d/sidereal.py
index 95387c5..8df0bcc 100644
--- a/pymap3d/sidereal.py
+++ b/pymap3d/sidereal.py
@@ -44,9 +44,9 @@ def datetime2sidereal(time: datetime,
     usevallado = usevallado or Time is None
     if usevallado:
         jd = juliandate(str2dt(time))
-# %% Greenwich Sidereal time RADIANS
+#  Greenwich Sidereal time RADIANS
         gst = julian2sidereal(jd)
-# %% Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME
+#  Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME
         tsr = gst + lon_radians
     else:
         tsr = Time(time).sidereal_time(kind='apparent',
@@ -122,7 +122,7 @@ def julian2sidereal(Jdate: float) -> float:
 
     tsr = np.empty(jdate.size)
     for i, jd in enumerate(jdate):
-        # %% Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1
+        #  Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1
         tUT1 = (jd - 2451545.0) / 36525.
 
         # Eqn. 3-47 p. 188
diff --git a/pymap3d/util.py b/pymap3d/util.py
new file mode 100644
index 0000000..48a856c
--- /dev/null
+++ b/pymap3d/util.py
@@ -0,0 +1,41 @@
+from numpy import arctan2, hypot, cos, sin, pi
+from typing import Tuple
+
+__all__ = ['cart2pol', 'pol2cart', 'cart2sph', 'sph2cart', 'wrapToPi', 'wrapTo2Pi']
+
+def cart2pol(x: float, y: float) -> Tuple[float, float]:
+    '''Transform Cartesian to polar coordinates'''
+    return arctan2(y, x), hypot(x, y)
+
+def pol2cart(theta: float, rho: float) -> Tuple[float, float]:
+    '''Transform polar to Cartesian coordinates'''
+    return  rho * cos(theta), rho * sin(theta)
+
+def cart2sph(x: float, y: float, z: float) -> Tuple[float, float, float]:
+    '''Transform Cartesian to spherical coordinates'''
+    hxy = hypot(x, y)
+    r = hypot(hxy, z)
+    el = arctan2(z, hxy)
+    az = arctan2(y, x)
+    return az, el, r
+
+def sph2cart(az: float, el: float, r: float) -> Tuple[float, float, float]:
+    '''Transform spherical to Cartesian coordinates'''
+    rcos_theta = r * cos(el)
+    x = rcos_theta * cos(az)
+    y = rcos_theta * sin(az)
+    z = r * sin(el)
+    return x, y, z
+
+def wrapToPi(x: float) -> float:
+    '''wraps angles to [-pi, pi)'''
+    return ((x + pi) % (2 * pi)) - pi
+
+def wrapTo2Pi(x: float) -> float:
+    '''wraps angles to [0, 2*pi]'''
+    return x % (2 * pi)
+
+
+
+
+
diff --git a/pymap3d/vallado.py b/pymap3d/vallado.py
index 04e9673..2b92676 100644
--- a/pymap3d/vallado.py
+++ b/pymap3d/vallado.py
@@ -68,7 +68,7 @@ def azel2radec(az_deg: float, el_deg: float,
     el = radians(el)
     lat = radians(lat)
     lon = radians(lon)
-# %% Vallado "algorithm 28" p 268
+#  Vallado "algorithm 28" p 268
     dec = arcsin(sin(el) * sin(lat) + cos(el) * cos(lat) * cos(az))
 
     lha = arctan2(-(sin(az) * cos(el)) / cos(dec),
@@ -135,11 +135,11 @@ def radec2azel(ra_deg: float, dec_deg: float,
     lon = radians(lon)
 
     lst = datetime2sidereal(time, lon)  # RADIANS
-# %% Eq. 4-11 p. 267 LOCAL HOUR ANGLE
+#  Eq. 4-11 p. 267 LOCAL HOUR ANGLE
     lha = lst - ra
-# %% #Eq. 4-12 p. 267
+#  #Eq. 4-12 p. 267
     el = arcsin(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha))
-# %% combine Eq. 4-13 and 4-14 p. 268
+#  combine Eq. 4-13 and 4-14 p. 268
     az = arctan2(-sin(lha) * cos(dec) / cos(el),
                  (sin(dec) - sin(el) * sin(lat)) / (cos(el) * cos(lat)))
 
diff --git a/pymap3d/vincenty.py b/pymap3d/vincenty.py
index 5ba72c6..ef247c7 100644
--- a/pymap3d/vincenty.py
+++ b/pymap3d/vincenty.py
@@ -86,7 +86,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
     """
     if ell is None:
         ell = Ellipsoid()
-# %% prepare inputs
+#  prepare inputs
     lat1 = np.atleast_1d(Lat1)
     lat2 = np.atleast_1d(Lat2)
     lon1 = np.atleast_1d(Lon1)
@@ -102,21 +102,21 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
         if lat2.size == 1:
             lat2 = np.broadcast_to(lat2, lat1.shape)
             lon2 = np.broadcast_to(lon2, lon1.shape)
-# %% Input check:
+#  Input check:
     if ((abs(lat1) > 90) | (abs(lat2) > 90)).any():
         raise ValueError('Input latitudes must be in [-90, 90] degrees.')
-# %% Supply WGS84 earth ellipsoid axis lengths in meters:
+#  Supply WGS84 earth ellipsoid axis lengths in meters:
     a = ell.semimajor_axis
     b = ell.semiminor_axis
-# %% preserve true input latitudes:
+#  preserve true input latitudes:
     lat1tr = lat1.copy()
     lat2tr = lat2.copy()
-# %% convert inputs in degrees to np.radians:
+#  convert inputs in degrees to np.radians:
     lat1 = np.radians(lat1)
     lon1 = np.radians(lon1)
     lat2 = np.radians(lat2)
     lon2 = np.radians(lon2)
-# %% correct for errors at exact poles by adjusting 0.6 millimeters:
+#  correct for errors at exact poles by adjusting 0.6 millimeters:
     kidx = abs(pi / 2 - abs(lat1)) < 1e-10
     if kidx.any():
         lat1[kidx] = sign(lat1[kidx]) * (pi / 2 - (1e-10))
@@ -203,7 +203,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
 
     dist_m = (b * A * (sigma - deltasigma))
 
-# %% From point #1 to point #2
+#  From point #1 to point #2
     # correct sign of lambda for azimuth calcs:
     lamb = abs(lamb)
     kidx = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0
@@ -213,12 +213,12 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
     a12 = arctan2(numer, denom)
     kidx = a12 < 0
     a12[kidx] = a12[kidx] + 2 * pi
-    # %% from poles
+    #  from poles
     a12[lat1tr <= -90] = 0
     a12[lat1tr >= 90] = pi
     az = np.degrees(a12)
 
-# %% From point #2 to point #1
+#  From point #2 to point #1
     # correct sign of lambda for azimuth calcs:
     lamb = abs(lamb)
     kidx = sign(sin(lon1 - lon2)) * sign(sin(lamb)) < 0
@@ -228,7 +228,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
     a21 = arctan2(numer, denom)
     kidx = a21 < 0
     a21[kidx] = a21[kidx] + 2 * pi
-    # %% backwards from poles:
+    #  backwards from poles:
     a21[lat2tr >= 90] = pi
     a21[lat2tr <= -90] = 0.
     backaz = np.degrees(a21)