def bovy_dens2d(X,**kwargs):
    """
    NAME:

       bovy_dens2d

    PURPOSE:

       plot a 2d density with optional contours

    INPUT:

       first argument is the density

       matplotlib.pyplot.imshow keywords (see http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.imshow)

       xlabel - (raw string!) x-axis label, LaTeX math mode, no $s needed

       ylabel - (raw string!) y-axis label, LaTeX math mode, no $s needed

       xrange

       yrange

       noaxes - don't plot any axes

       overplot - if True, overplot

       colorbar - if True, add colorbar

       shrink= colorbar argument: shrink the colorbar by the factor (optional)

       conditional - normalize each column separately (for probability densities, i.e., cntrmass=True)

       gcf=True does not start a new figure (does change the ranges and labels)

       Contours:
       
       justcontours - if True, only draw contours

       contours - if True, draw contours (10 by default)

       levels - contour-levels

       cntrmass - if True, the density is a probability and the levels are probability masses contained within the contour

       cntrcolors - colors for contours (single color or array)

       cntrlabel - label the contours

       cntrlw, cntrls - linewidths and linestyles for contour

       cntrlabelsize, cntrlabelcolors,cntrinline - contour arguments

       cntrSmooth - use ndimage.gaussian_filter to smooth before contouring

       onedhists - if True, make one-d histograms on the sides

       onedhistcolor - histogram color

       retAxes= return all Axes instances

       retCont= return the contour instance

    OUTPUT:

       plot to output device, Axes instances depending on input

    HISTORY:

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

    """
    overplot= kwargs.pop('overplot',False)
    gcf= kwargs.pop('gcf',False)
    if not overplot and not gcf:
        pyplot.figure()
    xlabel= kwargs.pop('xlabel',None)
    ylabel= kwargs.pop('ylabel',None)
    zlabel= kwargs.pop('zlabel',None)
    if 'extent' in kwargs:
        extent= kwargs.pop('extent')
    else:
        xlimits= kwargs.pop('xrange',[0,X.shape[1]])
        ylimits= kwargs.pop('yrange',[0,X.shape[0]])
        extent= xlimits+ylimits
    if not 'aspect' in kwargs:
        kwargs['aspect']= (xlimits[1]-xlimits[0])/float(ylimits[1]-ylimits[0])
    noaxes= kwargs.pop('noaxes',False)
    justcontours= kwargs.pop('justcontours',False)
    if ('contours' in kwargs and kwargs['contours']) or \
            'levels' in kwargs or justcontours or \
            ('cntrmass' in kwargs and kwargs['cntrmass']):
        contours= True
    else:
        contours= False
    kwargs.pop('contours',None)
    if 'levels' in kwargs:
        levels= kwargs['levels']
        kwargs.pop('levels')
    elif contours:
        if 'cntrmass' in kwargs and kwargs['cntrmass']:
            levels= sc.linspace(0.,1.,_DEFAULTNCNTR)
        elif True in sc.isnan(sc.array(X)):
            levels= sc.linspace(sc.nanmin(X),sc.nanmax(X),_DEFAULTNCNTR)
        else:
            levels= sc.linspace(sc.amin(X),sc.amax(X),_DEFAULTNCNTR)
    cntrmass= kwargs.pop('cntrmass',False)
    conditional= kwargs.pop('conditional',False)
    cntrcolors= kwargs.pop('cntrcolors','k')
    cntrlabel= kwargs.pop('cntrlabel',False)
    cntrlw= kwargs.pop('cntrlw',None)
    cntrls= kwargs.pop('cntrls',None)
    cntrSmooth= kwargs.pop('cntrSmooth',None)
    cntrlabelsize= kwargs.pop('cntrlabelsize',None)
    cntrlabelcolors= kwargs.pop('cntrlabelcolors',None)
    cntrinline= kwargs.pop('cntrinline',None)
    retCumImage= kwargs.pop('retCumImage',False)
    cb= kwargs.pop('colorbar',False)
    shrink= kwargs.pop('shrink',None)
    onedhists= kwargs.pop('onedhists',False)
    onedhistcolor= kwargs.pop('onedhistcolor','k')
    retAxes= kwargs.pop('retAxes',False)
    retCont= kwargs.pop('retCont',False)
    if onedhists:
        if overplot or gcf: fig= pyplot.gcf()
        else: fig= pyplot.figure()
        nullfmt   = NullFormatter()         # no labels
        # definitions for the axes
        left, width = 0.1, 0.65
        bottom, height = 0.1, 0.65
        bottom_h = left_h = left+width
        rect_scatter = [left, bottom, width, height]
        rect_histx = [left, bottom_h, width, 0.2]
        rect_histy = [left_h, bottom, 0.2, height]
        axScatter = pyplot.axes(rect_scatter)
        axHistx = pyplot.axes(rect_histx)
        axHisty = pyplot.axes(rect_histy)
        # no labels
        axHistx.xaxis.set_major_formatter(nullfmt)
        axHistx.yaxis.set_major_formatter(nullfmt)
        axHisty.xaxis.set_major_formatter(nullfmt)
        axHisty.yaxis.set_major_formatter(nullfmt)
        fig.sca(axScatter)
    ax=pyplot.gca()
    ax.set_autoscale_on(False)
    if conditional:
        plotthis= X/sc.tile(sc.sum(X,axis=0),(X.shape[1],1))
    else:
        plotthis= X
    if not justcontours:
        out= pyplot.imshow(plotthis,extent=extent,**kwargs)
    if not overplot:
        pyplot.axis(extent)
        _add_axislabels(xlabel,ylabel)
        _add_ticks()
    #Add colorbar
    if cb and not justcontours:
        if shrink is None:
            shrink= sc.amin([float(kwargs.pop('aspect',1.))*0.87,1.])
        CB1= pyplot.colorbar(out,shrink=shrink)
        if not zlabel is None:
            if zlabel[0] != '$':
                thiszlabel=r'$'+zlabel+'$'
            else:
                thiszlabel=zlabel
            CB1.set_label(thiszlabel)
    if contours or retCumImage:
        aspect= kwargs.get('aspect',None)
        origin= kwargs.get('origin',None)
        if cntrmass:
            #Sum from the top down!
            plotthis[sc.isnan(plotthis)]= 0.
            sortindx= sc.argsort(plotthis.flatten())[::-1]
            cumul= sc.cumsum(sc.sort(plotthis.flatten())[::-1])/sc.sum(plotthis.flatten())
            cntrThis= sc.zeros(sc.prod(plotthis.shape))
            cntrThis[sortindx]= cumul
            cntrThis= sc.reshape(cntrThis,plotthis.shape)
        else:
            cntrThis= plotthis
        if contours:
            if not cntrSmooth is None:
                cntrThis= ndimage.gaussian_filter(cntrThis,cntrSmooth,
                                                  mode='nearest')
            cont= pyplot.contour(cntrThis,levels,colors=cntrcolors,
                                 linewidths=cntrlw,extent=extent,aspect=aspect,
                                 linestyles=cntrls,origin=origin)
            if cntrlabel:
                pyplot.clabel(cont,fontsize=cntrlabelsize,
                              colors=cntrlabelcolors,
                              inline=cntrinline)
    if noaxes:
        ax.set_axis_off()
    #Add onedhists
    if not onedhists:
        if retCumImage:
            return cntrThis
        elif retAxes:
            return pyplot.gca()
        elif retCont:
            return cont
        elif justcontours:
            return cntrThis
        else:
            return out
    histx= sc.nansum(X.T,axis=1)*m.fabs(ylimits[1]-ylimits[0])/X.shape[1] #nansum bc nan is *no dens value*
    histy= sc.nansum(X.T,axis=0)*m.fabs(xlimits[1]-xlimits[0])/X.shape[0]
    histx[sc.isnan(histx)]= 0.
    histy[sc.isnan(histy)]= 0.
    dx= (extent[1]-extent[0])/float(len(histx))
    axHistx.plot(sc.linspace(extent[0]+dx,extent[1]-dx,len(histx)),histx,
                 drawstyle='steps-mid',color=onedhistcolor)
    dy= (extent[3]-extent[2])/float(len(histy))
    axHisty.plot(histy,sc.linspace(extent[2]+dy,extent[3]-dy,len(histy)),
                 drawstyle='steps-mid',color=onedhistcolor)
    axHistx.set_xlim( axScatter.get_xlim() )
    axHisty.set_ylim( axScatter.get_ylim() )
    axHistx.set_ylim( 0, 1.2*sc.amax(histx))
    axHisty.set_xlim( 0, 1.2*sc.amax(histy))
    if retCumImage:
        return cntrThis
    elif retAxes:
        return (axScatter,axHistx,axHisty)
    elif justcontours:
        return cntrThis
    else:
        return out
NAME:

       bovy_dens2d

    PURPOSE:

       plot a 2d density with optional contours

    INPUT:

       first argument is the density

       matplotlib.pyplot.imshow keywords (see http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.imshow)

       xlabel - (raw string!) x-axis label, LaTeX math mode, no $s needed

       ylabel - (raw string!) y-axis label, LaTeX math mode, no $s needed

       xrange

       yrange

       noaxes - don't plot any axes

       overplot - if True, overplot

       colorbar - if True, add colorbar

       shrink= colorbar argument: shrink the colorbar by the factor (optional)

       conditional - normalize each column separately (for probability densities, i.e., cntrmass=True)

       gcf=True does not start a new figure (does change the ranges and labels)

       Contours:
       
       justcontours - if True, only draw contours

       contours - if True, draw contours (10 by default)

       levels - contour-levels

       cntrmass - if True, the density is a probability and the levels are probability masses contained within the contour

       cntrcolors - colors for contours (single color or array)

       cntrlabel - label the contours

       cntrlw, cntrls - linewidths and linestyles for contour

       cntrlabelsize, cntrlabelcolors,cntrinline - contour arguments

       cntrSmooth - use ndimage.gaussian_filter to smooth before contouring

       onedhists - if True, make one-d histograms on the sides

       onedhistcolor - histogram color

       retAxes= return all Axes instances

       retCont= return the contour instance

    OUTPUT:

       plot to output device, Axes instances depending on input

    HISTORY:

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