docs: def potinf(self, x, y, aq=None):
          '''Can be called with only one x,y value'''
          if aq is None: aq = self.model.aq.find_aquifer_data(x, y)
          rv = np.zeros((self.nparam, aq.naq))
          if aq == self.aq:
              pot = np.zeros(aq.naq)
              pot[:] = self.bessel.potbeslsho(float(x), float(y), self.z1, self.z2, aq.lab, 0,
                                              aq.ilap, aq.naq)
              rv[:] = self.aq.coef[self.layers] * pot
          return rv