docs: def fit_first_and_second_harmonics(phi, intensities):
      """
      Fit the first and second harmonic function values to a set of
      (angle, intensity) pairs.
  
      This function is used to compute corrections for ellipse fitting:
  
      .. math::
  
          f(phi) = y0 + a1*\\sin(phi) + b1*\\cos(phi) + a2*\\sin(2*phi) +
                   b2*\\cos(2*phi)
  
      Parameters
      ----------
      phi : float or `~numpy.ndarray`
          The angle(s) along the elliptical path, going towards the positive
          y axis, starting coincident with the position angle. That is, the
          angles are defined from the semimajor axis that lies in
          the positive x quadrant.
      intensities : `~numpy.ndarray`
          The intensities measured along the elliptical path, at the
          angles defined by the ``phi`` parameter.
  
      Returns
      -------
      y0, a1, b1, a2, b2 : float
          The fitted harmonic coefficent values.
      """
  
      a1 = b1 = a2 = b2 = 1.
  
      def optimize_func(x):
          return first_and_second_harmonic_function(
              phi, np.array([x[0], x[1], x[2], x[3], x[4]])) - intensities
  
      return _least_squares_fit(optimize_func, [np.mean(intensities), a1, b1,
                                                a2, b2])