;***********************************************************************
;  This program creates a 2d histogram in irregularly spaced bins
;***********************************************************************

pro histogram_2d_ireg_bins,x_data,y_data,x_bins,$
                           y_bins,data_freq,data_counts

;************************************************************************
;  INPUTS
;************************************************************************

;  Bins passed in
;x_bins
;y_bins

;************************************************************************
;  Determine if the bins are increasing or decreasing
;************************************************************************

if y_bins[1] gt y_bins[0] then yinc=1
if y_bins[1] lt y_bins[0] then yinc=-1
if x_bins[1] gt x_bins[0] then xinc=1
if x_bins[1] lt x_bins[0] then xinc=-1

;************************************************************************
;  Set up the 2D histogram array
;************************************************************************

data_counts=make_array(n_elements(x_bins),n_elements(y_bins),/float,value=0) 
 
;************************************************************************  
; Build the 2D histogram 
;************************************************************************

for i=0,n_elements(x_bins)-2 do begin
  for j=0,n_elements(y_bins)-2 do begin

    if yinc eq 1 and xinc eq 1 then begin
      result=where(x_data ge x_bins[i] and $
                   x_data lt x_bins[i+1] and $
	  	   y_data ge y_bins[j] and $
	  	   y_data le y_bins[j+1],count) 
    endif else if yinc eq -1 and xinc eq 1 then begin
       result=where(x_data ge x_bins[i] and $
                    x_data lt x_bins[i+1] and $
	  	    y_data le y_bins[j] and $
	  	    y_data gt y_bins[j+1],count) 
    endif else if yinc eq 1 and xinc eq -1 then begin
       result=where(x_data le x_bins[i] and $
                    x_data gt x_bins[i+1] and $
	  	    y_data ge y_bins[j] and $
	  	    y_data lt y_bins[j+1],count)  
    endif else if yinc eq -1 and xinc eq -1 then begin
       result=where(x_data ge x_bins[i] and $
                    x_data lt x_bins[i+1] and $
	  	    y_data ge y_bins[j] and $
	  	    y_data lt y_bins[j+1],count)  
    endif    

    data_counts[i,j]=count
    ;print,x_bins[i],x_bins[i+1],y_bins[j],y_bins[j+1],count
    if count gt 0 then begin
      ;print,x_data[result]
      ;print,y_data[result]
    endif
  endfor  
endfor  

;  Calculate the frequency array
if total(data_counts) gt 0 then begin
  data_freq=float(data_counts)/float(total(data_counts))
endif else begin
  data_freq=data_counts
endelse

end