code: def terra(latitude, longitude, elevation, gast):
      
      zero = zeros_like(gast)
      sinphi = sin(latitude)
      cosphi = cos(latitude)
      c = 1.0 / sqrt(cosphi * cosphi +
                     sinphi * sinphi * one_minus_flattening_squared)
      s = one_minus_flattening_squared * c
      ach = earth_radius_au * c + elevation
      ash = earth_radius_au * s + elevation
  
      # Compute local sidereal time factors at the observer's longitude.
  
      stlocl = 15.0 * DEG2RAD * gast + longitude
      sinst = sin(stlocl)
      cosst = cos(stlocl)
  
      # Compute position vector components in kilometers.
  
      ac = ach * cosphi
      acsst = ac * sinst
      accst = ac * cosst
      pos = array((accst, acsst, zero + ash * sinphi))
  
      # Compute velocity vector components in kilometers/sec.
  
      vel = ANGVEL * DAY_S * array((-acsst, accst, zero))
  
      return pos, vel
docs: Compute the position and velocity of a terrestrial observer.
  
      `latitude` - Latitude in radians.
      `longitude` - Longitude in radians.
      `elevation` - Elevation above sea level in au.
      `gast` - Hours of Greenwich Apparent Sidereal Time (can be an array).
  
      The return value is a tuple of two 3-vectors `(pos, vel)` in the
      dynamical reference system (the true equator and equinox of date)
      whose components are measured in au with respect to the center of
      the Earth.