docs: def _compute_errors(self):
          """
          Compute parameter errors based on the diagonal of the covariance
          matrix of the four harmonic coefficients for harmonics n=1 and
          n=2.
          """
  
          try:
              coeffs = fit_first_and_second_harmonics(self.sample.values[0],
                                                      self.sample.values[2])
              covariance = coeffs[1]
              coeffs = coeffs[0]
              model = first_and_second_harmonic_function(self.sample.values[0],
                                                         coeffs)
              residual_rms = np.std(self.sample.values[2] - model)
              errors = np.diagonal(covariance) * residual_rms
  
              eps = self.sample.geometry.eps
              pa = self.sample.geometry.pa
  
              # parameter errors result from direct projection of
              # coefficient errors.  These showed to be the error estimators
              # that best convey the errors measured in Monte Carlo
              # experiments (see Busko 1996; ASPC 101, 139).
              ea = abs(errors[2] / self.grad)
              eb = abs(errors[1] * (1. - eps) / self.grad)
              self.x0_err = np.sqrt((ea * np.cos(pa))**2 + (eb * np.sin(pa))**2)
              self.y0_err = np.sqrt((ea * np.sin(pa))**2 + (eb * np.cos(pa))**2)
              self.ellip_err = (abs(2. * errors[4] * (1. - eps) / self.sma /
                                    self.grad))
              if (abs(eps) > np.finfo(float).resolution):
                  self.pa_err = (abs(2. * errors[3] * (1. - eps) / self.sma /
                                     self.grad / (1. - (1. - eps)**2)))
              else:
                  self.pa_err = 0.
          except Exception:    # we want to catch everything
              self.x0_err = self.y0_err = self.pa_err = self.ellip_err = 0.