code: def first_and_second_harmonic_function(phi, c):
      
  
      return (c[0] + c[1]*np.sin(phi) + c[2]*np.cos(phi) + c[3]*np.sin(2*phi) +
              c[4]*np.cos(2*phi))
docs: Compute the harmonic function value used to calculate the
      corrections for ellipse fitting.
  
      This function includes simultaneously both the first and second
      order harmonics:
  
      .. math::
  
          f(phi) = c[0] + c[1]*\\sin(phi) + c[2]*\\cos(phi) +
                   c[3]*\\sin(2*phi) + c[4]*\\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.
      c : `~numpy.ndarray` of shape (5,)
          Array containing the five harmonic coefficients.
  
      Returns
      -------
      result : float or `~numpy.ndarray`
          The function value(s) at the given input angle(s).