put ellipsoid in own file

explicit ellipsoid properties

enable limited fallback for Python 2.7-3.4

py26

py26 ci

.travis.yml

@@ -17,6 +17,14 @@ matrix:
     python: 3.8-dev
     install: pip install -e .[tests]
     script: pytest tests/test_math.py
+  - os: linux
+    python: 3.4
+    install: pip install -e .[tests]
+    script: pytest tests/test_py27.py
+  - os: linux
+    python: 2.6
+    install: pip install -e .[tests]
+    script: pytest tests/test_py27.py
   - os: linux
     python: 3.7
     install: pip install -e .[tests,full,cov]

README.md

@@ -26,15 +26,15 @@ Includes some relevant
 
 ## Prerequisites
 
-Python ≥ 3.4
+Python ≥ 2.6  (full features require Python ≥ 3.5)
 
 or
 
 PyPy3
 
 References to Numpy and AstroPy are *optional*, algorithms from Vallado and Meeus are used if AstroPy is not present.
 PyMap3D is regularly tested with Python ≥ 3.5.
-Limited Python 3.4 support is available for systems using MicroPython or other cases where a current Python version isn't available.
+Limited Python 2.6, 2.7 and 3.4 support is available for systems using MicroPython or other cases where a current Python version isn't available.
 
 
 ## Install

pymap3d/__init__.py

@@ -30,7 +30,18 @@
 * Matlab / GNU Octave: [Matmap3D](https://github.com/scivision/matmap3d)
 * Fortran: [Maptran3D](https://github.com/scivision/maptran3d)
 """
-try:
+import sys
+
+if sys.version_info >= (3, 5):
+    try:
+        import numpy
+        import dateutil  # noqa: F401
+    except ImportError:
+        numpy = None
+else:
+    numpy = None
+
+if numpy is not None:
     from .timeconv import str2dt  # noqa: F401
     from .azelradec import radec2azel, azel2radec  # noqa: F401
     from .eci import eci2ecef, ecef2eci  # noqa: F401
@@ -41,5 +52,5 @@
     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
-except ImportError:  # pure Python only
+else:  # pure Python only
     from .math import *  # type: ignore  # noqa: F401, F403

pymap3d/ecef.py

@@ -9,89 +9,9 @@
     tau = 2 * pi
 
 from .eci import eci2ecef
+from .ellipsoid import Ellipsoid
 
-__all__ = ['Ellipsoid', 'geodetic2ecef', 'ecef2geodetic', 'ecef2enuv', 'ecef2enu', 'enu2uvw', 'uvw2enu', 'eci2geodetic', 'enu2ecef']
-
-
-class Ellipsoid:
-    """
-    generate reference ellipsoid parameters
-
-    https://en.wikibooks.org/wiki/PROJ.4#Spheroid
-
-    https://tharsis.gsfc.nasa.gov/geodesy.html
-
-    https://nssdc.gsfc.nasa.gov/planetary/factsheet/index.html
-    """
-
-    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"""
-            self.a = 6378137.  # semi-major axis [m]
-            self.f = 1 / 298.2572235630  # flattening
-            self.b = self.a * (1 - self.f)  # semi-minor axis
-        elif model == 'grs80':
-            """https://en.wikipedia.org/wiki/GRS_80"""
-            self.a = 6378137.  # semi-major axis [m]
-            self.f = 1 / 298.257222100882711243  # flattening
-            self.b = self.a * (1 - self.f)  # semi-minor axis
-        elif model == 'clrk66':  # Clarke 1866
-            self.a = 6378206.4  # semi-major axis [m]
-            self.b = 6356583.8  # semi-minor axis
-            self.f = -(self.b / self.a - 1)
-        elif model == 'mars':  # https://tharsis.gsfc.nasa.gov/geodesy.html
-            self.a = 3396900
-            self.b = 3376097.80585952
-            self.f = 1 / 163.295274386012
-        elif model == 'moon':
-            self.a = 1738000.
-            self.b = self.a
-            self.f = 0.
-        elif model == 'venus':
-            self.a = 6051000.
-            self.b = self.a
-            self.f = 0.
-        elif model == 'pluto':
-            self.a = 1187000.
-            self.b = self.a
-            self.f = 0.
-        else:
-            raise NotImplementedError('{} model not implemented, let us know and we will add it (or make a pull request)'.format(model))
-
-
-def get_radius_normal(lat_radians: float, ell: Ellipsoid = None) -> float:
-    """
-    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()
-
-    a = ell.a
-    b = ell.b
-
-    return a**2 / sqrt(a**2 * cos(lat_radians)**2 + b**2 * sin(lat_radians)**2)
+__all__ = ['geodetic2ecef', 'ecef2geodetic', 'ecef2enuv', 'ecef2enu', 'enu2uvw', 'uvw2enu', 'eci2geodetic', 'enu2ecef']
 
 
 def geodetic2ecef(lat: float, lon: float, alt: float,
@@ -139,12 +59,12 @@ def geodetic2ecef(lat: float, lon: float, alt: float,
             raise ValueError('-90 <= lat <= 90')
 
     # radius of curvature of the prime vertical section
-    N = get_radius_normal(lat, ell)
+    N = ell.semimajor_axis**2 / sqrt(ell.semimajor_axis**2 * cos(lat)**2 + ell.semiminor_axis**2 * sin(lat)**2)
     # Compute cartesian (geocentric) coordinates given  (curvilinear) geodetic
     # coordinates.
     x = (N + alt) * cos(lat) * cos(lon)
     y = (N + alt) * cos(lat) * sin(lon)
-    z = (N * (ell.b / ell.a)**2 + alt) * sin(lat)
+    z = (N * (ell.semiminor_axis / ell.semimajor_axis)**2 + alt) * sin(lat)
 
     return x, y, z
 
@@ -189,7 +109,7 @@ def ecef2geodetic(x: float, y: float, z: float,
 
     r = sqrt(x**2 + y**2 + z**2)
 
-    E = sqrt(ell.a**2 - ell.b**2)
+    E = sqrt(ell.semimajor_axis**2 - ell.semiminor_axis**2)
 
     # eqn. 4a
     u = sqrt(0.5 * (r**2 - E**2) + 0.5 * sqrt((r**2 - E**2)**2 + 4 * E**2 * z**2))
@@ -203,21 +123,21 @@ def ecef2geodetic(x: float, y: float, z: float,
         Beta = arctan(huE / u * z / hypot(x, y))
 
     # eqn. 13
-    eps = ((ell.b * u - ell.a * huE + E**2) * sin(Beta)) / (ell.a * huE * 1 / cos(Beta) - E**2 * cos(Beta))
+    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
-    lat = arctan(ell.a / ell.b * tan(Beta))
+    lat = arctan(ell.semimajor_axis / ell.semiminor_axis * tan(Beta))
 
     lon = arctan2(y, x)
 
     # eqn. 7
-    alt = hypot(z - ell.b * sin(Beta),
-                Q - ell.a * cos(Beta))
+    alt = hypot(z - ell.semiminor_axis * sin(Beta),
+                Q - ell.semimajor_axis * cos(Beta))
 
     # inside ellipsoid?
     with np.errstate(invalid='ignore'):
-        inside = x**2 / ell.a**2 + y**2 / ell.a**2 + z**2 / ell.b**2 < 1
+        inside = x**2 / ell.semimajor_axis**2 + y**2 / ell.semimajor_axis**2 + z**2 / ell.semiminor_axis**2 < 1
     if isinstance(inside, np.ndarray):
         alt[inside] = -alt[inside]
     elif inside:

pymap3d/ellipsoid.py

@@ -0,0 +1,73 @@
+
+
+class Ellipsoid:
+    """
+    generate reference ellipsoid parameters
+
+    https://en.wikibooks.org/wiki/PROJ.4#Spheroid
+
+    https://nssdc.gsfc.nasa.gov/planetary/factsheet/index.html
+
+    as everywhere else in this program, distance units are METERS
+    """
+
+    # def __init__(self, model: str = 'wgs84'):
+    def __init__(self, model='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"""
+            self.semimajor_axis = 6378137.
+            self.flattening = 1 / 298.2572235630
+            self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+        elif model == 'wgs72':
+            self.semimajor_axis = 6378135.
+            self.flattening = 1 / 298.26
+            self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+        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)
+        elif model == 'clarke1866':
+            self.semimajor_axis = 6378206.4
+            self.semiminor_axis = 6356583.8
+            self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1)
+        elif model == 'mars':
+            """
+            https://tharsis.gsfc.nasa.gov/geodesy.html
+            """
+            self.semimajor_axis = 3396900.
+            self.semiminor_axis = 3376097.80585952
+            self.flattening = 1 / 163.295274386012
+        elif model == 'moon':
+            self.semimajor_axis = 1738000.
+            self.semiminor_axis = self.semimajor_axis
+            self.flattening = 0.
+        elif model == 'venus':
+            self.semimajor_axis = 6051000.
+            self.semiminor_axis = self.semimajor_axis
+            self.flattening = 0.
+        elif model == 'jupiter':
+            self.semimajor_axis = 71492000.
+            self.flattening = 1 / 15.415446277169725
+            self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1)
+        elif model == 'io':
+            """
+            https://doi.org/10.1006/icar.1998.5987
+            """
+            self.semimajor_axis = 1829.7
+            self.flattening = 1 / 131.633093525179
+            self.semiminor_axis = self.semimajor_axis * (1 - self.flattening)
+        elif model == 'pluto':
+            self.semimajor_axis = 1187000.
+            self.semiminor_axis = self.semimajor_axis
+            self.flattening = 0.
+        else:
+            raise NotImplementedError('{} model not implemented, let us know and we will add it (or make a pull request)'.format(model))

pymap3d/los.py

@@ -54,9 +54,9 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float,
 
     tilt = np.asarray(tilt)
 
-    a = ell.a
-    b = ell.a
-    c = ell.b
+    a = ell.semimajor_axis
+    b = ell.semimajor_axis
+    c = ell.semiminor_axis
 
     el = tilt - 90. if deg else tilt - pi / 2
 

pymap3d/lox.py

@@ -37,7 +37,7 @@ def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True):
     if ell is None:
         ell = Ellipsoid()
 
-    f = ell.f  # flattening of ellipsoid
+    f = ell.flattening  # flattening of ellipsoid
 
     if deg is True:
         lat = np.deg2rad(lat)
@@ -96,8 +96,8 @@ def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True):
     if ell is None:
         ell = Ellipsoid()
 
-    a = ell.a
-    f = ell.f  # flattening of ellipsoid
+    a = ell.semimajor_axis
+    f = ell.flattening  # flattening of ellipsoid
 
     e2 = f * (2 - f)  # eccentricity-squared
 

pymap3d/math/aer.py

@@ -1,14 +1,13 @@
 """ transforms involving AER: azimuth, elevation, slant range"""
 
-from .ecef import ecef2enu, geodetic2ecef, ecef2geodetic, enu2uvw, Ellipsoid
+from .ecef import ecef2enu, geodetic2ecef, ecef2geodetic, enu2uvw
 from .enu import geodetic2enu, aer2enu, enu2aer
 
 __all__ = ['aer2ecef', 'ecef2aer', 'geodetic2aer', 'aer2geodetic']
 
 
-def ecef2aer(x: float, y: float, z: float,
-             lat0: float, lon0: float, h0: float,
-             ell: Ellipsoid = None, deg: bool = True):
+# def ecef2aer(x: float, y: float, z: float, lat0: float, lon0: float, h0: float, ell: Ellipsoid = None, deg: bool = True):
+def ecef2aer(x, y, z, lat0, lon0, h0, ell=None, deg=True):
     """
     gives azimuth, elevation and slant range from an Observer to a Point with ECEF coordinates.
 
@@ -48,9 +47,8 @@ def ecef2aer(x: float, y: float, z: float,
     return enu2aer(xEast, yNorth, zUp, deg=deg)
 
 
-def geodetic2aer(lat: float, lon: float, h: float,
-                 lat0: float, lon0: float, h0: float,
-                 ell: Ellipsoid = None, deg: bool = True):
+# def geodetic2aer(lat: float, lon: float, h: float, lat0: float, lon0: float, h0: float, ell: Ellipsoid = None, deg: bool = True):
+def geodetic2aer(lat, lon, h, lat0, lon0, h0, ell=None, deg=True):
     """
     gives azimuth, elevation and slant range from an Observer to a Point with geodetic coordinates.
 
@@ -89,10 +87,8 @@ def geodetic2aer(lat: float, lon: float, h: float,
     return enu2aer(e, n, u, deg=deg)
 
 
-def aer2geodetic(az: float, el: float, srange: float,
-                 lat0: float, lon0: float, h0: float,
-                 ell: Ellipsoid = None,
-                 deg: bool = True):
+# def aer2geodetic(az: float, el: float, srange: float, lat0: float, lon0: float, h0: float, ell: Ellipsoid = None, deg: bool = True):
+def aer2geodetic(az, el, srange, lat0, lon0, h0, ell=None, deg=True):
     """
     gives geodetic coordinates of a point with az, el, range
     from an observer at lat0, lon0, h0
@@ -133,9 +129,8 @@ def aer2geodetic(az: float, el: float, srange: float,
     return ecef2geodetic(x, y, z, ell=ell, deg=deg)
 
 
-def aer2ecef(az: float, el: float, srange: float,
-             lat0: float, lon0: float, alt0: float,
-             ell: Ellipsoid = None, deg: bool = True):
+# def aer2ecef(az: float, el: float, srange: float, lat0: float, lon0: float, alt0: float, ell: Ellipsoid = None, deg: bool = True):
+def aer2ecef(az, el, srange, lat0, lon0, alt0, ell=None, deg=True):
     """
     converts target azimuth, elevation, range from observer at lat0,lon0,alt0 to ECEF coordinates.