; $Id: merge_arrays.pro,v 1.3 2005/02/09 18:13:38 dimeo Exp $
;###############################################################################
;
; NAME:
;  MERGE_ARRAYS
;
; PURPOSE:
;  Means to merge two arrays with (possibly) overlapping independent variables.
;
; CATEGORY:
;  DAVE, TAS, Data Reduction
;
; AUTHOR:
;   Robert M. Dimeo, Ph.D.
;   NIST Center for Neutron Research
;   100 Bureau Drive
;   Gaithersburg, MD 20899
;   Phone: (301) 975-8135
;   E-mail: robert.dimeo@nist.gov
;   http://www.ncnr.nist.gov/staff/dimeo
;
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;
;;;;;;;;;;;;;;
function merge_arrays,x,y,dy,uniq_index = uniq_index,straight = straight
if n_params() ne 3 then return,-1
xsort = sort(x)
x = x[xsort]
y = y[xsort] & dy = dy[xsort]
var_y = dy^2
uniq_index = uniq(x)
x_uniq = x[uniq_index]
y_uniq = y[uniq_index]
dy_uniq = dy[uniq_index]

nu = n_elements(uniq_index)
nx = n_elements(x)
; Calculate an array of the same length as the unique array whose zero
; positions indicate the occurrence of a multiply-occurring value.
du = [uniq_index[0]+1,uniq_index[1:nu-1] - uniq_index[0:nu-2]]
where_multiple = where(du ne 1,count_multiple)
overlap = 0L

if n_elements(straight) eq 0 then begin
; Sum using Poisson statistics
; If there are overlapping data points then add them using Poisson addition.
   if count_multiple gt 0 then begin
      for i = 0,count_multiple-1 do begin
         x_same = where(x eq (x_uniq)[where_multiple[i]])
         y_num = total((y[x_same]/dy[x_same])^2)
         y_den = total(y[x_same]/(dy[x_same])^2)
         y_bar = y_num/y_den
         var_y_bar = y_bar/total(y[x_same]/dy[x_same]^2)
         y_uniq[where_multiple[i]] = y_bar
         dy_uniq[where_multiple[i]] = sqrt(var_y_bar)
         overlap = overlap + n_elements(x_same)
      endfor
   endif
endif else begin
; Use straight addition
   if count_multiple gt 0 then begin
      for i = 0,count_multiple-1 do begin
         x_same = where(x eq (x_uniq)[where_multiple[i]])
         y_sum = total(y[x_same])
         y_uniq[where_multiple[i]] = y_sum
         dy_uniq[where_multiple[i]] = sqrt(y_sum)
         overlap = overlap + n_elements(x_same)
      endfor
   endif
endelse
print,overlap
return,{x:x_uniq,y:y_uniq,dy:dy_uniq}
end
;;;;;;;;;;
function scat_fun,x
area = 500.0 & center = 0.0 & fwhm = 1.0 & bg = 0.1*area
sig = fwhm/2.354
yp = (area/sqrt(2.0*!dpi*sig^2))*exp(-0.5*((x-center)/sig)^2) + bg
y = fltarr(n_elements(x))
for i = 0,n_elements(x)-1 do begin
   y[i] = randomn(s,1,poisson = yp[i])
   dum = temporary(s)
endfor
return,y
end
;;;;;;;;;;
pro test_merge
!except = 0
xlo = -5.0 & xhi = 5.0 & nx = 25
dx = (xhi-xlo)/(nx-1.0)
x1 = xlo+dx*findgen(nx)
x2 = x1+0.5*dx
x1 = [x1,0.0]
x2 = [x2,0.0]
y1 = scat_fun(x1) & y2 = scat_fun(x2)
dy1 = sqrt(y1) & dy2 = sqrt(y2)
xstr = merge_arrays([x1,x2],[y1,y2],[dy1,dy2])
x = xstr.x & y = xstr.y & dy = xstr.dy


device,decomposed = 0
loadct,0,/silent
white = 255 & black = 0
window,0
plot,x1,y1,psym = 4,xtitle = '!3E (meV)',ytitle = 'Intensity',title = 'Original Data Sets'
errplot,x1,y1-dy1,y1+dy1,width = 0.0
;window,1
oplot,x2,y2,psym = 5
errplot,x2,y2-dy2,y2+dy2,width = 0.0
window,2
plot,x,y,psym = 4,xtitle = '!3E (meV)',ytitle = 'Intensity',title = '"Summed" Data Sets
errplot,x,y-dy,y+dy,width = 0.0
end