def _sampleV_preoptimized(self,R,z,maxVT):
        """
        NAME:

           _sampleV_preoptimized

        PURPOSE:

           sample a radial, azimuthal, and vertical velocity at R,z;
           R,z can be an array of positions maxVT is already optimized

        INPUT:

           R - Galactocentric distance (can be Quantity)

           z - height (can be Quantity)
           
           maxVT - an array of pre-optimized maximum vT at corresponding R,z

        OUTPUT:

           a numpy array containing the sampled velocity, (vR, vT, vz),
           where each row correspond to the row of (R,z)

        HISTORY:

           2018-08-09 - Written - Samuel Wong (University of Toronto)

        """
        length = numpy.size(R)
        out= numpy.empty((length,3)) #Initialize output
        #Determine the maximum of the velocity distribution
        maxVR= numpy.zeros(length)
        maxVz= numpy.zeros(length)
        logmaxVD= self(R,maxVR,maxVT,z,maxVz,log=True,use_physical=False)
        #Now rejection-sample
        #Intiialize boolean index of position remaining to be sampled
        remain_indx = numpy.full(length,True)
        while numpy.any(remain_indx):
            nmore= numpy.sum(remain_indx)
            propvR= numpy.random.normal(size=nmore)*2.*self._sr
            propvT= numpy.random.normal(size=nmore)*2.*self._sr+maxVT[remain_indx]
            propvz= numpy.random.normal(size=nmore)*2.*self._sz
            VDatprop= self(R[remain_indx],propvR,propvT,z[remain_indx],propvz,
                           log=True, use_physical=False)-logmaxVD[remain_indx]
            VDatprop-= -0.5*(propvR**2./4./self._sr**2.+
                             propvz**2./4./self._sz**2.+
                             (propvT-maxVT[remain_indx])**2./4./self._sr**2.)
            accept_indx= (VDatprop > numpy.log(numpy.random.random(size=nmore)))
            vR_accept= propvR[accept_indx]
            vT_accept= propvT[accept_indx]
            vz_accept= propvz[accept_indx]
            #Get the indexing of rows of output array that need to be updated
            #with newly accepted velocity
            to_change= numpy.copy(remain_indx)
            to_change[remain_indx]= accept_indx
            out[to_change]= numpy.stack((vR_accept,vT_accept,vz_accept), axis = 1)
            #Removing accepted sampled from remain index
            remain_indx[remain_indx]= ~accept_indx
        return out
NAME:

           _sampleV_preoptimized

        PURPOSE:

           sample a radial, azimuthal, and vertical velocity at R,z;
           R,z can be an array of positions maxVT is already optimized

        INPUT:

           R - Galactocentric distance (can be Quantity)

           z - height (can be Quantity)
           
           maxVT - an array of pre-optimized maximum vT at corresponding R,z

        OUTPUT:

           a numpy array containing the sampled velocity, (vR, vT, vz),
           where each row correspond to the row of (R,z)

        HISTORY:

           2018-08-09 - Written - Samuel Wong (University of Toronto)