def plotPotentials(Pot,rmin=0.,rmax=1.5,nrs=21,zmin=-0.5,zmax=0.5,nzs=21,
                   phi=None,xy=False,t=0.,effective=False,Lz=None,
                   ncontours=21,savefilename=None,aspect=None,
                   justcontours=False,levels=None,cntrcolors=None):
        """
        NAME:

           plotPotentials

        PURPOSE:

           plot a set of potentials

        INPUT:

           Pot - Potential or list of Potential instances

           rmin= minimum R (can be Quantity) [xmin if xy]

           rmax= maximum R (can be Quantity) [ymax if xy]

           nrs= grid in R

           zmin= minimum z (can be Quantity) [ymin if xy]

           zmax= maximum z (can be Quantity) [ymax if xy]

           nzs= grid in z

           phi= (None) azimuth to use for non-axisymmetric potentials

           t= (0.) time to use to evaluate potential

           xy= (False) if True, plot the potential in X-Y

           effective= (False) if True, plot the effective potential Phi + Lz^2/2/R^2

           Lz= (None) angular momentum to use for the effective potential when effective=True

           justcontours= (False) if True, just plot contours

           levels= (None) contours to plot

           ncontours - number of contours when levels is None

           cntrcolors= (None) colors of the contours (single color or array with length ncontours)

           savefilename= save to or restore from this savefile (pickle)

        OUTPUT:

           plot to output device

        HISTORY:

           2010-07-09 - Written - Bovy (NYU)

        """
        Pot= flatten(Pot)
        if _APY_LOADED:
            if hasattr(Pot,'_ro'):
                tro= Pot._ro
            else:
                tro= Pot[0]._ro
            if isinstance(rmin,units.Quantity):
                rmin= rmin.to(units.kpc).value/tro
            if isinstance(rmax,units.Quantity):
                rmax= rmax.to(units.kpc).value/tro
            if isinstance(zmin,units.Quantity):
                zmin= zmin.to(units.kpc).value/tro
            if isinstance(zmax,units.Quantity):
                zmax= zmax.to(units.kpc).value/tro
        if not savefilename == None and os.path.exists(savefilename):
            print("Restoring savefile "+savefilename+" ...")
            savefile= open(savefilename,'rb')
            potRz= pickle.load(savefile)
            Rs= pickle.load(savefile)
            zs= pickle.load(savefile)
            savefile.close()
        else:
            if effective and Lz is None:
                raise RuntimeError("When effective=True, you need to specify Lz=")
            Rs= nu.linspace(rmin,rmax,nrs)
            zs= nu.linspace(zmin,zmax,nzs)
            potRz= nu.zeros((nrs,nzs))
            for ii in range(nrs):
                for jj in range(nzs):
                    if xy:
                        R,phi,z= bovy_coords.rect_to_cyl(Rs[ii],zs[jj],0.)
                    else:
                        R,z= Rs[ii], zs[jj]
                    potRz[ii,jj]= evaluatePotentials(Pot,nu.fabs(R),
                                                     z,phi=phi,t=t,
                                                     use_physical=False)
                if effective:
                    potRz[ii,:]+= 0.5*Lz**2/Rs[ii]**2.
            if not savefilename == None:
                print("Writing savefile "+savefilename+" ...")
                savefile= open(savefilename,'wb')
                pickle.dump(potRz,savefile)
                pickle.dump(Rs,savefile)
                pickle.dump(zs,savefile)
                savefile.close()
        if aspect is None:
            aspect=.75*(rmax-rmin)/(zmax-zmin)
        if xy:
            xlabel= r'$x/R_0$'
            ylabel= r'$y/R_0$'
        else:
            xlabel=r"$R/R_0$"
            ylabel=r"$z/R_0$"
        if levels is None:
            levels= nu.linspace(nu.nanmin(potRz),nu.nanmax(potRz),ncontours)
        if cntrcolors is None:
            cntrcolors= 'k'
        return plot.bovy_dens2d(potRz.T,origin='lower',cmap='gist_gray',contours=True,
                                xlabel=xlabel,ylabel=ylabel,
                                aspect=aspect,
                                xrange=[rmin,rmax],
                                yrange=[zmin,zmax],
                                cntrls='-',
                                justcontours=justcontours,
                                levels=levels,cntrcolors=cntrcolors)
NAME:

           plotPotentials

        PURPOSE:

           plot a set of potentials

        INPUT:

           Pot - Potential or list of Potential instances

           rmin= minimum R (can be Quantity) [xmin if xy]

           rmax= maximum R (can be Quantity) [ymax if xy]

           nrs= grid in R

           zmin= minimum z (can be Quantity) [ymin if xy]

           zmax= maximum z (can be Quantity) [ymax if xy]

           nzs= grid in z

           phi= (None) azimuth to use for non-axisymmetric potentials

           t= (0.) time to use to evaluate potential

           xy= (False) if True, plot the potential in X-Y

           effective= (False) if True, plot the effective potential Phi + Lz^2/2/R^2

           Lz= (None) angular momentum to use for the effective potential when effective=True

           justcontours= (False) if True, just plot contours

           levels= (None) contours to plot

           ncontours - number of contours when levels is None

           cntrcolors= (None) colors of the contours (single color or array with length ncontours)

           savefilename= save to or restore from this savefile (pickle)

        OUTPUT:

           plot to output device

        HISTORY:

           2010-07-09 - Written - Bovy (NYU)