docs: def potinflayers(self, x, y, layers, aq=None):
          '''Returns array of size (len(layers),nparam)
          only used in building equations'''
          if aq is None: aq = self.model.aq.find_aquifer_data(x, y)
          pot = self.potinf(x, y, aq)  # nparam rows, naq cols
          rv = np.sum(pot[:,np.newaxis,:] * aq.eigvec, 2).T  # Transopose as the first axes needs to be the number of layers
          return rv[layers,:]