docs: def potinf(self, x, y, aq=None):
          '''
          linesink 0, order 0, layer[0]
                      order 0, layer[1]
                      ...
                      order 1, layer[0]
                      order 1, layer[1]
                      ...
          linesink 1, order 0, layer[0]
                      order 0, layer[1]
                      ...
                      order 1, layer[0]
                      order 1, layer[1]
                      ...
          '''    
          if aq is None: aq = self.model.aq.find_aquifer_data(x, y)
          rv = np.zeros((self.nls, self.lslist[0].nparam, aq.naq))
          if aq in self.aq:
              for i, ls in enumerate(self.lslist):
                  rv[i] = ls.potinf(x, y, aq)
          rv.shape = (self.nparam, aq.naq)
          return rv