;$Id: hscn_addrun.pro,v 1.3 2007/06/27 19:22:38 azuah Exp $
;###############################################################################
;
; NAME:
;  HSCN_ADDRUN
;
; PURPOSE:
;  Combines fixed-window scan data files together.
;
; CATEGORY:
;  DAVE, HFBS, 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.
;
;###############################################################################
;+
; NAME:
;       HSCN_ADDRUN
;
; PURPOSE:
;
;   Duplicate the functionality of the FORTRAN procedure named ADDRUN originally
;   written by Dan Neumann at the NIST Center for Neutron Research.  This function
;   is a very flexible means of combining data sets together through addition,
;   merging, and subtraction.  The quantity returned is a named structure whose
;   fields contain the result of the operations the user specified.  If the
;   algorithm failed for some reason (e.g. input error) then the function returns
;   -1 and the MSG keyword provides information in the form of string output as
;   to the reason of the failure.  Note that this is not an exact duplicate of
;   Neumann's ADDRUN function in syntax but in functionality.
;
; 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
;
; CATEGORY:
;
;       Data reduction, DAVE project
;
; CALLING SEQUENCE:
;
;   SOUT = ADDRUN( STRUCT_ARRAY,            $
;            MERGE = merge,            $
;            ADD = add,             $
;            POISSON = poisson,         $
;            INTERP = interp,          $
;            SUBTRACT = subtract,         $
;            BG_SCALE_FACTOR = bg_scale_factor, $
;            MSG = msg              )
;
; PARAMETERS
;
;   STRUCT_ARRAY: array of named structures, each one defined as follows:
;
;      ARRAY = {  DATASTRUCT,         $
;            PX:      PTR_NEW(),    $
;            PY:      PTR_NEW(),    $
;            PDY:    PTR_NEW(),    $
;            PMON:     PTR_NEW(),     $
;            PDMON:   PTR_NEW()     }
;
;      PX:   pointer to the independent variable
;      PY:   pointer to the dependent variable
;      PDY: pointer to the uncertainty in the dependent variable
;      PMON:   pointer to the monitor spectrum
;      PDMON:  pointer to the uncertainty in the monitor spectrum
;
; KEYWORDS
;
;   MERGE:      Keyword that, if set, merges each of the data sets together.  If there
;            are overlapping (identical) points in the domain then these points
;            are simply summed together in the returned data structure.
;   ADD:      Keyword that, if set, returns the sum of the data sets.
;   POISSON:    Keyword that specifies if the "ADDED" data are to be summed in
;            accordance with Poisson statistics.  The resulting values are
;            a "Poisson average" of the raw counts.  If not set then the data
;            are summed in a straight manner (Y = Y1 + Y2 + ...)
;   SUBTRACT:     Keyword that, if set, returns the difference between the
;            first data structure and second data structure.
;   INTERP:      Keyword that specifies if the data are to be interpolated
;            if they do not share identical values for the independent
;            variable.  A linear interpolation is used if set.  This option
;            is used for the ADD and SUBTRACT functionality.
;   BG_SCALE_FACTOR:  Scale factor which is multiplied with the background
;              data prior to the SUBTRACT operation (def: 1.0)
;   MSG:      Output keyword that provides a string array text message
;            indicating if the algorithm failed.
;
; COMMON BLOCKS:
;
;   None
;
; REQUIREMENTS:
;
;   DATASTRUCT__DEFINE.PRO
;   CMSET_OP.PRO
;
; MODIFICATION HISTORY:
;
;       Written by Rob Dimeo, July 8, 2003
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function hscn_addrun,   struct_array,            $
           merge = merge,            $
           add = add,               $
           poisson = poisson,          $
           interp = interp,           $
           subtract = subtract,        $
           bg_scale_factor = bg_scale_factor,   $
           msg = msg
!except = 0
msg = ""
if n_elements(struct_array) eq 0 then return,-1
if n_elements(interp) eq 0 then interp = 0
if n_elements(add) eq 0 then add = 0
if n_elements(subtract) eq 0 then subtract = 0
if n_elements(poisson) eq 0 then poisson = 0
if n_elements(merge) eq 0 then merge = 0
if n_elements(bg_scale_factor) eq 0 then bg_scale_factor = 1.0
n_struct = n_elements(struct_array)

if (subtract ne 0) or (add ne 0) then begin
; Perform the operations that are common to both.  These operations
; involve comparing the "sameness" of the independent variable

   same_array = intarr(n_struct)
   same_array[0] = 1
   same = 0
; Get out the minimum and maximum values for the independent variable and
; store these in XLO and XHI...
   xlo = -1.e55 & xhi = -xlo
   for i = 0,n_struct-1 do begin
     x = *(struct_array[i]).px
     if i gt 0 then begin
      if ((cmset_op(x,'AND',xold))[0] eq -1) and $
         (interp eq 0) then begin
         msg =   'At least two of the data sets do not have' + $
            ' the same independent variables.'
         return,-1
      endif
      if ((cmset_op(x,'AND',xold))[0] eq -1) then same = 0
      if ((cmset_op(x,'AND',xold))[0] ne -1) and $
         n_elements(cmset_op(x,'AND',xold)) eq n_elements(x) then same = 1
      same_array[i] = same
     endif
     xold = x
     xmin = min(x,max = xmax)
     xlo = xlo > xmin
     xhi = xhi < xmax
   endfor
; Use the first data set to specify the output independent variable
   x = *(struct_array[0]).px
   ok = where((x ge xlo) and (x le xhi),count_ok)
   if count_ok le 1 then begin
     msg = 'At least one of the data sets does not have enough'+$
         ' overlapping data points'
     return,-1
   endif
   x_out = x[ok]
   y_out = (*(struct_array[0]).py)[ok]
   dy_out = (*(struct_array[0]).pdy)[ok]
   mon_out = (*(struct_array[0]).pmon)[ok]
   dmon_out = (*(struct_array[0]).pdmon)[ok]
; Create the big arrays that will hold the data sets with the
; "common" independent variable
   x_array = fltarr(n_struct,count_ok)
   y_array = fltarr(n_struct,count_ok)
   dy_array = fltarr(n_struct,count_ok)
   mon_array = fltarr(n_struct,count_ok)
   dmon_array = fltarr(n_struct,count_ok)
   if (total(same_array) eq n_struct) then begin  ; all points overlap!
     for i = 0,n_struct-1 do begin
      x_array[i,*] = (*(struct_array[i]).px)[ok]
      y_array[i,*] = (*(struct_array[i]).py)[ok]
      dy_array[i,*] = (*(struct_array[i]).pdy)[ok]
      mon_array[i,*] = (*(struct_array[i]).pmon)[ok]
      dmon_array[i,*] = (*(struct_array[i]).pdmon)[ok]
     endfor
   endif else begin
     if (total(same_array) ne n_struct) and $
      interp eq 1 then begin
; We are assuming that the data sets are to be linearly interpolated
; at this stage and now do it...
      x_array[0,*] = (*(struct_array[0]).px)[ok]
      y_array[0,*] = (*(struct_array[0]).py)[ok]
      dy_array[0,*] = (*(struct_array[0]).pdy)[ok]
      mon_array[0,*] = (*(struct_array[0]).pmon)[ok]
      dmon_array[0,*] = (*(struct_array[0]).pdmon)[ok]
      for i = 1,n_struct-1 do begin
         xo = (*(struct_array[i]).px)
         yo = (*(struct_array[i]).py)
         dyo = (*(struct_array[i]).pdy)
         mono = (*(struct_array[i]).pmon)
         dmono = (*(struct_array[i]).pdmon)
         xi = x_out
         yi = interpol(yo,xo,xi)
         dyi = interpol(dyo,xo,xi)
         ymoni = interpol(mono,xo,xi)
         dymoni = interpol(dmono,xo,xi)
         x_array[i,*] = xi
         y_array[i,*] = yi
         dy_array[i,*] = dyi
         mon_array[i,*] = ymoni
         dmon_array[i,*] = dymoni
      endfor
     endif else begin
      msg = 'Points do not overlap.'
      return,-1
     endelse
   endelse
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  ADD            ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Consider the add case first
if add then begin

   if poisson eq 0 then begin
; STRAIGHT ADDITION
     xout = reform(x_array[0,*])
     yout = total(y_array,1)
     dy_out = sqrt(total(dy_array^2,1))
     mon_out = total(mon_array,1)
     dmon_out = sqrt(total(dmon_array^2,1))

   endif else begin ; Poisson combination of data
; POISSON ADDITION

     xout = reform(x_array[0,*])
; Calculate the resulting y
     ux = 1+bytarr(n_struct)
     sf = ux#(mon_array[0,*])/mon_array
     num = total((y_array/dy_array)^2,1)
     den = total((y_array*sf)/(dy_array*sf)^2,1)
     yout = num/den
; Calculate the resulting dy
     dy_out = sqrt(yout/den)
     mon_out = reform(mon_array[0,*])
     dmon_out = reform(dmon_array[0,*])

   endelse

endif ; end ADD

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  SUBTRACT        ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if subtract then begin
   if n_struct gt 2 then begin
     msg = 'Only two files can be specified at a time for subtraction'
     return,-1
   endif
   xout = reform(x_array[0,*])
   sf = mon_array[0,*]/mon_array[1,*]
   y1 = (y_array[0,*])
   y2 = sf*(y_array[1,*])
   dy1 = (dy_array[0,*])
   dy2 = sf*(dy_array[1,*])
   yout = y1-bg_scale_factor*y2
   dy_out = sqrt(dy1^2+(bg_scale_factor*dy2)^2)
   mon_out = mon_array[0,*]
   dmon_out = mon_array[1,*]
endif ; end SUBTRACT

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  MERGE         ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if merge then begin
; First append all of the data sets into common arrays
   x = *(struct_array[0]).px
   y = *(struct_array[0]).py
   dy = *(struct_array[0]).pdy
   mon = *(struct_array[0]).pmon
   dmon = *(struct_array[0]).pdmon
   for i = 1,n_struct-1 do begin
     x1 = *(struct_array[i]).px
     y1 = *(struct_array[i]).py
     dy1 = *(struct_array[i]).pdy
     mon1 = *(struct_array[i]).pmon
     dmon1 = *(struct_array[i]).pdmon
; Are there any overlapping points? If so, add them together.
     overlap_values = cmset_op(x,'AND',x1,count = count_overlap)
     if count_overlap gt 0 then begin
         ;; RTA - logical error present when x has one element
         ;;       only. Add code to deal with this separately.
         nx = n_elements(x)
         if (nx eq 1) then begin
             y = y + y1
             mon = mon + mon1
             dy = sqrt(dy^2 + dy1^2)
         endif else begin
             for j = 0,count_overlap-1 do begin
                 ;; Use the rule for combining like quantities (i.e. combination by
                 ;; "least-variance" for Gaussian deviates)
                 x_index = where(x eq overlap_values[j])
                 x1_index = where(x1 eq overlap_values[j])
                 sf = mon[x_index[0]]/mon1[x1_index[0]]
                 y[x_index[0]] = y[x_index[0]] + y1[x1_index[0]]
                 mon[x_index[0]] = mon[x_index[0]]+mon1[x1_index[0]]
                 dy[x_index[0]] = sqrt((dy[x_index[0]])^2 + (dy1[x1_index[0]])^2)
             endfor
         endelse
     endif
; Now put the points in x1 that do not overlap x into x
     exclusive_values = cmset_op(x,'XOR',x1,count = count_exclusive)
     if count_exclusive gt 0 then begin
      for j = 0,count_exclusive-1 do begin
         x1_index = where(x1 eq exclusive_values[j],count_ok)
         if count_ok then begin
           x = [x,x1[x1_index[0]]]
           y = [y,y1[x1_index[0]]]
           dy = [dy,dy1[x1_index[0]]]
           mon = [mon,mon1[x1_index[0]]]
           dmon = [dmon,dmon1[x1_index[0]]]
         endif
      endfor
     endif
   endfor
   xout = x
   yout = y
   dy_out = dy
   mon_out = mon
   dmon_out = dmon
endif ; end MERGE


; Fill up the outgoing structure with the appropriate values
struct   = {datastruct}
struct.px = ptr_new(xout)
struct.py = ptr_new(yout)
struct.pdy = ptr_new(dy_out)
struct.pmon = ptr_new(mon_out)
struct.pdmon = ptr_new(dmon_out)

return,struct
end
