diff --git a/pymap3d/__init__.py b/pymap3d/__init__.py
index 1483e7a..6781172 100644
--- a/pymap3d/__init__.py
+++ b/pymap3d/__init__.py
@@ -18,7 +18,17 @@
 
 deg : bool = True means degrees. False = radians.
 
-Most functions accept NumPy arrays of any shape.
+Most functions accept NumPy arrays of any shape, as well as compatible data types
+including AstroPy, Pandas and Xarray that have Numpy-like data properties.
+For clarity, we omit all these types in the docs, and just specify the scalar type.
+
+Other languages
+---------------
+
+Companion packages exist for:
+
+* Matlab / GNU Octave: [Matmap3D](https://github.com/scivision/matmap3d)
+* Fortran: [Maptran3D](https://github.com/scivision/maptran3d)
 """
 from .timeconv import str2dt  # noqa: F401
 from .azelradec import radec2azel, azel2radec  # noqa: F401
diff --git a/pymap3d/aer.py b/pymap3d/aer.py
index c0bf472..b005f25 100644
--- a/pymap3d/aer.py
+++ b/pymap3d/aer.py
@@ -18,28 +18,31 @@ def ecef2aer(x: float, y: float, z: float,
     Parameters
     ----------
 
-    x : float
-    y : float
-    z : float
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
     lat0 : float
-           Observer geodetic latitude
+        Observer geodetic latitude
     lon0 : float
-           Observer geodetic longitude
+        Observer geodetic longitude
     h0 : float
          observer altitude above geodetic ellipsoid (meters)
     ell : Ellipsoid, optional
-          reference ellipsoid
+        reference ellipsoid
     deg : bool, optional
-          degrees input/output  (False: radians in/out)
+        degrees input/output  (False: radians in/out)
 
     Returns
     -------
-    az : float
-         azimuth (degrees)  [0,360)
-    el : float
-         elevation (degrees/radians)  [0,90]
-    srange : float
-         slant range [meters] [0,Infinity)
+    az : float or numpy.ndarray of float
+         azimuth to target
+    el : float or numpy.ndarray of float
+         elevation to target
+    srange : float or numpy.ndarray of float
+         slant range [meters]
     """
     xEast, yNorth, zUp = ecef2enu(x, y, z, lat0, lon0, h0, ell, deg=deg)
 
@@ -56,16 +59,16 @@ def geodetic2aer(lat: float, lon: float, h: float,
     Parameters
     ----------
 
-    lat : float
-           target geodetic latitude
-    lon : float
-           target geodetic longitude
-    h : float
-         target altitude above geodetic ellipsoid (meters)
+    lat : float or numpy.ndarray of float
+        target geodetic latitude
+    lon : float or numpy.ndarray of float
+        target geodetic longitude
+    h : float or numpy.ndarray of float
+        target altitude above geodetic ellipsoid (meters)
     lat0 : float
-           Observer geodetic latitude
+        Observer geodetic latitude
     lon0 : float
-           Observer geodetic longitude
+        Observer geodetic longitude
     h0 : float
          observer altitude above geodetic ellipsoid (meters)
     ell : Ellipsoid, optional
@@ -75,12 +78,12 @@ def geodetic2aer(lat: float, lon: float, h: float,
 
     Returns
     -------
-    az : float
-         azimuth (degrees)  [0,360)
-    el : float
-         elevation (degrees/radians)  [0,90]
-    srange : float
-         slant range [meters] [0,Infinity)
+    az : float or numpy.ndarray of float
+         azimuth
+    el : float or numpy.ndarray of float
+         elevation
+    srange : float or numpy.ndarray of float
+         slant range [meters]
     """
     e, n, u = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
 
@@ -95,15 +98,14 @@ def aer2geodetic(az: float, el: float, srange: float,
     gives geodetic coordinates of a point with az, el, range
     from an observer at lat0, lon0, h0
 
-
     Parameters
     ----------
-    az : float
-         azimuth (degrees)  [0,360)
-    el : float
-         elevation (degrees/radians)  [0,90]
-    srange : float
-         slant range [meters] [0,Infinity)
+    az : float or numpy.ndarray of float
+         azimuth to target
+    el : float or numpy.ndarray of float
+         elevation to target
+    srange : float or numpy.ndarray of float
+         slant range [meters]
     lat0 : float
            Observer geodetic latitude
     lon0 : float
@@ -120,11 +122,11 @@ def aer2geodetic(az: float, el: float, srange: float,
 
     In reference ellipsoid system:
 
-    lat : float
+    lat : float or numpy.ndarray of float
           geodetic latitude
-    lon : float
+    lon : float or numpy.ndarray of float
           geodetic longitude
-    alt : float
+    alt : float or numpy.ndarray of float
           altitude above ellipsoid  (meters)
     """
     x, y, z = aer2ecef(az, el, srange, lat0, lon0, h0, ell=ell, deg=deg)
@@ -157,11 +159,11 @@ def eci2aer(eci: Tuple[float, float, float],
     Returns
     -------
     az : float
-         azimuth (degrees)  [0,360)
+         azimuth to target
     el : float
-         elevation (degrees/radians)  [0,90]
+         elevation to target
     srange : float
-         slant range [meters] [0,Infinity)
+         slant range [meters]
     """
     ecef = np.atleast_2d(eci2ecef(eci, t, useastropy))
 
@@ -177,12 +179,12 @@ def aer2eci(az: float, el: float, srange: float,
 
     Parameters
     ----------
-    az : float
-         azimuth (degrees)  [0,360)
-    el : float
-         elevation (degrees/radians)  [0,90]
-    srange : float
-         slant range [meters] [0,Infinity)
+    az : float or numpy.ndarray of float
+         azimuth to target
+    el : float or numpy.ndarray of float
+         elevation to target
+    srange : float or numpy.ndarray of float
+         slant range [meters]
     lat0 : float
            Observer geodetic latitude
     lon0 : float
@@ -201,9 +203,12 @@ def aer2eci(az: float, el: float, srange: float,
 
     Earth Centered Inertial x,y,z
 
-    x : float
-    y : float
-    z : float
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
     """
     x, y, z = aer2ecef(az, el, srange, lat0, lon0, h0, ell, deg)
 
@@ -218,12 +223,12 @@ def aer2ecef(az: float, el: float, srange: float,
 
     Parameters
     ----------
-    az : float
-         azimuth (degrees)  [0,360)
-    el : float
-         elevation (degrees/radians)  [0,90]
-    srange : float
-         slant range [meters] [0,Infinity)
+    az : float or numpy.ndarray of float
+         azimuth to target
+    el : float or numpy.ndarray of float
+         elevation to target
+    srange : float or numpy.ndarray of float
+         slant range [meters]
     lat0 : float
            Observer geodetic latitude
     lon0 : float
@@ -240,9 +245,13 @@ def aer2ecef(az: float, el: float, srange: float,
 
     ECEF (Earth centered, Earth fixed)  x,y,z
 
-    x : float
-    y : float
-    z : float
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
+
 
     Notes
     ------
diff --git a/pymap3d/azelradec.py b/pymap3d/azelradec.py
index bee9b02..6c1b8ab 100644
--- a/pymap3d/azelradec.py
+++ b/pymap3d/azelradec.py
@@ -22,9 +22,9 @@ def azel2radec(az_deg: float, el_deg: float,
 
     Parameters
     ----------
-    az_deg : float
+    az_deg : float or numpy.ndarray of float
          azimuth [degrees clockwize from North]
-    el_deg : float
+    el_deg : float or numpy.ndarray of float
              elevation [degrees above horizon (neglecting aberration)]
     lat_deg : float
               observer latitude [-90, 90]
@@ -37,9 +37,9 @@ def azel2radec(az_deg: float, el_deg: float,
 
     Returns
     -------
-    ra_deg : float
+    ra_deg : float or numpy.ndarray of float
          ecliptic right ascension (degress)
-    dec_deg : float
+    dec_deg : float or numpy.ndarray of float
          ecliptic declination (degrees)
     """
 
@@ -64,9 +64,9 @@ def radec2azel(ra_deg: float, dec_deg: float,
 
     Parameters
     ----------
-    ra_deg : float
+    ra_deg : float or numpy.ndarray of float
          ecliptic right ascension (degress)
-    dec_deg : float
+    dec_deg : float or numpy.ndarray of float
          ecliptic declination (degrees)
     lat_deg : float
               observer latitude [-90, 90]
@@ -79,9 +79,9 @@ def radec2azel(ra_deg: float, dec_deg: float,
 
     Returns
     -------
-    az_deg : float
+    az_deg : float or numpy.ndarray of float
              azimuth [degrees clockwize from North]
-    el_deg : float
+    el_deg : float or numpy.ndarray of float
              elevation [degrees above horizon (neglecting aberration)]
     """
 
diff --git a/pymap3d/ecef.py b/pymap3d/ecef.py
index 6c42c8a..aa07545 100644
--- a/pymap3d/ecef.py
+++ b/pymap3d/ecef.py
@@ -24,9 +24,14 @@ class Ellipsoid:
     https://nssdc.gsfc.nasa.gov/planetary/factsheet/index.html
     """
 
-    def __init__(self, model: str = 'wgs84') -> None:
+    def __init__(self, model: str = 'wgs84'):
         """
         feel free to suggest additional ellipsoids
+
+        Parameters
+        ----------
+        model : str
+                name of ellipsoid
         """
         if model == 'wgs84':
             """https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84"""
@@ -63,7 +68,23 @@ def __init__(self, model: str = 'wgs84') -> None:
 
 
 def get_radius_normal(lat_radians: float, ell: Ellipsoid = None) -> float:
-    """ Compute normal radius of planetary body"""
+    """
+    Compute normal radius of planetary body
+
+    Parameters
+    ----------
+
+    lat_radians : float
+        latitude in radians
+    ell : Ellipsoid, optional
+        reference ellipsoid
+
+    Returns
+    -------
+
+    radius : float
+        normal radius (meters)
+    """
     if ell is None:
         ell = Ellipsoid()
 
@@ -81,11 +102,11 @@ def geodetic2ecef(lat: float, lon: float, alt: float,
     Parameters
     ----------
 
-    lat : float
+    lat : float or numpy.ndarray of float
            target geodetic latitude
-    lon : float
+    lon : float or numpy.ndarray of float
            target geodetic longitude
-    h : float
+    h : float or numpy.ndarray of float
          target altitude above geodetic ellipsoid (meters)
     ell : Ellipsoid, optional
           reference ellipsoid
@@ -98,12 +119,12 @@ def geodetic2ecef(lat: float, lon: float, alt: float,
 
     ECEF (Earth centered, Earth fixed)  x,y,z
 
-    x : float
-        target x ECEF coordinate
-    y : float
-        target y ECEF coordinate
-    z : float
-        target z ECEF coordinate
+    x : float or numpy.ndarray of float
+        target x ECEF coordinate (meters)
+    y : float or numpy.ndarray of float
+        target y ECEF coordinate (meters)
+    z : float or numpy.ndarray of float
+        target z ECEF coordinate (meters)
     """
     if ell is None:
         ell = Ellipsoid()
@@ -138,12 +159,12 @@ def ecef2geodetic(x: float, y: float, z: float,
 
     Parameters
     ----------
-    x : float
-        target x ECEF coordinate
-    y : float
-        target y ECEF coordinate
-    z : float
-        target z ECEF coordinate
+    x : float or numpy.ndarray of float
+        target x ECEF coordinate (meters)
+    y : float or numpy.ndarray of float
+        target y ECEF coordinate (meters)
+    z : float or numpy.ndarray of float
+        target z ECEF coordinate (meters)
     ell : Ellipsoid, optional
           reference ellipsoid
     deg : bool, optional
@@ -151,11 +172,11 @@ def ecef2geodetic(x: float, y: float, z: float,
 
     Returns
     -------
-    lat : float
+    lat : float or numpy.ndarray of float
            target geodetic latitude
-    lon : float
+    lon : float or numpy.ndarray of float
            target geodetic longitude
-    h : float
+    h : float or numpy.ndarray of float
          target altitude above geodetic ellipsoid (meters)
 
     based on:
@@ -213,12 +234,12 @@ def ecef2enuv(u: float, v: float, w: float,
 
     Parameters
     ----------
-    u : float
-        target x ECEF coordinate
-    v : float
-        target y ECEF coordinate
-    w : float
-        target z ECEF coordinate
+    u : float or numpy.ndarray of float
+        target x ECEF coordinate (meters)
+    v : float or numpy.ndarray of float
+        target y ECEF coordinate (meters)
+    w : float or numpy.ndarray of float
+        target z ECEF coordinate (meters)
     lat0 : float
            Observer geodetic latitude
     lon0 : float
@@ -230,12 +251,12 @@ def ecef2enuv(u: float, v: float, w: float,
 
     Returns
     -------
-    uEast : float
-        target east ENU coordinate
-    vNorth : float
-        target north ENU coordinate
-    wUp : float
-        target up ENU coordinate
+    uEast : float or numpy.ndarray of float
+        target east ENU coordinate (meters)
+    vNorth : float or numpy.ndarray of float
+        target north ENU coordinate (meters)
+    wUp : float or numpy.ndarray of float
+        target up ENU coordinate (meters)
 
     """
     if deg:
@@ -258,12 +279,12 @@ def ecef2enu(x: float, y: float, z: float,
 
     Parameters
     ----------
-    x : float
-        target x ECEF coordinate
-    y : float
-        target y ECEF coordinate
-    z : float
-        target z ECEF coordinate
+    x : float or numpy.ndarray of float
+        target x ECEF coordinate (meters)
+    y : float or numpy.ndarray of float
+        target y ECEF coordinate (meters)
+    z : float or numpy.ndarray of float
+        target z ECEF coordinate (meters)
     lat0 : float
            Observer geodetic latitude
     lon0 : float
@@ -277,12 +298,12 @@ def ecef2enu(x: float, y: float, z: float,
 
     Returns
     -------
-    East : float
-        target east ENU coordinate
-    North : float
-        target north ENU coordinate
-    Up : float
-        target up ENU coordinate
+    East : float or numpy.ndarray of float
+        target east ENU coordinate (meters)
+    North : float or numpy.ndarray of float
+        target north ENU coordinate (meters)
+    Up : float or numpy.ndarray of float
+        target up ENU coordinate (meters)
 
     """
     x0, y0, z0 = geodetic2ecef(lat0, lon0, h0, ell, deg=deg)
@@ -292,6 +313,25 @@ def ecef2enu(x: float, y: float, z: float,
 
 def enu2uvw(east: float, north: float, up: float,
             lat0: float, lon0: float, deg: bool = True) -> Tuple[float, float, float]:
+    """
+    Parameters
+    ----------
+
+    e1 : float or numpy.ndarray of float
+        target east ENU coordinate (meters)
+    n1 : float or numpy.ndarray of float
+        target north ENU coordinate (meters)
+    u1 : float or numpy.ndarray of float
+        target up ENU coordinate (meters)
+
+    Results
+    -------
+
+    u : float or numpy.ndarray of float
+    v : float or numpy.ndarray of float
+    w : float or numpy.ndarray of float
+    """
+
     if deg:
         lat0 = radians(lat0)
         lon0 = radians(lon0)
@@ -307,6 +347,25 @@ def enu2uvw(east: float, north: float, up: float,
 
 def uvw2enu(u: float, v: float, w: float,
             lat0: float, lon0: float, deg: bool = True) -> Tuple[float, float, float]:
+    """
+    Parameters
+    ----------
+
+    u : float or numpy.ndarray of float
+    v : float or numpy.ndarray of float
+    w : float or numpy.ndarray of float
+
+
+    Results
+    -------
+
+    East : float or numpy.ndarray of float
+        target east ENU coordinate (meters)
+    North : float or numpy.ndarray of float
+        target north ENU coordinate (meters)
+    Up : float or numpy.ndarray of float
+        target up ENU coordinate (meters)
+    """
     if deg:
         lat0 = radians(lat0)
         lon0 = radians(lon0)
@@ -326,7 +385,7 @@ def eci2geodetic(eci: np.ndarray, t: datetime,
 
     Parameters
     ----------
-    eci : tuple
+    eci : tuple of float
           [meters] Nx3 target ECI location (x,y,z)
     t : datetime.datetime, float
           length N vector of datetime OR greenwich sidereal time angle [radians].
@@ -360,16 +419,35 @@ def enu2ecef(e1: float, n1: float, u1: float,
     """
     ENU to ECEF
 
-    inputs:
-     e1, n1, u1 (meters)   east, north, up
-     observer: lat0, lon0, h0 (degrees/radians,degrees/radians, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
+    Parameters
+    ----------
 
+    e1 : float or numpy.ndarray of float
+        target east ENU coordinate (meters)
+    n1 : float or numpy.ndarray of float
+        target north ENU coordinate (meters)
+    u1 : float or numpy.ndarray of float
+        target up ENU coordinate (meters)
+    lat0 : float
+           Observer geodetic latitude
+    lon0 : float
+           Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
 
-    output
-    ------
-    x,y,z  [meters] target ECEF location                         [0,Infinity)
+
+    Results
+    -------
+    x : float or numpy.ndarray of float
+        target x ECEF coordinate (meters)
+    y : float or numpy.ndarray of float
+        target y ECEF coordinate (meters)
+    z : float or numpy.ndarray of float
+        target z ECEF coordinate (meters)
     """
     x0, y0, z0 = geodetic2ecef(lat0, lon0, h0, ell, deg=deg)
     dx, dy, dz = enu2uvw(e1, n1, u1, lat0, lon0, deg=deg)
diff --git a/pymap3d/eci.py b/pymap3d/eci.py
index 7fa0631..1806daf 100644
--- a/pymap3d/eci.py
+++ b/pymap3d/eci.py
@@ -13,16 +13,25 @@ def eci2ecef(eci: np.ndarray,
              time: datetime,
              useastropy: bool = True) -> np.ndarray:
     """
-     Observer => Point
-
-    input
-    -----
-    eci [meters] Nx3 target ECI location (x,y,z)                    [0,Infinity)
-    t  time (datetime.datetime)   time of obsevation (UTC)
-
-    output
-    ------
-    x,y,z  [meters] target ECEF location                             [0,Infinity)
+    Observer => Point  ECI  =>  ECEF
+
+    Parameters
+    ----------
+    eci : tuple of float
+        Nx3 target ECI location (x,y,z) [meters]
+    time : datetime.datetime
+        time of obsevation (UTC)
+    useastropy : bool, optional
+        use AstroPy for conversion
+
+    Results
+    -------
+    x : float
+        target x ECEF coordinate
+    y : float
+        target y ECEF coordinate
+    z : float
+        target z ECEF coordinate
     """
     useastropy = useastropy and Time
 
@@ -53,17 +62,26 @@ def ecef2eci(ecef: np.ndarray,
              time: datetime,
              useastropy: bool = True) -> np.ndarray:
     """
-    Point => Point
+    Point => Point   ECEF => ECI
 
     input
     -----
-    ecef:  Nx3  x,y,z  (meters)
-    time:  datetime.datetime
-
-
-    output
-    ------
-    eci  x,y,z (meters)
+    ecef : tuple of float
+        Nx3 target ECEF location (x,y,z) [meters]
+    time : datetime.datetime
+        time of observation
+    useastropy : bool, optional
+        use AstroPy for conversion
+
+
+    Results
+    -------
+    x : float
+        target x ECI coordinate
+    y : float
+        target y ECI coordinate
+    z : float
+        target z ECI coordinate
     """
     useastropy = useastropy and Time
 
@@ -90,12 +108,19 @@ def ecef2eci(ecef: np.ndarray,
 
 
 def _rottrip(ang: np.ndarray) -> np.ndarray:
+    """
+    transformation matrix
+
+    Parameters
+    ----------
+
+    ang : N x 3 numpy.ndarray
+        angle to transform (radians)
+    """
     ang = ang.squeeze()
     if ang.size > 1:
         raise ValueError('only one angle allowed at a time')
-    """ported from:
-    https://github.com/dinkelk/astrodynamics/blob/master/rot3.m
-    """
+
     return np.array([[np.cos(ang), np.sin(ang), 0],
                      [-np.sin(ang), np.cos(ang), 0],
                      [0, 0, 1]])
diff --git a/pymap3d/enu.py b/pymap3d/enu.py
index caf0b89..8d67216 100644
--- a/pymap3d/enu.py
+++ b/pymap3d/enu.py
@@ -17,15 +17,27 @@ def enu2aer(e: np.ndarray, n: np.ndarray, u: np.ndarray, deg: bool = True) -> Tu
     """
     ENU to Azimuth, Elevation, Range
 
-    ## Inputs
-
-    * e,n,u  [meters]  East, north, up                                [0,Infinity)
-    * deg    degrees input/output  (False: radians in/out)
-
-    ## Outputs
+    Parameters
+    ----------
+
+    e : float or np.ndarray of float
+        ENU East coordinate (meters)
+    n : float or np.ndarray of float
+        ENU North coordinate (meters)
+    u : float or np.ndarray of float
+        ENU Up coordinate (meters)
+    deg : bool, optional
+        degrees input/output  (False: radians in/out)
+
+    Results
+    -------
 
-    * azimuth, elevation (degrees/radians)                             [0,360),[0,90]
-    * slant range [meters]                                             [0,Infinity)
+    azimuth : float or np.ndarray of float
+        azimuth to rarget
+    elevation : float or np.ndarray of float
+        elevation to target
+    srange : float or np.ndarray of float
+        slant range [meters]
     """
     # 1 millimeter precision for singularity
 
@@ -52,16 +64,27 @@ def enu2aer(e: np.ndarray, n: np.ndarray, u: np.ndarray, deg: bool = True) -> Tu
 
 def aer2enu(az: float, el: float, srange: float, deg: bool = True) -> Tuple[float, float, float]:
     """
-    input:
-    ------
-    azimuth, elevation (degrees/radians)                             [0,360),[0,90]
-    slant range [meters]                                             [0,Infinity)
-    deg    degrees input/output  (False: radians in/out)
-
-    output:
-    -------
-    e,n,u   East, North, Up [m]
-
+    Azimuth, Elevation, Slant range to target to East, north, Up
+
+    Parameters
+    ----------
+    azimuth : float or np.ndarray of float
+            azimuth clockwise from north (degrees)
+    elevation : float or np.ndarray of float
+        elevation angle above horizon, neglecting aberattions (degrees)
+    srange : float or np.ndarray of float
+        slant range [meters]
+    deg : bool, optional
+        degrees input/output  (False: radians in/out)
+
+    Returns
+    --------
+    e : float or np.ndarray of float
+        East ENU coordinate (meters)
+    n : float or np.ndarray of float
+        North ENU coordinate (meters)
+    u : float or np.ndarray of float
+        Up ENU coordinate (meters)
     """
     if deg:
         el = radians(el)
@@ -80,19 +103,36 @@ def enu2geodetic(e: float, n: float, u: float,
                  lat0: float, lon0: float, h0: float,
                  ell=None, deg: bool = True) -> Tuple[float, float, float]:
     """
-
-    input
-    -----
-    e,n,u   East, North, Up [m]
-    Observer: lat0, lon0, h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-
-    output:
+    East, North, Up to target to geodetic coordinates
+
+    Parameters
+    ----------
+    e : float or np.ndarray of float
+        East ENU coordinate (meters)
+    n : float or np.ndarray of float
+        North ENU coordinate (meters)
+    u : float or np.ndarray of float
+        Up ENU coordinate (meters)
+    lat0 : float
+           Observer geodetic latitude
+    lon0 : float
+           Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+
+    Results
     -------
-    target: lat,lon, h  (degrees/radians,degrees/radians, meters)
-
+    lat : float or np.ndarray of float
+          geodetic latitude
+    lon : float or np.ndarray of float
+          geodetic longitude
+    alt : float or np.ndarray of float
+          altitude above ellipsoid  (meters)
     """
 
     x, y, z = enu2ecef(e, n, u, lat0, lon0, h0, ell, deg=deg)
@@ -104,17 +144,34 @@ def geodetic2enu(lat: float, lon: float, h: float,
                  lat0: float, lon0: float, h0: float,
                  ell=None, deg: bool = True) -> Tuple[float, float, float]:
     """
-    input
-    -----
-    target: lat,lon, h
-    Observer: lat0, lon0, h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-
-    output:
+    Parameters
+    ----------
+    lat : float or np.ndarray of float
+          target geodetic latitude
+    lon : float or np.ndarray of float
+          target geodetic longitude
+    h : float or np.ndarray of float
+          target altitude above ellipsoid  (meters)
+    lat0 : float
+           Observer geodetic latitude
+    lon0 : float
+           Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+
+    Results
     -------
-    e,n,u   East, North, Up [m]
+    e : float or np.ndarray of float
+        East ENU
+    n : float or np.ndarray of float
+        North ENU
+    u : float or np.ndarray of float
+        Up ENU
     """
     x1, y1, z1 = geodetic2ecef(lat, lon, h, ell, deg=deg)
     x2, y2, z2 = geodetic2ecef(lat0, lon0, h0, ell, deg=deg)
diff --git a/pymap3d/haversine.py b/pymap3d/haversine.py
index 366d420..718e7e1 100644
--- a/pymap3d/haversine.py
+++ b/pymap3d/haversine.py
@@ -2,7 +2,8 @@
 """
 Michael Hirsch
 
-Note: decimal points on constants made 0 difference in %timeit execution time
+Note:
+decimal points on constants made 0 difference in `%timeit` execution time
 
 The Meeus algorithm is about 9.5% faster than Astropy/Vicenty on my PC,
 and gives virtually identical result
@@ -25,8 +26,28 @@
 def anglesep_meeus(lon0: float, lat0: float,
                    lon1: float, lat1: float, deg: bool = True) -> float:
     """
-    inputs: DEGREES  (right ascension, declination Meeus p. 109)
+    Parameters
+    ----------
 
+    lon0 : float or numpy.ndarray of float
+        longitude of first point
+    lat0 : float or numpy.ndarray of float
+        latitude of first point
+    lon1 : float or numpy.ndarray of float
+        longitude of second point
+    lat1 : float or numpy.ndarray of float
+        latitude of second point
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+
+    sep_rad : float or numpy.ndarray of float
+        angular separation
+
+
+    Meeus p. 109
 
     from "Astronomical Algorithms" by Jean Meeus Ch. 16 p. 111 (16.5)
     gives angular distance in degrees between two rightAscension,Declination
@@ -34,8 +55,6 @@ def anglesep_meeus(lon0: float, lat0: float,
 
     Meeus haversine method is stable all the way to exactly 0 deg.
 
-    assumes degrees input, degrees output
-
     either the arrays must be the same size, or one of them must be a scalar
     """
 
@@ -57,7 +76,25 @@ def anglesep_meeus(lon0: float, lat0: float,
 def anglesep(lon0: float, lat0: float,
              lon1: float, lat1: float, deg: bool = True) -> float:
     """
-    inputs: DEGREES
+    Parameters
+    ----------
+
+    lon0 : float
+        longitude of first point
+    lat0 : float
+        latitude of first point
+    lon1 : float
+        longitude of second point
+    lat1 : float
+        latitude of second point
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+
+    sep_rad : float or numpy.ndarray of float
+        angular separation
 
     For reference, this is from astropy astropy/coordinates/angle_utilities.py
     Angular separation between two points on a sphere.
@@ -81,9 +118,21 @@ def anglesep(lon0: float, lat0: float,
 
 def haversine(theta: float) -> float:
     """
-    Compute haversine of angle theta (radians)
+    Compute haversine
+
+    Parameters
+    ----------
+
+    theta : float
+        angle (radians)
+
+    Results
+    -------
+
+    htheta : float
+        haversine of `theta`
 
-    http://en.wikipedia.org/wiki/Haversine
+    https://en.wikipedia.org/wiki/Haversine
     Meeus p. 111
     """
     return (1 - cos(theta)) / 2.
diff --git a/pymap3d/los.py b/pymap3d/los.py
index b2f7302..8380451 100644
--- a/pymap3d/los.py
+++ b/pymap3d/los.py
@@ -9,38 +9,48 @@
 
 
 def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
-                   ell=Ellipsoid(), deg: bool = True) -> Tuple[float, float, float]:
+                   ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
     Calculates line-of-sight intersection with Earth (or other ellipsoid) surface from above surface / orbit
 
-    ## Inputs
-
-    lat0, lon0
-    : latitude and longitude of starting point
-
-    h0
-    : altitude of starting point in meters
-
-    az
-    : azimuth angle of line-of-sight, clockwise from North
-
-    tilt
-    : tilt angle of line-of-sight with respect to local vertical (nadir = 0)
-
-    ## Outputs
-
-    lat, lon
-    : latitude and longitude where the line-of-sight intersects with the Earth ellipsoid
-
-    d
-    : slant range in meters from the starting point to the intersect point
+    Parameters
+    ----------
+
+    lat0 : float
+           observer geodetic latitude
+    lon0 : float
+           observer geodetic longitude
+    h0 : float
+        observer altitude (meters)  Must be non-negative since this function doesn't consider terrain
+    az : float or numpy.ndarray of float
+        azimuth angle of line-of-sight, clockwise from North
+    tilt : float or numpy.ndarray of float
+        tilt angle of line-of-sight with respect to local vertical (nadir = 0)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
+    -------
+
+    lat0 : float or numpy.ndarray of float
+           geodetic latitude where the line-of-sight intersects with the Earth ellipsoid
+    lon0 : float or numpy.ndarray of float
+           geodetic longitude where the line-of-sight intersects with the Earth ellipsoid
+    d : float or numpy.ndarray of float
+        slant range (meters) from starting point to intersect point
 
     Values will be NaN if the line of sight does not intersect.
 
     Algorithm based on https://medium.com/@stephenhartzell/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6 Stephen Hartzell
     """
+
     if (np.asarray(h0) < 0).any():
-        raise ValueError('Altitude   [0, Infinity)')
+        raise ValueError('Intersection calculation requires altitude  [0, Infinity)')
+
+    if ell is None:
+        ell = Ellipsoid()
 
     tilt = np.asarray(tilt)
 
@@ -61,13 +71,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
     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 dedd949..995a10c 100644
--- a/pymap3d/lox.py
+++ b/pymap3d/lox.py
@@ -7,23 +7,24 @@ def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True):
     """
      computes isometric latitude of a point on an ellipsoid
 
-     ## Inputs
+     Parameters
+     ----------
 
-     lat
-     : latitude (degrees/radians)
+     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)
 
-     ell
-     :reference ellipsoid
+     Returns
+     -------
 
-     deg
-     : degrees input/output  (False: radians in/out)
+     isolat : float or numpy.ndarray of float
+         isometric latiude
 
-     ## Outputs
-
-     isolat
-     : isometric latiude (degrees/radians)
-
-     ## Notes
+     Notes
+     -----
 
      Isometric latitude is an auxiliary latitude proportional to the spacing
      of parallels of latitude on an ellipsoidal mercator projection.
@@ -61,23 +62,25 @@ def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True):
     """
     computes the meridian distance on an ellipsoid *from the equator* to a latitude.
 
-    ## Inputs
-
-    lat
-    : latitude (degrees/radians)
+    Parameters
+    ----------
 
-    ell
-    : reference ellipsoid
+    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)
 
-    deg
-    : degrees input/output  (False: radians in/out)
+    Results
+    -------
 
-    ## Outputs
+    mdist : float or numpy.ndarray of float
+         meridian distance (degrees/radians)
 
-    mdist
-    : meridian distance (degrees/radians)
 
-    ## Notes
+    Notes
+    -----
 
     Formula given Baeschlin, C.F., 1948,
     "Lehrbuch Der Geodasie", Orell Fussli Verlag, Zurich, pp.47-50.
@@ -130,30 +133,29 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float,
     computes the arc length and azimuth of the loxodrome
     between two points on the surface of the reference ellipsoid
 
-    ## Inputs
-
-    lat1
-    : GEODETIC latitude of first point (degrees/radians)
-
-    lon1
-    : longitude of first point (degrees/radians)
-
-    lat2, lon2
-    : second point (degrees/radians)
-
-    ell
-    : reference ellipsoid
-
-    deg
-    : degrees input/output  (False: radians in/out)
-
-    ## Outputs
-
-    lox_s
-    : distance along loxodrome
-
-    az12
-    : azimuth of loxodrome (degrees/radians)
+    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
+    ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+    deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+    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
diff --git a/pymap3d/ned.py b/pymap3d/ned.py
index 8ddc992..e8b733e 100644
--- a/pymap3d/ned.py
+++ b/pymap3d/ned.py
@@ -1,59 +1,105 @@
 """ Transforms involving NED North East Down """
 from .enu import geodetic2enu, aer2enu, enu2aer
-from .ecef import ecef2geodetic, ecef2enuv, ecef2enu, enu2ecef
+from .ecef import ecef2geodetic, ecef2enuv, ecef2enu, enu2ecef, Ellipsoid
 from typing import Tuple
 
 
-def aer2ned(az: float, elev: float, slantRange: float, deg: bool = True) -> Tuple[float, float, float]:
+def aer2ned(az: float, elev: float, slantRange: float,
+            deg: bool = True) -> Tuple[float, float, float]:
     """
-    input
-    -----
-    azimuth, elevation (degrees/radians)                             [0,360),[0,90]
-    slant range [meters]                                             [0,Infinity)
-    deg    degrees input/output  (False: radians in/out)
+    converts azimuth, elevation, range to target from observer to North, East, Down
 
-    output:
+    Parameters
+    -----------
+
+    az : float or numpy.ndarray of float
+         azimuth
+    elev : float or numpy.ndarray of float
+         elevation
+    slantRange : float or numpy.ndarray of float
+         slant range [meters]
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
     -------
-    n,e,d  North,east, down [m]
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
     """
     e, n, u = aer2enu(az, elev, slantRange, deg=deg)
 
     return n, e, -u
 
 
-def ned2aer(n: float, e: float, d: float, deg: bool = True) -> Tuple[float, float, float]:
+def ned2aer(n: float, e: float, d: float,
+            deg: bool = True) -> Tuple[float, float, float]:
     """
-    Observer => Point
+    converts North, East, Down to azimuth, elevation, range
+
+    Parameters
+    ----------
 
-    input
-    -----
-    n,e,d  [meters]  North,east, down                                [0,Infinity)
-    deg    degrees input/output  (False: radians in/out)
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
+    deg : bool, optional
+        degrees input/output  (False: radians in/out)
 
-    output: AER
-    ------
-    azimuth, elevation (degrees/radians)                             [0,360),[0,90]
-    slant range [meters]                                             [0,Infinity)
+    Results
+    -------
+
+    az : float or numpy.ndarray of float
+         azimuth
+    elev : float or numpy.ndarray of float
+         elevation
+    slantRange : float or numpy.ndarray of float
+         slant range [meters]
     """
     return enu2aer(e, n, -d, deg=deg)
 
 
 def ned2geodetic(n: float, e: float, d: float,
                  lat0: float, lon0: float, h0: float,
-                 ell=None, deg: bool = True) -> Tuple[float, float, float]:
+                 ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
-
-    input
-    -----
-    n,e,d   North, east, down (meters)
-    Observer: lat0, lon0, h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-
-    output:
+    Converts North, East, Down to target latitude, longitude, altitude
+
+    Parameters
+    ----------
+
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
+    lat0 : float
+        Observer geodetic latitude
+    lon0 : float
+        Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
     -------
-    target: lat,lon, h  (degrees/radians,degrees/radians, meters)
+
+    lat : float
+        target geodetic latitude
+    lon : float
+        target geodetic longitude
+    h : float
+        target altitude above geodetic ellipsoid (meters)
 
     """
     x, y, z = enu2ecef(e, n, -d, lat0, lon0, h0, ell, deg=deg)
@@ -63,36 +109,78 @@ def ned2geodetic(n: float, e: float, d: float,
 
 def ned2ecef(n: float, e: float, d: float,
              lat0: float, lon0: float, h0: float,
-             ell=None, deg: bool = True) -> Tuple[float, float, float]:
+             ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
-    Observer => Point
-
-    input
-    -----
-    n,e,d  [meters]  North,east, down                                [0,Infinity)
-    Observer: lat0, lon0, h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-    output:
-    ------
-    ECEF  x,y,z  (meters)
+    North, East, Down to target ECEF coordinates
+
+    Parameters
+    ----------
+
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
+    lat0 : float
+        Observer geodetic latitude
+    lon0 : float
+        Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
+    -------
+
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
     """
     return enu2ecef(e, n, -d, lat0, lon0, h0, ell, deg=deg)
 
 
-def ecef2ned(x, y, z, lat0, lon0, h0, ell=None, deg=True) -> Tuple[float, float, float]:
+def ecef2ned(x: float, y: float, z: float,
+             lat0: float, lon0: float, h0: float,
+             ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
-    input
-    -----
-    x,y,z  [meters] target ECEF location                             [0,Infinity)
-    Observer: lat0, lon0, h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-    output:
+    Convert ECEF x,y,z to North, East, Down
+
+    Parameters
+    ----------
+
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
+    lat0 : float
+        Observer geodetic latitude
+    lon0 : float
+        Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
     -------
-    n,e,d  North,east, down [m]
+
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
 
     """
     e, n, u = ecef2enu(x, y, z, lat0, lon0, h0, ell, deg=deg)
@@ -100,31 +188,81 @@ def ecef2ned(x, y, z, lat0, lon0, h0, ell=None, deg=True) -> Tuple[float, float,
     return n, e, -u
 
 
-def geodetic2ned(lat, lon, h, lat0, lon0, h0, ell=None, deg=True) -> Tuple[float, float, float]:
+def geodetic2ned(lat: float, lon: float, h: float,
+                 lat0: float, lon0: float, h0: float,
+                 ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]:
     """
-    input
-    -----
-    target: lat,lon  (degrees/radians)
-    h (altitude, meters)
-    Observer: lat0, lon0 (degrees/radians)
-    h0 (altitude, meters)
-    ell    reference ellipsoid
-    deg    degrees input/output  (False: radians in/out)
-
-
-    output:
+    convert latitude, longitude, altitude of target to North, East, Down from observer
+
+    Parameters
+    ----------
+
+    lat : float
+        target geodetic latitude
+    lon : float
+        target geodetic longitude
+    h : float
+        target altitude above geodetic ellipsoid (meters)
+    lat0 : float
+        Observer geodetic latitude
+    lon0 : float
+        Observer geodetic longitude
+    h0 : float
+         observer altitude above geodetic ellipsoid (meters)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+
+    Results
     -------
-    n,e,d  North,east, down [m]
+
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
     """
     e, n, u = geodetic2enu(lat, lon, h, lat0, lon0, h0, ell, deg=deg)
 
     return n, e, -u
 
 
-def ecef2nedv(u, v, w, lat0, lon0, deg=True) -> Tuple[float, float, float]:
+def ecef2nedv(x: float, y: float, z: float,
+              lat0: float, lon0: float,
+              deg: bool = True) -> Tuple[float, float, float]:
     """
     for VECTOR between two points
+
+    Parameters
+    ----------
+    x : float or numpy.ndarray of float
+        ECEF x coordinate (meters)
+    y : float or numpy.ndarray of float
+        ECEF y coordinate (meters)
+    z : float or numpy.ndarray of float
+        ECEF z coordinate (meters)
+    lat0 : float
+        Observer geodetic latitude
+    lon0 : float
+        Observer geodetic longitude
+    deg : bool, optional
+          degrees input/output  (False: radians in/out)
+
+    Results
+    -------
+
+    (Vector)
+
+    n : float or numpy.ndarray of float
+        North NED coordinate (meters)
+    e : float or numpy.ndarray of float
+        East NED coordinate (meters)
+    d : float or numpy.ndarray of float
+        Down NED coordinate (meters)
     """
-    e, n, u = ecef2enuv(u, v, w, lat0, lon0, deg=deg)
+    e, n, u = ecef2enuv(x, y, z, lat0, lon0, deg=deg)
 
     return n, e, -u
diff --git a/pymap3d/sidereal.py b/pymap3d/sidereal.py
index 333385b..95387c5 100644
--- a/pymap3d/sidereal.py
+++ b/pymap3d/sidereal.py
@@ -25,14 +25,21 @@ def datetime2sidereal(time: datetime,
     """
     Convert ``datetime`` to sidereal time
 
-    :algorithm: D. Vallado Fundamentals of Astrodynamics and Applications
+    from D. Vallado "Fundamentals of Astrodynamics and Applications"
 
 
-    time
-        Python datetime
+    time : datetime.datetime
+        time to convert
+    lon_radians : float
+        longitude (radians)
+    usevallado : bool, optional
+        use vallado instead of AstroPy (default is Vallado)
 
-    lon
-        longitude in RADIANS
+    Results
+    -------
+
+    tsr : float
+        Sidereal time
     """
     usevallado = usevallado or Time is None
     if usevallado:
@@ -54,6 +61,18 @@ def juliandate(time: datetime) -> float:
 
     from D.Vallado Fundamentals of Astrodynamics and Applications p.187
      and J. Meeus Astronomical Algorithms 1991 Eqn. 7.1 pg. 61
+
+    Parameters
+    ----------
+
+    time : datetime.datetime
+        time to convert
+
+    Results
+    -------
+
+    jd : float
+        Julian date
     """
 
     times = np.atleast_1d(time)
@@ -84,11 +103,18 @@ def julian2sidereal(Jdate: float) -> float:
 
     D. Vallado Ed. 4
 
-    input:
+    Parameters
+    ----------
 
-    juliandate
+    Jdate: float
         Julian centuries from J2000.0
 
+    Results
+    -------
+
+    tsr : float
+        Sidereal time
+
     """
 
     jdate = np.atleast_1d(Jdate)
diff --git a/pymap3d/timeconv.py b/pymap3d/timeconv.py
index bb14129..f605e26 100644
--- a/pymap3d/timeconv.py
+++ b/pymap3d/timeconv.py
@@ -11,9 +11,16 @@ def str2dt(time: datetime) -> np.ndarray:
     """
     Converts times in string or list of strings to datetime(s)
 
-    ## Output
+    Parameters
+    ----------
+
+    time : str or datetime.datetime or numpy.datetime64
+
+    Results
+    -------
+
+    t : datetime.datetime
 
-    datetime
     """
     if isinstance(time, datetime):
         return time
diff --git a/pymap3d/vallado.py b/pymap3d/vallado.py
index f661c94..5122dc7 100644
--- a/pymap3d/vallado.py
+++ b/pymap3d/vallado.py
@@ -20,32 +20,34 @@ def azel2radec(az_deg: float, el_deg: float,
                lat_deg: float, lon_deg: float,
                time: datetime) -> Tuple[float, float]:
     """
-    `azel2radec` converts azimuth, elevation to right ascension, declination
+    converts azimuth, elevation to right ascension, declination
 
-    ## Inputs
+    Parameters
+    ----------
 
-    az_deg
-    : Numpy ndarray of azimuth to point [degrees]
+    az_deg : float or numpy.ndarray of float
+        azimuth (clockwise) to point [degrees]
 
-    el_deg
-    : Numpy ndarray of elevation to point [degrees]
+    el_deg : float or numpy.ndarray of float
+        elevation above horizon to point [degrees]
 
-    lat_deg
-    : scalar observer WGS84 latitude [degrees]
+    lat_deg : float
+        observer WGS84 latitude [degrees]
 
-    lon_deg
-    : scalar observer WGS84 longitude [degrees]
+    lon_deg : float
+        observer WGS84 longitude [degrees]
 
-    time
-    : time of observation
+    time : datetime.datetime
+        time of observation
 
-    ## Outputs
+    Results
+    -------
 
-    ra_deg
-    : Numpy ndarray of right ascension values [degrees]
+    ra_deg : float or numpy.ndarray of float
+        right ascension to target [degrees]
 
-    dec_deg
-    : Numpy ndarray of declination values [degrees]
+    dec_deg : float or numpy.ndarray of float
+        declination of target [degrees]
 
     from D.Vallado Fundamentals of Astrodynamics and Applications
     p.258-259
@@ -84,32 +86,34 @@ def radec2azel(ra_deg: float, dec_deg: float,
                lat_deg: float, lon_deg: float,
                time: datetime) -> Tuple[float, float]:
     """
-    `radec2azel` converts right ascension, declination to azimuth, elevation
+    converts right ascension, declination to azimuth, elevation
 
-    ## Inputs
+    Parameters
+    ----------
 
-    ra_deg
-    : Numpy ndarray of right ascension values [degrees]
+    ra_deg : float or numpy.ndarray of float
+        right ascension to target [degrees]
 
-    dec_deg
-    : Numpy ndarray of declination values [degrees]
+    dec_deg : float or numpy.ndarray of float
+        declination to target [degrees]
 
-    lat_deg
-    : scalar observer WGS84 latitude [degrees]
+    lat_deg : float
+        observer WGS84 latitude [degrees]
 
-    lon_deg
-    : scalar observer WGS84 longitude [degrees]
+    lon_deg : float
+        observer WGS84 longitude [degrees]
 
-    time
-    : time of observation
+    time : datetime.datetime
+        time of observation
 
-    ## Outputs
+    Results
+    -------
 
-    az_deg
-    : Numpy ndarray of azimuth to point [degrees]
+    az_deg : float or numpy.ndarray of float
+        azimuth clockwise from north to point [degrees]
 
-    el_deg
-    : Numpy ndarray of elevation to point [degrees]
+    el_deg : float or numpy.ndarray of float
+        elevation above horizon to point [degrees]
 
 
     from D. Vallado "Fundamentals of Astrodynamics and Applications "
diff --git a/pymap3d/vincenty.py b/pymap3d/vincenty.py
index 46ffe69..2b3d77e 100644
--- a/pymap3d/vincenty.py
+++ b/pymap3d/vincenty.py
@@ -13,37 +13,40 @@
 
 def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = None) -> Tuple[float, float, float]:
     """
-    Using the WGS-84 Earth ellipsoid, compute the distance between two points
+    Using the reference ellipsoid, compute the distance between two points
     within a few millimeters of accuracy, compute forward azimuth,
     and compute backward azimuth, all using a vectorized version of
+    Vincenty's algorithm:
 
-    Vincenty's algorithm::
-
-        s = vdist(lat1,lon1,lat2,lon2)
-        s,a12 = vdist(lat1,lon1,lat2,lon2)
-        s,a12,a21 = vdist(lat1,lon1,lat2,lon2)
-
-    ## Inputs
-
-    s
-   : distance in meters (inputs may be scalars, vectors, or matrices)
+    ```python
+    s,a12,a21 = vdist(lat1, lon1, lat2, lon2, ell)
+    ```
 
-    a12
-    : azimuth in degrees from first point to second point (forward)
+    Parameters
+    ----------
 
-    a21
-    : azimuth in degrees from second point to first point (backward)
+    Lat1 : float or numpy.ndarray of float
+        Geodetic latitude of first point (degrees)
+    Lon1 : float or numpy.ndarray of float
+        Geodetic longitude of first point (degrees)
+    Lat2 : float or numpy.ndarray of float
+        Geodetic latitude of second point (degrees)
+    Lon2 : float or numpy.ndarray of float
+        Geodetic longitude of second point (degrees)
+    ell : Ellipsoid, optional
+          reference ellipsoid
 
-    (Azimuths are in degrees clockwise from north.)
+    Results
+    -------
 
-    lat1
-    : GEODETIC latitude of first point (degrees)
+    s : float or numpy.ndarray of float
+        distance (meters)
+    a12 : float or numpy.ndarray of float
+        azimuth (degrees) clockwise from first point to second point (forward)
+    a21 : float or numpy.ndarray of float
+        azimuth (degrees) clockwise from second point to first point (backward)
 
-    lon1
-    : longitude of first point (degrees)
 
-    lat2, lon2
-    : second point(s) (degrees)
 
     Original algorithm source:
     T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid
@@ -52,7 +55,8 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
     Available at: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
 
     Notes:
-     1. lat1,lon1,lat2,lon2 can be any (identical) size/shape. Outputs  will have the same size and shape.
+
+     1. lat1,lon1,lat2,lon2 can be any (identical) size/shape. Outputs will have the same size and shape.
      2. Error correcting code, convergence failure traps, antipodal corrections, polar error corrections, WGS84 ellipsoid parameters, testing, and comments: Michael Kleder, 2004.
      3. Azimuth implementation (including quadrant abiguity resolution) and code vectorization, Michael Kleder, Sep 2005.
      4. Vectorization is convergence sensitive; that is, quantities which have already converged to within tolerance are not recomputed during subsequent iterations (while other quantities are still converging).
@@ -61,7 +65,8 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
      7. Distance failures for points exactly at the poles are eliminated by moving the points by 0.6 millimeters.
      8. The Vincenty distance algorithm was transcribed verbatim by Peter Cederholm, August 12, 2003. It was modified and translated to English by Michael Kleder. Mr. Cederholm's website is http://www.plan.aau.dk/~pce/
      9. Distances agree with the Mapping Toolbox, version 2.2 (R14SP3) with a max relative difference of about 5e-9, except when the two points are nearly antipodal, and except when one point is near the equator and the two longitudes are nearly 180 degrees apart. This function (vdist) is more accurate in such cases.
-        For example, note this difference (as of this writing)::
+        For example, note this difference (as of this writing):
+
         ```python
         vdist(0.2,305,15,125)
         ```
@@ -73,11 +78,11 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N
         ```
 
         > 0
+
      10. Azimuths FROM the north pole (either forward starting at the north pole or backward when ending at the north pole) are set to 180 degrees by convention.
          Azimuths FROM the south pole are set to 0 degrees by convention.
      11. Azimuths agree with the Mapping Toolbox, version 2.2 (R14SP3) to within about a hundred-thousandth of a degree, except when traversing to or from a pole, where the convention for this function is described in (10), and except in the cases noted above in (9).
      12. No warranties; use at your own risk.
-
     """
     if ell is None:
         ell = Ellipsoid()
@@ -237,45 +242,39 @@ def vreckon(Lat1: float, Lon1: float, Rng: float, Azim: float,
     This is the Vincenty "forward" solution.
 
     Computes points at a specified azimuth and range in an ellipsoidal earth.
-    Using the WGS-84 Earth ellipsoid, travel a given distance along a given azimuth starting at a given initial point,
+    Using the reference ellipsoid, travel a given distance along a given azimuth starting at a given initial point,
     and return the endpoint within a few millimeters of accuracy, using Vincenty's algorithm.
 
-    ## Usage
-
     ```python
-    lat2,lon2 = vreckon(lat1, lon1, rng, azim)
+    lat2, lon2 = vreckon(lat1, lon1, rng, azim)
     ```
 
-    Transmits ellipsoid definition (either as [a,b] or [a,f]) as fifth argument ELLIPSOID
-
-
-    ## Inputs
-
-    lat1
-    : inital latitude (degrees)
+    Parameters
+    ----------
 
-    lon1
-    : initial longitude (degrees)
+    Lat1 : float
+        inital geodetic latitude (degrees)
+    Lon1 : float
+        initial geodetic longitude (degrees)
+    Rng : float
+        ground distance (meters)
+    Azim : float
+        intial azimuth (degrees) clockwide from north.
+    ell : Ellipsoid, optional
+          reference ellipsoid
 
-    rng
-    : distance (meters). Scalar or a vector. Latter case computes a series of circles (or arc circles, see azim) centered on X,Y (which are scalars)
+    Results
+    -------
 
-    azim
-    : intial azimuth (degrees). "azim" is a scalar or vector
+    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
 
-    ellipsoid
-    : two-element ellipsoid vector. Either [a b] or [a f] If omitted, defaults to WGS-84
 
-
-    ## Outputs
-
-    lat2, lon2
-    : second point (degrees)
-
-    a21
-    : reverse azimuth (degrees), at final point facing back toward the intial point
-
-    :Original algorithm: T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations", Survey Review, vol. 23, no. 176, April 1975, pp 88-93. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
+    Original algorithm: T. Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations", Survey Review, vol. 23, no. 176, April 1975, pp 88-93. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
 
     Notes:
 
@@ -296,7 +295,6 @@ def vreckon(Lat1: float, Lon1: float, Rng: float, Azim: float,
 
     Added ellipsoid and vectorized whenever possible. Also, lon2 is always converted to the [-180 180] interval.
     Joaquim Luis
-
     """
 
     lat1 = np.atleast_1d(Lat1)
@@ -418,29 +416,36 @@ def vreckon(Lat1: float, Lon1: float, Rng: float, Azim: float,
 def track2(lat1: float, lon1: float, lat2: float, lon2: float,
            ell: Ellipsoid = None, npts: int = 100, deg: bool = True):
     """
-     computes great circle tracks starting at the point lat1, lon1 and ending at lat2, lon2
-
-     input
-     -----
-     lat1
-     GEODETIC latitude of first point (degrees/radians)
-
-     lon1
-     longitude of first point (degrees/radians)
-
-     lat2, lon2
-     second point (degrees/radians)
-
-     ell    reference ellipsoid
-     npts   number of points (default is 100)
-     deg    degrees input/output  (False: radians in/out)
-
-     output
-     ------
-     lats, lons  latitudes and longitudes of points along great circle track
-
-     Based on code posted to the GMT mailing list in Dec 1999 by Jim Levens and by Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
-     """
+    computes great circle tracks starting at the point lat1, lon1 and ending at lat2, lon2
+
+    Parameters
+    ----------
+
+    Lat1 : float or numpy.ndarray of float
+        Geodetic latitude of first point (degrees)
+    Lon1 : float or numpy.ndarray of float
+        Geodetic longitude of first point (degrees)
+    Lat2 : float or numpy.ndarray of float
+        Geodetic latitude of second point (degrees)
+    Lon2 : float or numpy.ndarray of float
+        Geodetic longitude of second point (degrees)
+    ell : Ellipsoid, optional
+          reference ellipsoid
+    npts : int, optional
+        number of points (default is 100)
+    deg : bool, optional
+        degrees input/output  (False: radians in/out)
+
+    Results
+    -------
+
+    lats : numpy.ndarray of float
+        latitudes of points along track
+    lons : numpy.ndarray of float
+        longitudes of points along track
+
+    Based on code posted to the GMT mailing list in Dec 1999 by Jim Levens and by Jeff Whitaker <jeffrey.s.whitaker@noaa.gov>
+    """
 
     if ell is None:
         ell = Ellipsoid()