code: def compute_pixels(orb, sgeom, times, rpy=(0.0, 0.0, 0.0)):
      
      if isinstance(orb, (list, tuple)):
          tle1, tle2 = orb
          orb = Orbital("mysatellite", line1=tle1, line2=tle2)
  
      # get position and velocity for each time of each pixel
      pos, vel = orb.get_position(times, normalize=False)
  
      # now, get the vectors pointing to each pixel
      vectors = sgeom.vectors(pos, vel, *rpy)
  
      # compute intersection of lines (directed by vectors and passing through
      # (0, 0, 0)) and ellipsoid. Derived from:
      # http://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
  
      # do the computation between line and ellipsoid (WGS 84)
      # NB: AAPP uses GRS 80...
      centre = -pos
      a__ = 6378.137  # km
      # b__ = 6356.75231414 # km, GRS80
      b__ = 6356.752314245  # km, WGS84
      radius = np.array([[1 / a__, 1 / a__, 1 / b__]]).T
      shape = vectors.shape
  
      xr_ = vectors.reshape([3, -1]) * radius
      cr_ = centre.reshape([3, -1]) * radius
      ldotc = np.einsum("ij,ij->j", xr_, cr_)
      lsq = np.einsum("ij,ij->j", xr_, xr_)
      csq = np.einsum("ij,ij->j", cr_, cr_)
  
      d1_ = (ldotc - np.sqrt(ldotc ** 2 - csq * lsq + lsq)) / lsq
  
      # return the actual pixel positions
      return vectors * d1_.reshape(shape[1:]) - centre
docs: Compute cartesian coordinates of the pixels in instrument scan.