def disvecinf(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((2, self.nparam, aq.naq))
        if aq == self.aq:
            qxqyrv = rv.reshape((2, self.order + 1, self.nlayers,
                                 aq.naq))  # clever way of using a reshaped rv here
            qxqy = np.zeros((2 * (self.order + 1), aq.naq))
            qxqy[:, :] = self.bessel.disbesldv(float(x), float(y), self.z1, self.z2, aq.lab,
                                               self.order, aq.ilap, aq.naq)
            qxqyrv[0, :] = self.aq.coef[self.layers] * qxqy[:self.order + 1,
                                                       np.newaxis, :]
            qxqyrv[1, :] = self.aq.coef[self.layers] * qxqy[self.order + 1:,
                                                       np.newaxis, :]
        return rv