add fix for elevation at the poles

Author: jamesgrimmett

skyfield/toposlib.py

@@ -144,6 +144,7 @@ def __init__(self, name, radius_m, inverse_flattening):
         self._one_minus_flattening_squared = omf * omf
         f = 1.0 / inverse_flattening
         self._e2 = 2.0*f - f*f
+        self.semi_minor_axis = Distance(m=radius_m * (1.0 - f))
 
     def latlon(self, latitude_degrees, longitude_degrees, elevation_m=0.0,
                cls=GeographicPosition):
@@ -211,7 +212,10 @@ def height_of(self, position):
 
         """
         xyz_au, x, y, aC, R, lat = self._compute_latitude(position)
-        height_au = R / cos(lat) - aC
+        if R > 0.0:
+            height_au = R / cos(lat) - aC
+        else:
+            height_au = abs(xyz_au[2]) - self.semi_minor_axis.au
         return Distance(height_au)
 
     def geographic_position_of(self, position):
@@ -225,7 +229,10 @@ def geographic_position_of(self, position):
         """
         xyz_au, x, y, aC, R, lat = self._compute_latitude(position)
         lon = (arctan2(y, x) - pi) % tau - pi
-        height_au = R / cos(lat) - aC
+        if R > 0.0:
+            height_au = R / cos(lat) - aC
+        else:
+            height_au = abs(xyz_au[2]) - self.semi_minor_axis.au
         return GeographicPosition(
             latitude=Angle(radians=lat),
             longitude=Angle(radians=lon),