Fix the cubic coefficient in gmst()

pyorbital/astronomy.py

@@ -50,7 +50,7 @@ def gmst(utc_time):
     """
     ut1 = jdays2000(utc_time) / 36525.0
     theta = 67310.54841 + ut1 * (876600 * 3600 + 8640184.812866 + ut1 *
-                                 (0.093104 - ut1 * 6.2 * 10e-6))
+                                 (0.093104 - ut1 * 6.2e-6))
     return np.deg2rad(theta / 240.0) % (2 * np.pi)
 
 

pyorbital/orbital.py

@@ -145,6 +145,7 @@ def get_observer_look_from_cartesian_position(utc_time, lon, lat, alt, pos_x, po
     # Due to rounding top_z can be larger than rg_ (when el_ ~ 90).
     top_z_divided_by_rg_ = top_z_divided_by_rg_.clip(max=1)
     el_ = np.arcsin(top_z_divided_by_rg_)
+    az_ = az_ + (np.pi - az_) * (top_z_divided_by_rg_ == 1)
     return np.rad2deg(az_), np.rad2deg(el_)
 
 

pyorbital/tests/test_astronomy.py

@@ -130,6 +130,15 @@ def test_sun_zenith_angle_integer_array(self, dtype):
         assert np.issubdtype(result.dtype, np.floating)
         np.testing.assert_allclose(result, expected, atol=1e-8)
 
+    def test_gmst_far_from_j2000(self):
+        """Test GMST against the AIAA 2006 gstime polynomial two centuries from J2000."""
+        utc_time = dt.datetime(2200, 1, 1, 12, 0)
+        t = astr.jdays2000(utc_time) / 36525.0
+        theta = (-6.2e-6 * t ** 3 + 0.093104 * t ** 2 + (876600 * 3600 + 8640184.812866) * t
+                 + 67310.54841)
+        expected = np.deg2rad(theta / 240.0) % (2 * np.pi)
+        assert astr.gmst(utc_time) == pytest.approx(expected, abs=1e-10)
+
     def test_sun_earth_distance_correction(self):
         """Test the sun-earth distance correction."""
         utc_time = dt.datetime(2022, 6, 15, 12, 0, 0)