handle vector input for height computation, move to new function

Author: jamesgrimmett

skyfield/tests/test_topos.py

@@ -260,7 +260,7 @@ def test_intersection_from_pole(ts):
     error_degrees = abs(earth_point.latitude.degrees - 90.0)
     error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
     assert error_mas < 0.1
-    assert earth_point.elevation.m < 0.01
+    assert earth_point.elevation.m < 0.1
 
 def test_intersection_from_equator(ts):
     t = ts.utc(2018, 1, 19, 14, 37, 55)
@@ -275,7 +275,7 @@ def test_intersection_from_equator(ts):
     error_degrees = abs(earth_point.longitude.degrees - 0.0)
     error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
     assert error_mas < 0.1
-    assert earth_point.elevation.m < 0.01
+    assert earth_point.elevation.m < 0.1
 
 def test_limb_intersection_points(ts):
     t = ts.utc(2018, 1, 19, 14, 37, 55)
@@ -307,23 +307,23 @@ def test_limb_intersection_points(ts):
 
     error_degrees = abs(intersection_bottom.latitude.degrees + 90.0)
     assert error_degrees < 0.1
-    assert intersection_bottom.elevation.m < 0.01
+    assert intersection_bottom.elevation.m < 0.1
 
     error_degrees = abs(intersection_top.latitude.degrees - 90.0)
     assert error_degrees < 0.1
-    assert intersection_top.elevation.m < 0.01
+    assert intersection_top.elevation.m < 0.1
 
     error_degrees = abs(intersection_left.latitude.degrees - 0.0)
     assert error_degrees < 0.1
 
     error_degrees = abs(intersection_left.longitude.degrees - left_limb_lon)
     assert error_degrees < 0.1
-    assert intersection_left.elevation.m < 0.01
+    assert intersection_left.elevation.m < 0.1
 
     error_degrees = abs(intersection_right.latitude.degrees - 0.0)
     assert error_degrees < 0.1
 
     error_degrees = abs(intersection_right.longitude.degrees - right_limb_lon)
     assert error_degrees < 0.1
-    assert intersection_right.elevation.m < 0.01
+    assert intersection_right.elevation.m < 0.1
 

skyfield/toposlib.py

@@ -1,6 +1,9 @@
 # -*- coding: utf-8 -*-
 
-from numpy import arctan2, array, array2string, cos, exp, sin, sqrt
+from numpy import (
+    arctan2, array, array2string, asarray, cos, exp, float64, sin, sqrt
+)
+
 from .constants import ANGVEL, DAY_S, RAD2DEG, T0, pi, tau
 from .earthlib import refract
 from .framelib import itrs
@@ -212,12 +215,18 @@ def height_of(self, position):
 
         """
         xyz_au, x, y, aC, R, lat = self._compute_latitude(position)
-        if R > 1.e-12:
-            height_au = R / cos(lat) - aC
-        else:
-            height_au = abs(xyz_au[2]) - self.semi_minor_axis.au
+        height_au = self._compute_height(R, lat, aC, xyz_au)
         return Distance(height_au)
 
+    def _compute_height(self, R, lat, aC, xyz_au):
+        height_au = asarray(R)
+        near_pole = height_au <= 1.e-12
+        if near_pole.any():
+            height_au[near_pole] = abs(xyz_au[2]) - self.semi_minor_axis.au
+        else:
+            height_au[~near_pole] = R / cos(lat) - aC
+        return float64(height_au)
+
     def geographic_position_of(self, position):
         """Return the `GeographicPosition` of a ``position``.
 
@@ -229,10 +238,7 @@ def geographic_position_of(self, position):
         """
         xyz_au, x, y, aC, R, lat = self._compute_latitude(position)
         lon = (arctan2(y, x) - pi) % tau - pi
-        if R > 1.e-12:
-            height_au = R / cos(lat) - aC
-        else:
-            height_au = abs(xyz_au[2]) - self.semi_minor_axis.au
+        height_au = self._compute_height(R, lat, aC, xyz_au)
         return GeographicPosition(
             latitude=Angle(radians=lat),
             longitude=Angle(radians=lon),
@@ -280,12 +286,13 @@ def intersection_of(self, position, direction):
         """Return the surface point intersected by a vector.
 
         The ``position`` and ``direction`` describe the tail and orientation
-        of the vector to intersect the ellipsoid, respectively. The vector must
-        be provided in ICRF coordinates with a ``.center``` of 399, the center
-        of the Earth. The direction should be an |xyz| unit vector. Returns a
-        `GeographicPosition` giving the geodetic ``latitude`` and ``longitude``
-        at the point that the ray intersects the surface of the Earth.
-    
+        of the vector to intersect the ellipsoid, respectively. The position
+        must be provided in ICRF coordinates with a ``.center``` of 399, the
+        center of the Earth. The direction should be an |xyz| unit vector in
+        the same coordinate frame. Returns a `GeographicPosition` giving the
+        geodetic ``latitude`` and ``longitude`` at the point that the ray
+        intersects the surface of the Earth.
+
         The main calculation implemented here is based on JPL's NAIF toolkit;
         https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/ellipses.html
         """
@@ -309,7 +316,7 @@ def intersection_of(self, position, direction):
 
         # Define a distortion matrix to map the ellipsoid into a unit sphere
         a = self.radius.au
-        c = a * (1.0 - 1.0 / self.inverse_flattening)
+        c = self.semi_minor_axis.au
         d_matrix = array(
             [[1.0 / a, 0.0, 0.0], [0.0, 1.0 / a, 0.0], [0.0, 0.0, 1.0 / c]]
             )
@@ -344,9 +351,9 @@ def intersection_of(self, position, direction):
             )
 
         # Retrieve latitude and longitude
-        intersection_latlon = self.geographic_position_of(intersection)
+        geo_intersection = self.geographic_position_of(intersection)
 
-        return intersection_latlon
+        return geo_intersection
 
 wgs84 = Geoid('WGS84', 6378137.0, 298.257223563)
 iers2010 = Geoid('IERS2010', 6378136.6, 298.25642)