code: def _matrix_sigma_eps(self, donor_catchments):
          
          p = len(donor_catchments)
          sigma = np.empty((p, p))
          for i in range(p):
              beta_i = self._beta(donor_catchments[i])
              n_i = donor_catchments[i].amax_records_end() - donor_catchments[i].amax_records_start() + 1
              for j in range(p):
                  beta_j = self._beta(donor_catchments[j])
                  n_j = donor_catchments[j].amax_records_end() - donor_catchments[j].amax_records_start() + 1
                  rho_ij = self._lnqmed_corr(donor_catchments[i], donor_catchments[j])
                  n_ij = min(donor_catchments[i].amax_records_end(), donor_catchments[j].amax_records_end()) - \
                         max(donor_catchments[i].amax_records_start(), donor_catchments[j].amax_records_start()) + 1
                  sigma[i, j] = 4 * beta_i * beta_j * n_ij / n_i / n_j * rho_ij
          return sigma
docs: Return sampling error coveriance matrix Sigma eta
  
          Methodology source: Kjeldsen & Jones 2009, eq 9
  
          :param donor_catchments: Catchments to use as donors
          :type donor_catchments: list of :class:`Catchment`
          :return: 2-Dimensional, symmetric covariance matrix
          :rtype: :class:`numpy.ndarray`