;***********************************************************************
;  This program creates a 2d histogram
;***********************************************************************

pro histogram_2d,x_data,y_data,x_start_bin,x_end_bin,x_dbin,x_bins,$
                 y_start_bin,y_end_bin,y_dbin,y_bins,data_freq,data_counts

;***********************************************************************
;  Data inputs
;***********************************************************************

;x_start_bin=0
;x_end_bin=25
;x_dbin=0.1  ;half bin width

;y_start_bin=0
;y_end_bin=25
;y_dbin=0.1  ;half bin width

;************************************************************************
;  Set up the bins
;************************************************************************

x_bins=x_start_bin
while x_bins[n_elements(x_bins)-1] lt x_end_bin $
  do x_bins=[x_bins,x_bins[n_elements(x_bins)-1]+(2.*x_dbin)]

y_bins=y_start_bin
while y_bins[n_elements(y_bins)-1] lt y_end_bin $
  do y_bins=[y_bins,y_bins[n_elements(y_bins)-1]+(2.*y_dbin)]

;************************************************************************
;  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)-1 do begin
  for j=0,n_elements(y_bins)-1 do begin
    result=where(x_data gt x_bins[i]-x_dbin and $
                 x_data le x_bins[i]+x_dbin and $
		             y_data gt y_bins[j]-y_dbin and $
		             y_data le y_bins[j]+y_dbin,count) 
    data_counts[i,j]=count
    ;print,x_bins[i]-x_dbin,x_bins[i]+x_dbin,y_bins[j]-y_dbin,y_bins[j]+y_dbin,count
    ;if count gt 0 then begin
    ;  print,x_bins[i]-x_dbin,x_bins[i]+x_dbin,y_bins[j]-y_dbin,y_bins[j]+y_dbin,count
    ;  print,x_data[result]
    ;  print,y_data[result]
    ;  stop
    ;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