code: def fit_upper_harmonic(phi, intensities, order):
      
  
      an = bn = 1.
  
      def optimize_func(x):
          return (x[0] + x[1]*np.sin(order*phi) + x[2]*np.cos(order*phi) -
                  intensities)
  
      return _least_squares_fit(optimize_func, [np.mean(intensities), an, bn])
docs: Fit upper harmonic function to a set of (angle, intensity) pairs.
  
      With ``order`` set to 3 or 4, the resulting amplitudes, divided by
      the semimajor axis length and local gradient, measure the deviations
      from perfect ellipticity.
  
      The harmonic function that is fit is:
  
      .. math::
          y(phi, order) = y0 + An*\\sin(order*phi) + Bn*\\cos(order*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.
      order : int
          The order of the harmonic to be fitted.
  
      Returns
      -------
      y0, An, Bn : float
          The fitted harmonic values.