code: def reverse_terra(xyz_au, gast, iterations=3):
      
      x, y, z = xyz_au
      R = sqrt(x*x + y*y)
  
      lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
      lat = arctan2(z, R)
  
      a = ERAD / AU_M
      f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
      e2 = 2.0*f - f*f
      i = 0
      while i < iterations:
          i += 1
          C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
          lat = arctan2(z + a * C * e2 * sin(lat), R)
      elevation_m = ((R / cos(lat)) - a * C) * AU_M
      return lat, lon, elevation_m
docs: Convert a geocentric (x,y,z) at time `t` to latitude and longitude.
  
      Returns a tuple of latitude, longitude, and elevation whose units
      are radians and meters.  Based on Dr. T.S. Kelso's quite helpful
      article "Orbital Coordinate Systems, Part III":
      https://www.celestrak.com/columns/v02n03/