docs: def elld2 (x0, y0, mjr, mnr, pa, x, y):
      """Given an 2D Gaussian expressed as an ellipse (major, minor, pa), compute a
      "squared distance parameter" such that
  
         z = exp (-0.5 * d2)
  
      Inputs:
  
      * x0: position of Gaussian center on x axis
      * y0: position of Gaussian center on y axis
      * mjr: major axis (sigma not FWHM) of the Gaussian
      * mnr: minor axis (sigma not FWHM) of the Gaussian
      * pa: position angle (from +x to +y) of the Gaussian, radians
      * x: x coordinates of the locations for which to evaluate d2
      * y: y coordinates of the locations for which to evaluate d2
  
      Returns: d2, distance parameter defined as above.
  
      x0, y0, mjr, and mnr may be in any units so long as they're consistent. x
      and y may be arrays (of the same shape), in which case d2 will be an array
      as well.
  
      """
      _ellcheck (mjr, mnr, pa)
  
      dx, dy = x - x0, y - y0
      c, s = np.cos (pa), np.sin (pa)
      a = c * dx + s * dy
      b = -s * dx + c * dy
      return (a / mjr)**2 + (b / mnr)**2