def pvz(self,vz,R,z,gl=True,ngl=_DEFAULTNGL2,
            nsigma=4.,vTmax=1.5,
            _return_actions=False,_jr=None,_lz=None,_jz=None,
            _return_freqs=False,
            _rg=None,_kappa=None,_nu=None,_Omega=None,
            _sigmaR1=None):
        """
        NAME:

           pvz

        PURPOSE:

           calculate the marginalized vz probability at this location (NOT normalized by the density)

        INPUT:

           vz - vertical velocity (can be Quantity)

           R - radius (can be Quantity)

           z - height (can be Quantity)

           gl - use Gauss-Legendre integration (True, currently the only option)

           ngl - order of Gauss-Legendre integration

           nsigma - sets integration limits to [-1,+1]*nsigma*sigma_R(R) for integration over vR (default: 4)

           vTmax - sets integration limits to [0,vTmax] for integration over vT (default: 1.5)

        OUTPUT:

           p(vz,R,z)

        HISTORY:

           2012-12-22 - Written - Bovy (IAS)

        """
        if _sigmaR1 is None:
            sigmaR1= self._sr*numpy.exp((self._refr-R)/self._hsr)
        else:
            sigmaR1= _sigmaR1
        if gl:
            if ngl % 2 == 1:
                raise ValueError("ngl must be even")
            #Use Gauss-Legendre integration for all
            if ngl == _DEFAULTNGL:
                glx, glw= self._glxdef, self._glwdef
                glx12, glw12= self._glxdef12, self._glwdef12
            elif ngl == _DEFAULTNGL2:
                glx, glw= self._glxdef2, self._glwdef2
                glx12, glw12= self._glxdef, self._glwdef
            else:
                glx, glw= numpy.polynomial.legendre.leggauss(ngl)
                glx12, glw12= numpy.polynomial.legendre.leggauss(ngl//2)
            #Evaluate everywhere
            if isinstance(self._aA,(actionAngle.actionAngleAdiabatic,
                                    actionAngle.actionAngleAdiabaticGrid)):
                vRgl= (glx+1.)
                vRglw= glw
                vRfac= nsigma*sigmaR1  #2 x integration over [0,nsigma*sigmaR1]
            else:
                vRgl= (glx12+1.)
                vRgl= list(vRgl)
                vRgl.extend(-(glx12+1.))
                vRgl= numpy.array(vRgl)
                vRglw= glw12
                vRglw= list(vRglw)
                vRglw.extend(glw12)
                vRglw= numpy.array(vRglw)
                vRfac = 0.5*nsigma*sigmaR1   #integration over [-nsigma*sigmaR1,0] and [0,nsigma*sigmaR1]
            vTgl= vTmax/2.*(glx+1.)
            vTfac= 0.5 * vTmax  #integration over [0.,vTmax]
            #Tile everything
            vTgl= numpy.tile(vTgl,(ngl,1)).T
            vRgl= numpy.tile(vRgl,(ngl,1))
            vTglw= numpy.tile(glw,(ngl,1)).T #also tile weights
            vRglw= numpy.tile(vRglw,(ngl,1))
            #If inputs are arrays, tile
            if isinstance(R,numpy.ndarray):
                nR= len(R)
                R= numpy.tile(R,(ngl,ngl,1)).T.flatten()
                z= numpy.tile(z,(ngl,ngl,1)).T.flatten()
                vz= numpy.tile(vz,(ngl,ngl,1)).T.flatten()
                vTgl= numpy.tile(vTgl,(nR,1,1)).flatten()
                vRgl= numpy.tile(vRgl,(nR,1,1)).flatten()
                vTglw= numpy.tile(vTglw,(nR,1,1))
                vRglw= numpy.tile(vRglw,(nR,1,1))
                scalarOut= False
            else:
                R= R+numpy.zeros(ngl*ngl)
                z= z+numpy.zeros(ngl*ngl)
                vz= vz+numpy.zeros(ngl*ngl)
                nR= 1
                scalarOut= True
                vRgl= vRgl.flatten()
            vRgl*= numpy.tile(nsigma*sigmaR1/2.,(ngl,ngl,1)).T.flatten()
            #evaluate
            if _jr is None and _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R,
                                                                 vRgl.flatten(),
                                                                 vTgl.flatten(),
                                                                 z,
                                                                 vz,
                                                                 log=True,
                                                                 _return_actions=True,
                                                                 _return_freqs=True,
                                                                 use_physical=False)
                logqeval= numpy.reshape(logqeval,(nR,ngl*ngl))
            elif not _jr is None and not _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                           rg=_rg,kappa=_kappa,nu=_nu,
                                           Omega=_Omega,
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(nR,ngl*ngl))
            elif not _jr is None and _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self((_jr,_lz,_jz),
                                           log=True,
                                           _return_actions=True,
                                           _return_freqs=True,
                                           use_physical=False)
                logqeval= numpy.reshape(logqeval,(nR,ngl*ngl))
            elif _jr is None and not _rg is None:
                logqeval, jr, lz, jz, rg, kappa, nu, Omega= self(R,
                                             vRgl.flatten(),
                                             vTgl.flatten(),
                                             z,
                                             vz,
                                             rg=_rg,kappa=_kappa,nu=_nu,
                                             Omega=_Omega,
                                             log=True,
                                             _return_actions=True,
                                             _return_freqs=True,
                                             use_physical=False)
                logqeval= numpy.reshape(logqeval,(nR,ngl*ngl))
            vRglw= numpy.reshape(vRglw,(nR,ngl*ngl))
            vTglw= numpy.reshape(vTglw,(nR,ngl*ngl))
            if scalarOut:
                result= numpy.sum(numpy.exp(logqeval)*vTglw*vRglw,axis=1)[0]*vRfac*vTfac
            else:
                result= numpy.sum(numpy.exp(logqeval)*vTglw*vRglw,axis=1)*vRfac*vTfac
            if _return_actions and _return_freqs:
                return (result,
                        jr,lz,jz,
                        rg, kappa, nu, Omega)          
            elif _return_freqs:
                return (result,
                        rg, kappa, nu, Omega)
            elif _return_actions:
                return (result,
                        jr,lz,jz)
            else:
                return result
NAME:

           pvz

        PURPOSE:

           calculate the marginalized vz probability at this location (NOT normalized by the density)

        INPUT:

           vz - vertical velocity (can be Quantity)

           R - radius (can be Quantity)

           z - height (can be Quantity)

           gl - use Gauss-Legendre integration (True, currently the only option)

           ngl - order of Gauss-Legendre integration

           nsigma - sets integration limits to [-1,+1]*nsigma*sigma_R(R) for integration over vR (default: 4)

           vTmax - sets integration limits to [0,vTmax] for integration over vT (default: 1.5)

        OUTPUT:

           p(vz,R,z)

        HISTORY:

           2012-12-22 - Written - Bovy (IAS)