docs: def potentiallayers(self, x, y, layers, aq=None):
          '''Returns array of size len(layers)
          only used in building equations'''
          if aq is None: aq = self.model.aq.find_aquifer_data(x, y)
          pot = np.sum(self.potential(x, y, aq) * aq.eigvec, 1 )
          return pot[layers]