code: def sampleVRVT(self,R,n=1,nsigma=None,target=True):
          
          #Determine where the max of the v-distribution is using asymmetric drift
          maxVR= 0.
          maxVT= optimize.brentq(_vtmaxEq,0.,R**self._beta+0.2,(R,self))
          maxVD= self(Orbit([R,maxVR,maxVT]))
          #Now rejection-sample
          if nsigma == None:
              nsigma= _NSIGMA
          out= []
          if target:
              sigma= math.sqrt(self.targetSigma2(R,use_physical=False))
          else:
              sigma= math.sqrt(self.sigma2(R,use_physical=False))
          while len(out) < n:
              #sample
              vrg, vtg= nu.random.normal(), nu.random.normal()
              propvR= vrg*nsigma*sigma
              propvT= vtg*nsigma*sigma/self._gamma+maxVT
              VDatprop= self(Orbit([R,propvR,propvT]))
              if VDatprop/maxVD > nu.random.uniform()*nu.exp(-0.5*(vrg**2.+vtg**2.)): #accept
                  out.append(sc.array([propvR,propvT]))
          return nu.array(out)
docs: NAME:
  
             sampleVRVT
  
          PURPOSE:
  
             sample a radial and azimuthal velocity at R
  
          INPUT:
  
             R - Galactocentric distance (can be Quantity)
  
             n= number of distances to sample
  
             nsigma= number of sigma to rejection-sample on
  
             target= if True, sample using the 'target' sigma_R rather than the actual sigma_R (default=True)
  
          OUTPUT:
  
             list of samples
  
          BUGS:
  
             should use the fact that vR and vT separate
  
          HISTORY:
  
             2011-03-24 - Written - Bovy (NYU)