FUNCTION densityplot,datax,datay,rangex=rangex,rangey=rangey,nbins=nbins
;
; Creates a density array (2D histogram) of a scatter plot.
; INPUT:    datax - array to be plotted along x-axis
;           datay - array to be plotted along y-axis
; OUTPUT:   2D density image of the scatter plot (h),
;           normalised to maximum
; KEYWORDS: rangex = xrange, [x1,x2], default minmax(datax)
;           rangey = yrange, [y1,y2], default minmax(datax)
;           nbins - number of bins per axis (default 50)
;
; METHOD: Derived from STATCONT.PRO from RJLIB/CUBELIB by Rob Rutten;
;         1D histogram converted to 2D
; NOTE:   Before running this function, a plot frame must be made,
;         to define the !x,y.crange. For example,
;         plot,[0,0],/nodata,xra=[0.5,1.1],yra=[0,3500],...

; keyword defaults
if (n_elements(rangex) eq 0) then rangex=minmax(datax)
if (n_elements(rangey) eq 0) then rangey=minmax(datay)
if (n_elements(nbins) eq 0) then nbins=50

dx=!x.crange[1]-!x.crange[0]
dy=!y.crange[1]-!y.crange[0]

; number of bins
nx=nbins
ny=nbins

; binsize in x and y, avoid discretization issues
bx=float(dx)/(nx-1.)
if (size(datax,/type) ne 4 and size(datax,/type) ne 5) then begin
  bx=fix(bx)
  nx=fix(dx/bx)+1
endif
by=float(dy)/(ny-1.)
if (size(datay,/type) ne 4 and size(datay,/type) ne 5) then begin
  by=fix(by)
  ny=fix(dy/by)+1
endif

; scale data
if (!x.crange[0] ne 0) then xd=datax-!x.crange[0] else xd=datax
if (bx ne 1) then xd/=bx
if (!y.crange[0] ne 0) then yd=datay-!y.crange[0] else yd=datay
if (by ne 1) then yd/=by

; create histogram array (1D)
h=long(nx)*fix(yd)+fix(xd)

; set points outside of !x.crange and !y.crange to -1
range=xd ge 0 and xd lt nx and yd ge 0 and yd lt ny
h=(temporary(h)+1)*temporary(range)-1

; compute histogram
h=histogram(h,min=0,max=long(nx)*long(ny)-1l)
h=reform(h,nx,ny,/overwrite)  ; convert to 2D

; normalisation to maximum
mh=max(h)
h=float(h)/max(h)
print,'Density maximum: ',mh

RETURN,h

END