def METHOD_NAME(mats, mean_matrix):
  """Given the sequence of matrices mats, calculate the variance-covariance
  matrix of the elements, using the known mean values for each elt in mat
  to avoid the approximation implied in taking the sample covariance"""
  # check input
  assert all([e.is_square() for e in mats])
  n = mats[0].n_rows()
  assert all([e.n_rows() == n for e in mats])
  # create an empty var-cov matrix
  covmat = flex.double(flex.grid(n**2,n**2), 0.0)
  for i in range(covmat.all()[0]):
    for j in range(covmat.all()[1]):
      resid_a = flex.double([m[i] - mean_matrix[i] for m in mats])
      resid_b = flex.double([m[j] - mean_matrix[j] for m in mats])
      covmat[i,j] = flex.mean(resid_a*resid_b)
  return covmat