docs: def potinf(self, x, y, aq=None):
          '''Can be called with only one x,y value
          Returns array(nparam, self.aq.naq) with order
          order 0, layer[0]
          order 0, layer[1]
          ...
          order 1, layer[0]
          order 1, layer[1]
          etc
          '''
          if aq is None: aq = self.model.aq.find_aquifer_data(x, y)
          rv = np.zeros((self.nparam, aq.naq))
          if aq == self.aq:
              # clever way of using a reshaped rv here
              potrv = rv.reshape((self.order + 1, self.nlayers, aq.naq))
              pot = np.zeros((self.order + 1, aq.naq))
              pot[:, :] = self.bessel.potbeslsv(float(x), float(y), self.z1, self.z2, aq.lab,
                                                self.order, aq.ilap, aq.naq)
              potrv[:] = self.aq.coef[self.layers] * pot[:, np.newaxis, :]
          return rv