; $Id: fans_output.pro,v 1.37 2007/10/18 20:55:18 azuah Exp $
;###############################################################################
;
; NAME:
;   fans_output
;
; PURPOSE:
;   Utility routines for FANS data reduction application defined in
;   the file fans_datreduc.pro. The routines in this file perform
;   output data calculations and saving to file in DAVE format.
;
; CATEGORY:
;   FANS Data Reduction
;
; AUTHOR:
;   Richard Tumanjong Azuah
;   NIST Center for Neutron Research
;   100 Bureau Drive, Gaithersburg, MD 20899
;   United States
;   azuah@nist.gov; (301) 9755604
;   June, 2002
;
; 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.
;
;###############################################################################



;===============================================================================
; fans_apdData
;
; PURPOSE:
;   Append sample, bkgd or fast bkgd datasets subject to the
;   requirement that the energy range of the data DO NOT COMPLETELY
;   OVERLAP. Any overlapping points in the higher energy-range data is
;   ignored.
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
;   dataType [in] - data type string: 'samp'|'bkgd'|'fast'
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_apdData, sPtr, dataType, err=err, emsg=emsg
err=0  & emsg=''

; get appropriate data; samp, bkgd or fast bkgd
switch dataType of

    'samp': begin
        nd = (*sPtr).sampCount
        if (nd le 1) then return
        dataPtr = (*sPtr).sDataPtr
        typeName = 'Sample'
        
        break
    end

    'bkgd': begin
        nd = (*sPtr).bkgdCount
        if (nd le 1) then return
        dataPtr = (*sPtr).bDataPtr
        typeName = 'Background'

        break
    end

    'fast': begin
        nd = (*sPtr).fastCount
        if (nd le 1) then return
        dataPtr = (*sPtr).fDataPtr
        typeName = 'Fast bkgd'

        break
    end

endswitch

; Get the min and max energy values of the dataset
minima = fltarr(nd) - 1
maxima = minima
for i = 0,nd-1 do begin
    mine = min( (*(*((*dataPtr)[i])).energyPtr), max=maxe)
    minima[i] = mine
    maxima[i] = maxe
endfor 
ord = sort(minima)              ; find order of increasing min energy


; make a copy of the dataset in the list with the smallest min
; energy. 
mergePtr = fans_initDataPtr((*dataPtr)[ord[0]])
e1Ptr = (*mergePtr).energyPtr  
mon1 = float((*mergePtr).monitor * (*mergePtr).prefactor)
;oldHist = (*(*mergePtr).historyPtr)
(*(*mergePtr).historyPtr)[0] = typeName+' data is a merger of the following file(s):'

; loop through and merge the datasets.
min1 = minima[ord[0]]  &  max1 = maxima[ord[0]]
for i = 1,nd-1 do begin
    j = ord[i]

    ;; ignore cases of complete overlap
    min2 = minima[j]  &  max2 = maxima[j]
    if ( (min1 le min2)  and (max1 ge max2) ) then continue ; total overlap - goto next data
    
    ;; identify partially overlapping regions. Only the upper bounds
    ;; of the overlap is important since, for FANS, it was decided
    ;; that the overlapping data for the higher energy scan should
    ;; simply be discarded.
    curPtr = (*dataPtr)[j]
    e2Ptr = (*curPtr).energyPtr
    ii = 0
    res = where( *e2Ptr gt max1)
    if (res[0] ge 0) then ii = res[0] ; ii is the *e2Ptr index of upper bound of overlap
    
    ;; now copy non-overlapping data from the higher energy scan to
    ;; the merged data
    ;; NOTE: the second data is scale to the monitor counts of the first
    ndets = (*mergePtr).ndets
    nspec1 = (*mergePtr).nspec
    nspec2 = (*curPtr).nspec
    mon2 = float((*mergePtr).monitor * (*mergePtr).prefactor)
    scale = mon1/mon2
    (*mergePtr).eFinal = max2
    (*mergePtr).nspec = (*mergePtr).nspec + nspec2 - ii
    nspec = (*mergePtr).nspec
    (*(*mergePtr).qxPtr) = [(*(*mergePtr).qxPtr),(*(*curPtr).qxPtr)[ii:nspec2-1]]
    (*(*mergePtr).qyPtr) = [(*(*mergePtr).qyPtr),(*(*curPtr).qyPtr)[ii:nspec2-1]]
    (*(*mergePtr).qzPtr) = [(*(*mergePtr).qzPtr),(*(*curPtr).qzPtr)[ii:nspec2-1]]
    (*(*mergePtr).tempPtr) = [(*(*mergePtr).tempPtr),(*(*curPtr).tempPtr)[ii:nspec2-1]]
    (*(*mergePtr).energyPtr) = [(*(*mergePtr).energyPtr),(*(*curPtr).energyPtr)[ii:nspec2-1]]
    (*(*mergePtr).timePtr) = [(*(*mergePtr).timePtr),(*(*curPtr).timePtr)[ii:nspec2-1] * scale]
    (*(*mergePtr).totCntsPtr) = [(*(*mergePtr).totCntsPtr),(*(*curPtr).totCntsPtr)[ii:nspec2-1] * scale]
    (*(*mergePtr).totErrsPtr) = [(*(*mergePtr).totErrsPtr),(*(*curPtr).totErrsPtr)[ii:nspec2-1] * scale]
    detCnts = fltarr(ndets,nspec)
    detErrs = fltarr(ndets,nspec)
    detCnts[0:ndets-1,0:nspec1-1] = (*(*mergePtr).detCntsPtr)
    detErrs[0:ndets-1,0:nspec1-1] = (*(*mergePtr).detErrsPtr)
    detCnts[0:ndets-1,nspec1:nspec-1] = (*(*curPtr).detCntsPtr)[*,ii:nspec2-1] * scale
    detErrs[0:ndets-1,nspec1:nspec-1] = (*(*curPtr).detErrsPtr)[*,ii:nspec2-1] * scale
    (*(*mergePtr).detCntsPtr) = detCnts
    (*(*mergePtr).detErrsPtr) = detErrs
    ;;(*(*mergePtr).detCntsPtr) = [(*(*mergePtr).detCntsPtr),(*(*curPtr).detCntsPtr)[*,ii:nspec2-1]]
    max1 = max2

    ;; update history pointer with merged data filename (index 1 of
    ;; the dataPtr.historyPtr)
    (*(*mergePtr).historyPtr) = [(*(*mergePtr).historyPtr),(*(*curPtr).historyPtr)[1]]
endfor

; Add merged data to list of data. Update info for appropriate data type
if (ptr_valid(mergePtr) eq 0) then return     

*dataPtr = [*dataPtr,mergePtr]  
switch dataType of 

    'samp': begin
        (*sPtr).sampIndex = nd
;        (*sPtr).sampCount = nd + 1
        (*sPtr).sampStatus = 2
        break
    end

    'bkgd': begin
        (*sPtr).bkgdIndex = nd
;        (*sPtr).bkgdCount = nd + 1
        (*sPtr).bkgdStatus = 2
        break
    end 

    'fast': begin
        (*sPtr).fastIndex = nd
;        (*sPtr).fastCount = nd + 1
        (*sPtr).fastStatus = 2
        break
    end 

endswitch 

end
;-------------------------------------------------------------------------------



;===============================================================================
; fans_sumData
;
; PURPOSE:
;   Calculates the sum or weighted average of all the current sample,
;   bkgd or fast bkgd data.
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
;   dataType [in] - data type string: 'samp'|'bkgd'|'fast'
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_sumData, sPtr, dataType, err=err, emsg=emsg
err=0  & emsg=''

; get appropriate data; samp or bkgd
switch dataType of

    'samp': begin
        nd = (*sPtr).sampCount
        if (nd le 1) then return
        dataPtr = (*sPtr).sDataPtr
        fnames = (*(*sPtr).sampListPtr)
        typeName = 'Sample '
        id = 0
        break
    end

    'bkgd': begin
        nd = (*sPtr).bkgdCount
        if (nd le 1) then return
        dataPtr = (*sPtr).bDataPtr
        fnames = (*(*sPtr).bkgdListPtr)
        typeName = 'Background '
        id = 1
        break
    end

    'fast': begin
        nd = (*sPtr).fastCount
        if (nd le 1) then return
        dataPtr = (*sPtr).fDataPtr
        fnames = (*(*sPtr).fastListPtr)
        typeName = 'Fast bkgd '
        id = 2
        break
    end

endswitch

; Get the min and max energy values of the dataset
minima = fltarr(nd) - 1
maxima = minima
for i = 0,nd-1 do begin
    mine = min( (*(*((*dataPtr)[i])).energyPtr), max=maxe)
    minima[i] = mine
    maxima[i] = maxe
endfor 
ord = sort(minima)              ; find order of increasing min energy

; make a copy of the first dataset in the list. Note that dataPtr
; points to an array of ptrs (each pointing to a dataset).
sumPtr = fans_initDataPtr((*dataPtr)[ord[0]])
monConstant = (*sumPtr).monitor * (*sumPtr).prefactor ; record that for first dataset
omega1 = (*(*sumPtr).energyPtr)
temp1 =  (*(*sumPtr).tempPtr)
(*(*sumPtr).historyPtr)[0] = typeName+'data is a sum of the following file(s):'

; Monitor requires special attention! The monitor count is made of a
; product of the mon cnt and the prefactor and it is not possible to
; keep track of these separately when files are averaged. So set monitor
; count to mon*pref and pref to 1. Subsequent data are scaled to the
; mon*pref value of the first dataset.
;(*sumPtr).monitor = (*sumPtr).monitor * (*sumPtr).prefactor
;(*sumPtr).prefactor = 1

; loop through and sum all the datasets
; only the time, energy, detCnts, totCnts and totErrs components need
; to be averaged! - this is done in the inner loop
min1 = minima[ord[0]]  &  max1 = maxima[ord[0]]
for i = 1,nd-1 do begin         
    j = ord[i]

    ;; if dataset barely overlap or not at all then skip the current
    ;; dataset and advise user to select the merge option instead
    min2 = minima[j]  &  max2 = maxima[j]
    ;; get the current dataset to be dealt with
    curPtr = (*dataPtr)[j]
    omega2 = (*(*curPtr).energyPtr)
    temp2 = (*(*curPtr).tempPtr)
    size1 = n_elements(omega1)
    size2 = n_elements(omega2)
;    mon1 = float((*sumPtr).monitor * (*sumPtr).prefactor)
;    mon2 = float((*curPtr).monitor * (*curPtr).prefactor)
;    mon_scale = mon1/mon2
;    mon_scale_sqrt = sqrt(mon_scale)
    monitor1 = (*(*sumPtr).monPtr)
    monitor2 = (*(*curPtr).monPtr)
    
    ;; Examine various exclusive scenarios

    ;---------------------------------------------------------------
    ;; Scenario 1 - min2 > min1 ==> there is a 'no overlap' region in
    ;; dataset1. Have to copy this data into a buffer
    ibeg1 = 0 & iend1=0 & ibeg2=0 & iend2=0
    if (min1 ne min2) then begin
        res = where(omega1 ge min2) ; ge is deliberate
        iend1 = (res[0] le 0)? size1-1 : res[0] - 1
        if (iend1 ge 0) then begin
            omega = omega1[ibeg1:iend1]
            times = (*(*sumPtr).timePtr)[ibeg1:iend1]
            temperature = (*(*sumPtr).tempPtr)[ibeg1:iend1]
            totCnts = (*(*sumPtr).totCntsPtr)[ibeg1:iend1]
;            totErrs = (*(*sumPtr).totErrsPtr)[ibeg1:iend1]
            detCnts = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1]
;            detErrs = (*(*sumPtr).detErrsPtr)[*,ibeg1:iend1]
            monitor = monitor1[ibeg1:iend1]
            ibeg1=iend1+1
            iend1 = iend1 + 1
        endif
    endif
    
    ;----------------------------------------------------------------------
    ;; Scenario 2 - the overlap region;
    if (max1 lt min2) then goto, scenario3 ; inadequate overlap
    res = where(omega1 gt max2) ; the 'gt' is deliberate!
    iend1 = (res[0] lt 0)? size1-1 : res[0]-1 ; upper bound of overlap for ds1
    if (iend1 lt ibeg1) then iend1 = ibeg1
    ibeg2 = 0                   ; lower bound of overlap for ds2
    res = where(omega2 gt max1) 
    iend2 = (res[0] lt 0)? size2-1 : res[0]-1 ; upper bound of overlap for ds2
    if (iend2 lt ibeg2) then iend2 = ibeg2
    if ((where(omega1[ibeg1:iend1] ne omega2[ibeg2:iend2]))[0] ge 0) then begin 
        ;; ==> in overlap region, omega1 != omega2
        ;; hence need to merge (interleave) the data

        ;; determine which of dataset has a smaller min omega and use
        ;; that ds as base for inserting or appending the other ds 
        if (omega1[ibeg1] le omega2[ibeg2]) then begin
            om = omega1[ibeg1:iend1]
            mn = monitor1[ibeg1:iend1]
            temp = temp1[ibeg1:iend1]
            t = (*(*sumPtr).timePtr)[ibeg1:iend1]
            x = (*(*sumPtr).totCntsPtr)[ibeg1:iend1]
;            w = (*(*sumPtr).totErrsPtr)[ibeg1:iend1]
            xd = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1]
;            wd = (*(*sumPtr).detErrsPtr)[*,ibeg1:iend1]
            om2 = omega2[ibeg2:iend2]
            mn2 = monitor2[ibeg2:iend2]
            tmp2 = temp2[ibeg2:iend2]
            t2 = (*(*curPtr).timePtr)[ibeg2:iend2] ;* mon_scale
            x2 = (*(*curPtr).totCntsPtr)[ibeg2:iend2] ;* mon_scale
;            w2 = (*(*curPtr).totErrsPtr)[ibeg2:iend2] * mon_scale ;mon_scale_sqrt
            xd2 = (*(*curPtr).detCntsPtr)[*,ibeg2:iend2] ;* mon_scale
;            wd2 = (*(*curPtr).detErrsPtr)[*,ibeg2:iend2] * mon_scale ;mon_scale_sqrt
            ibeg = 0
            iend = iend2 - ibeg2
            len = n_elements(om) - 1
            first_dset = 1 
        endif else begin
            om = omega2[ibeg2:iend2]
            mn = monitor2[ibeg2:iend2]
            temp = temp2[ibeg2:iend2]
            t = (*(*curPtr).timePtr)[ibeg2:iend2] ;* mon_scale
            x = (*(*curPtr).totCntsPtr)[ibeg2:iend2] ;* mon_scale
;            w = (*(*curPtr).totErrsPtr)[ibeg2:iend2] * mon_scale ;mon_scale_sqrt
            xd = (*(*curPtr).detCntsPtr)[*,ibeg2:iend2] ;* mon_scale
;            wd = (*(*curPtr).detErrsPtr)[*,ibeg2:iend2] * mon_scale ;mon_scale_sqrt
            om2 = omega1[ibeg1:iend1]
            mn2 = monitor1[ibeg1:iend1]
            tmp2 = temp1[ibeg1:iend1]
            t2 = (*(*sumPtr).timePtr)[ibeg1:iend1]
            x2 = (*(*sumPtr).totCntsPtr)[ibeg1:iend1]
;            w2 = (*(*sumPtr).totErrsPtr)[ibeg1:iend1]
            xd2 = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1]
;            wd2 = (*(*sumPtr).detErrsPtr)[*,ibeg1:iend1]
            ibeg = 0
            iend = iend1 - ibeg1
            len = n_elements(om) - 1 ;iend2 - ibeg2
            first_dset = 2
        endelse
        for k=ibeg,iend do begin
            ;; loop thro and locate pos of om2[k] in om

            res = where(om ge om2[k])
            if ( res[0] ge 0 ) then begin 
                if (om2[k] eq om[res[0]]) then begin
                    ;; ==> energy value of om2[k] matches with one in om
                    ;; therefore perform an average.
                    kk = res[0]
                    x3 = x[kk]
                    x4 = x2[k]
                    t3 = t[kk]
                    t4 = t2[k]
                    xd3 = xd[*,kk]
                    xd4 = xd2[*,k]

                    ;; 1) Average total counts
                    x[kk] = x[kk] + x2[k]
                    mn[kk] = mn[kk] + mn2[k]
                    ;w[kk] = 0.5 * sqrt(sum)
                    ;; 2) Average counting times
                    t[kk] = t[kk] + t2[k]
                    ;; 3) Average detector counts
                    xd[*,kk] = xd[*,kk] + xd2[*,k]
                    ;wd[*,kk] = 0.5 * sqrt(sum)

;                    ;; 1) Average total counts
;                    w3sq = (w[kk])^2
;                    w4sq = (w2[k])^2
;                    invwsq = (1.0/w3sq) + (1.0/w4sq)
;                    x[kk] = ((x3/w3sq) + (x4/w4sq)) / invwsq
;                    w[kk] = sqrt(1/invwsq)
;                    ;; 2) Average counting times
;                    t[kk]  = ((t3/w3sq) + (t4/w4sq)) / invwsq
;                    ;; 3) Average detector counts
;                    wd3sq = (wd[*,kk])^2
;                    wd4sq = (wd2[*,k])^2
;                    invwdsq = (1.0/wd3sq) + (1.0/wd4sq) 
;                    xd[*,kk] = ((xd3/wd3sq) + (xd4/wd4sq)) / invwdsq
;                    wd[*,kk] = sqrt(1/invwdsq) 


                    continue    ; move on to next for loop iteration
                endif
            endif  
            sz = size(xd)
            nw = sz[1]  
            ns = (sz[0] eq 1)? 1 : sz[2]
            tmpxd = fltarr(nw,ns+1)
            tmpwd = fltarr(nw,ns+1)
            if (res[0] lt 0) then begin
                ;; ==> need to _append_ ds2 values at end
                om = [om,om2[k]]
                temp = [temp,tmp2[k]]
                t = [t,t2[k]]
                x = [x,x2[k]]
                mn = [mn,mn2[k]]
;                w = [w,w2[k]]
                tmpxd[*,0:ns-1] = xd
;                tmpwd[*,0:ns-1] = wd
                tmpxd[*,ns] = xd2[*,k]
;                tmpwd[*,ns] = wd2[*,k]
                xd = tmpxd
                wd = tmpwd
            endif else begin
                ;; ==> need to _insert_ ds2 vals before res[0] index
                ;; in ds1
                kk = res[0]
                om = [om[0:kk-1],om2[k],om[kk:len]]
                temp = [temp[0:kk-1],tmp2[k],temp[kk:len]]
                t = [t[0:kk-1],t2[k],t[kk:len]]
                x = [x[0:kk-1],x2[k],x[kk:len]]
                mn = [mn[0:kk-1],mn2[k],mn[kk:len]]
;                w = [w[0:kk-1],w2[k],w[kk:len]]
                tmpxd[*,0:kk-1] = xd[*,0:kk-1]
;                tmpwd[*,0:kk-1] = wd[*,0:kk-1]
                tmpxd[*,kk] = xd2[*,k]
;                tmpwd[*,kk] = wd2[*,k]
                tmpxd[*,kk+1:ns] = xd[*,kk:len]
;                tmpwd[*,kk+1:ns] = wd[*,kk:len]
                xd = tmpxd
;                wd = tmpwd
            endelse
            len = n_elements(om) - 1 ; update len to reflect expanded buffers
        endfor
        iend1 = iend1 + 1
        iend2 = iend2 + 1
    endif else begin 
        ;; ==> in overlap region, omega1 == omega2
        ;; hence do the averaging in the overlap region
        om = omega1[ibeg1:iend1]
        mn = monitor1[ibeg1:iend1]
        mn2 = monitor2[ibeg2:iend2]
        temp = temp1[ibeg1:iend1]

        ;; 1) Average total counts
        buf1 = (*(*sumPtr).totCntsPtr)[ibeg1:iend1]
        buf2 = (*(*curPtr).totCntsPtr)[ibeg2:iend2]
        x =  buf1 + buf2
        mn = mn + mn2
        
        ;; 2) Average counting times
        buf1 = (*(*sumPtr).timePtr)[ibeg1:iend1]
        buf2 = (*(*curPtr).timePtr)[ibeg2:iend2]; * mon_scale
        t = buf1 + buf2

        ;; 3) Average detector counts
        buf1 = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1]
        buf2 = (*(*curPtr).detCntsPtr)[*,ibeg2:iend2]; * mon_scale
        xd = buf1 + buf2
        ;wd = 0.5 * sqrt(sum)

;        ;; 1) Average total counts
;        x1 = (*(*sumPtr).totCntsPtr)[ibeg1:iend1]
;        x2 = (*(*curPtr).totCntsPtr)[ibeg2:iend2] * mon_scale
;        w1sq = ((*(*sumPtr).totErrsPtr)[ibeg1:iend1])^2
;        w2sq = ((*(*curPtr).totErrsPtr)[ibeg2:iend2])^2 * (mon_scale)^2
;        invwsq = (1.0/w1sq) + (1.0/w2sq)
;        x = ((x1/w1sq) + (x2/w2sq)) / invwsq
;        w = sqrt(1/invwsq)
;        ;; 2) Average counting times
;        t1 = (*(*sumPtr).timePtr)[ibeg1:iend1]
;        t2 = (*(*curPtr).timePtr)[ibeg2:iend2] * mon_scale
;        t  = ((t1/w1sq) + (t2/w2sq)) / invwsq
;        ;; 3) Average detector counts
;        xd1 = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1]
;        xd2 = (*(*curPtr).detCntsPtr)[*,ibeg2:iend2] * mon_scale
;        wd1sq = ((*(*sumPtr).detErrsPtr)[*,ibeg1:iend1])^2
;        wd2sq = ((*(*curPtr).detErrsPtr)[*,ibeg2:iend2])^2 * (mon_scale)^2
;        invwdsq = (1.0/wd1sq) + (1.0/wd2sq) 
;        xd = ((xd1/wd1sq) + (xd2/wd2sq)) / invwdsq
;        wd = sqrt(1/invwdsq)

        iend2 = iend2 + 1
        iend1 = iend1 + 1
    endelse

    ;; include the om,t,x,w,xd,wd calculated using either interleaving or
    ;; averaging into the sum buffers
    if (ibeg1 eq 0) then begin
        omega = om
        monitor = mn
        temperature = temp
        times = t
        totCnts = x
;        totErrs = w
        detCnts = xd
;        detErrs = wd
    endif else begin
        ;; sum buffers already contain data from Scenario 1
        omega = [omega,om]
        monitor = [monitor,mn]
        temperature = [temperature,temp]
        times = [times,t]
        totCnts = [totCnts,x]
;        totErrs = [totErrs,w]
        if ((iend1 eq ibeg1) and (n_elements(om) eq 1)) then begin  
            ;; ie only a one data point overlap ==> xd,wd are 1D not 2D!
            sz = size(detCnts)
            nw = sz[1]
            ns = (sz[0] eq 1)? 1 : sz[2]
            tmpdetCnts = fltarr(nw,ns+1)
            tmpdetErrs = fltarr(nw,ns+1)
            tmpdetCnts[*,0:ns-1] = detCnts
;            tmpdetErrs[*,0:ns-1] = detErrs
            tmpdetCnts[*,ns] = xd
;            tmpdetErrs[*,ns] = wd
            detCnts = tmpdetCnts
;            detErrs = tmpdetErrs
        endif else begin
            sz1 = size(detCnts)
            sz2 = size(xd)
            nw = sz1[1]
            ns1 = (sz1[0] eq 1)? 1 : sz1[2] 
            ns = (sz2[0] eq 1)? 1 : sz2[2]
            ns = ns + ns1
            tmpdetCnts = fltarr(nw,ns)
;            tmpdetErrs = fltarr(nw,ns)
            tmpdetCnts[*,0:ns1-1] = detCnts
;            tmpdetErrs[*,0:ns1-1] = detErrs
            tmpdetCnts[*,ns1:ns-1] = xd
;            tmpdetErrs[*,ns1:ns-1] = wd
            detCnts = tmpdetCnts
;            detErrs = tmpdetErrs
        endelse
    endelse
 

    ;----------------------------------------------------------------------
    ;; Scenario 3 - beyond the overlap region - include any remaining data
    ;; Only ds1 OR ds2 could still contain unaccounted data; NOT
    ;; BOTH. However, both could be all used up already
    scenario3:
;    if (((iend1 eq 0) and (iend1 lt size1-1)) or (iend1 lt size1-1))then begin ; pick up the remaining data in ds1
    if ((iend1 eq 0) or (iend1 le size1-1)) then begin ; pick up the remaining data in ds1
        res = where(omega1[iend1] le omega, count)
        ibeg1 = (count lt 1)? iend1 : iend1 + 1
        iend1 = size1 - 1
        omega = [omega,omega1[ibeg1:iend1]]
        monitor = [monitor,monitor1[ibeg1:iend1]]
        temperature = [temperature,temp1[ibeg1:iend1]]
        times = [times,(*(*sumPtr).timePtr)[ibeg1:iend1]]
        totCnts = [totCnts,(*(*sumPtr).totCntsPtr)[ibeg1:iend1]]
;        totErrs = [totErrs,(*(*sumPtr).totErrsPtr)[ibeg1:iend1]]
        sz1 = size(detCnts)
        sz2 = size((*(*sumPtr).detCntsPtr)[*,ibeg1:iend1])
        nw = sz1[1]
        ns1 = (sz1[0] eq 1)? 1 : sz1[2] 
        ns = (sz2[0] eq 1)? 1 : sz2[2]
        ns = ns + ns1
        tmpdetCnts = fltarr(nw,ns)
;        tmpdetErrs = fltarr(nw,ns)
        tmpdetCnts[*,0:ns1-1] = detCnts
;        tmpdetErrs[*,0:ns1-1] = detErrs
        tmpdetCnts[*,ns1:ns-1] = (*(*sumPtr).detCntsPtr)[*,ibeg1:iend1] 
;        tmpdetErrs[*,ns1:ns-1] = (*(*sumPtr).detErrsPtr)[*,ibeg1:iend1]
        detCnts = tmpdetCnts
;        detErrs = tmpdetErrs
        iend1 = iend1 + 1
    endif else if ((iend2 eq 0) or (iend2 le size2-1)) then begin ; pick up the remaining data in ds2
        res = where(omega2[iend2] le omega, count)
        ibeg2 = (count lt 1)? iend2 : iend2 + 1
        iend2 = size2 - 1
        omega = [omega,omega2[ibeg2:iend2]]
        monitor = [monitor,monitor2[ibeg2:iend2]]
        temperature = [temperature,temp2[ibeg2:iend2]]
        times = [times,(*(*curPtr).timePtr)[ibeg2:iend2]] 
        totCnts = [totCnts,(*(*curPtr).totCntsPtr)[ibeg2:iend2]]
;        totErrs = [totErrs,(*(*curPtr).totErrsPtr)[ibeg2:iend2]*mon_scale]
        sz1 = size(detCnts)
        sz2 = size((*(*curPtr).detCntsPtr)[*,ibeg2:iend2])
        nw = sz1[1]
        ns1 = (sz1[0] eq 1)? 1 : sz1[2] 
        ns = (sz2[0] eq 1)? 1 : sz2[2]
        ns = ns + ns1
        tmpdetCnts = fltarr(nw,ns)
;        tmpdetErrs = fltarr(nw,ns)
        tmpdetCnts[*,0:ns1-1] = detCnts
;        tmpdetErrs[*,0:ns1-1] = detErrs
        tmpdetCnts[*,ns1:ns-1] = (*(*curPtr).detCntsPtr)[*,ibeg2:iend2]
;        tmpdetErrs[*,ns1:ns-1] = (*(*curPtr).detErrsPtr)[*,ibeg2:iend2]*mon_scale
        detCnts = tmpdetCnts
;        detErrs = tmpdetErrs
        iend2 = iend2 + 1
    endif 

    ;----------------------------------------------------------
    ;; Phew! Now, should have reached the end of both datasets
    ;; If not, then there is something wrong so inform user but do nothing
    if ( (iend1 le size1-1) or (iend2 le size2-1) ) then begin
        message = ['Processing data from: ', $
                   fnames[j], $
                   'It apeears not all the data was accounted for!']
        void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/information)
    endif


    ;;------------------------------------------------------------
    ;; replace sum data with calculations
    (*(*sumPtr).energyPtr) = omega
    (*(*sumPtr).timePtr) = times
    (*(*sumPtr).tempPtr) = temperature
    (*(*sumPtr).totCntsPtr) = totCnts
;    (*(*sumPtr).totErrsPtr) = totErrs
    (*(*sumPtr).detCntsPtr) = detCnts
;    (*(*sumPtr).detErrsPtr) = detErrs
    (*(*sumPtr).monPtr) = monitor
    (*sumPtr).nspec = n_elements(omega)
  
    ;; update historyPtr
    (*(*sumPtr).historyPtr) = [(*(*sumPtr).historyPtr),(*(*curPtr).historyPtr)[1]]

    ;; update certain quantities before the next loop
    omega1 = (*(*sumPtr).energyPtr)
    temp1 = (*(*sumPtr).tempPtr)
    min1 = min(omega1, max=max1)
endfor

; normalize spectra to a constant monitor
ndet = (size((*(*sumPtr).detCntsPtr)))[1]
factor = monConstant / (*(*sumPtr).monPtr)
(*(*sumPtr).totErrsPtr) = sqrt((*(*sumPtr).totCntsPtr)) * factor
(*(*sumPtr).totCntsPtr) = (*(*sumPtr).totCntsPtr) * factor
(*(*sumPtr).detErrsPtr) = sqrt((*(*sumPtr).detCntsPtr)) * ((fltarr(ndet)+1)#factor)
(*(*sumPtr).detCntsPtr) = (*(*sumPtr).detCntsPtr) * ((fltarr(ndet)+1)#factor)
(*(*sumPtr).timePtr) = (*(*sumPtr).timePtr) * factor

; Assign sum to the appropriate destination
*dataPtr = [*dataPtr,sumPtr]  
switch dataType of 
    'samp': begin               ; add the sumPtr to sample dataset
        (*sPtr).sampIndex = nd
        (*sPtr).sampStatus = 1
        break
    end
    'bkgd': begin               ; add the sumPtr to bkgd dataset
        (*sPtr).bkgdIndex = nd
        (*sPtr).bkgdStatus = 1
        break
    end 
    'fast': begin               ; add the sumPtr to fast bkgd dataset
        (*sPtr).fastIndex = nd
        (*sPtr).fastStatus = 1
        break
    end 
endswitch 

end
;-------------------------------------------------------------------------------



;===============================================================================
; fans_ApplySettings
;
; PURPOSE:
;   Whenever settings for treatment of multiple files are modified
;   from the settings dialog, this procedure is called to apply them.
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
;   old_eunit [in] - Flag indicating the existing energy units
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_ApplySettings, sPtr, old_eunit, err=err, emsg=emsg
err=0  & emsg=''

; if energy units have changed then apply it
new_eunit = (*(*sPtr).mfoPtr)[3]
if (new_eunit ne old_eunit) then begin
    ;; if new_eunit==0 then we are converting wavenumber-->meV
    ;; else converting meV-->wavenumber
    conv_factor = (new_eunit eq 0)? 1.0/((*sPtr).mev2wnos) : (*sPtr).mev2wnos 
    
    ;; convert energy units of sample data
    if ((*sPtr).sampCount gt 0) then begin
        for i=0, n_elements((*(*sPtr).sDataPtr))-1 do begin
            (*(*(*(*sPtr).sDataPtr)[i]).energyPtr) = (*(*(*(*sPtr).sDataPtr)[i]).energyPtr)*conv_factor
        endfor
    endif

    ;; convert energy units of bkgd data
    if ((*sPtr).bkgdCount gt 0) then begin
        for i=0, n_elements((*(*sPtr).bDataPtr))-1 do begin
            (*(*(*(*sPtr).bDataPtr)[i]).energyPtr) = (*(*(*(*sPtr).bDataPtr)[i]).energyPtr)*conv_factor
        endfor
    endif

    ;; convert energy units of fast bkgd data
    if ((*sPtr).fastCount gt 0) then begin
        for i=0, n_elements((*(*sPtr).fDataPtr))-1 do begin
            (*(*(*(*sPtr).fDataPtr)[i]).energyPtr) = (*(*(*(*sPtr).fDataPtr)[i]).energyPtr)*conv_factor
        endfor
        ;; force a fast bkgd fit to be performed
        ;(*sPtr).fitPresent = 0
        fans_Fastfit, sPtr=sPtr
    endif

endif

; Apply other required setting changes

; Deal with Sample data
if ((*sPtr).sampCount gt 1) then begin

    switch (*(*sPtr).mfoPtr)[0] of
        
        ;; data to be summed
        0: begin
            if ((*sPtr).sampStatus eq 1) then begin ; summed data exists
                (*sPtr).sampIndex = (*sPtr).sampCount ; - just point to it
                break
            endif 
            fans_sumData, sPtr, 'samp', err=err, emsg=emsg
            break
        end
        
        ;; data to be merged
        1: begin
            if ((*sPtr).sampStatus eq 2) then begin ; merged data exists
                (*sPtr).sampIndex = (*sPtr).sampCount ; - just point to it
                break
            endif 
            fans_apdData, sPtr, 'samp', err=err, emsg=emsg
            break
        end
        
        ;; data to be treated separately
        2: begin
            (*sPtr).sampIndex = (*sPtr).sampSepIndex
            break
        end
        
    endswitch 
    
endif 

; Deal with Bkgd data
if ((*sPtr).bkgdCount gt 1) then begin

    switch (*(*sPtr).mfoPtr)[1] of
        
        ;; data to be summed
        0: begin
            if ((*sPtr).bkgdStatus eq 1) then begin ; summed data exists
                (*sPtr).bkgdIndex = (*sPtr).bkgdCount ; - just point to it
                break
            endif 
            fans_sumData, sPtr, 'bkgd', err=err, emsg=emsg
            break
        end
        
        ;; data to be merged
        1: begin
            if ((*sPtr).bkgdStatus eq 2) then begin ; merged data exists
                (*sPtr).bkgdIndex = (*sPtr).bkgdCount ; - just point to it
                break
            endif 
            fans_apdData, sPtr, 'bkgd', err=err, emsg=emsg
            break
        end
        
        ;; data to be treated separately
        2: begin
            (*sPtr).bkgdIndex = (*sPtr).bkgdSepIndex
            break
        end
        
    endswitch 
    
endif 

; Deal with fast bkgd data
if ((*sPtr).fastCount gt 1) then begin

    switch (*(*sPtr).mfoPtr)[2] of
        
        ;; data to be summed
        0: begin
            if ((*sPtr).fastStatus eq 1) then begin ; merged data exists
                (*sPtr).fastIndex = (*sPtr).fastCount ; - just point to it
                break
            endif 
            fans_sumData, sPtr, 'fast', err=err, emsg=emsg
            break
        end
        
        ;; data to be merged
        1: begin
            if ((*sPtr).fastStatus eq 2) then begin ; merged data exists
                (*sPtr).fastIndex = (*sPtr).fastCount ; - just point to it
                break
            endif 
            fans_apdData, sPtr, 'fast', err=err, emsg=emsg
            break
        end
        
        ;; data to be treated separately
        2: begin
            (*sPtr).fastIndex = (*sPtr).fastSepIndex
            break
        end
        
    endswitch 
    
endif 


;; reset output calculation flag
(*sPtr).doOutpCalc = 1

;; Update plot and/or slider using appropriate data
widget_control, (*sPtr).sViewID, get_uvalue=sViewOn
widget_control, (*sPtr).bViewID, get_uvalue=bViewOn
widget_control, (*sPtr).fViewID, get_uvalue=fViewOn
widget_control, (*sPtr).oViewID, get_uvalue=oViewOn
case 1 of
    sViewOn: begin
        if ((*sPtr).sampCount lt 1) then return
        index = (*sPtr).sampIndex
        dataPtr = (*(*sPtr).sDataPtr)[index]
        updatePlotSlider, sPtr,'samp',dataPtr,err=err,emsg=emsg
    end

    bViewOn: begin
        if ((*sPtr).bkgdCount lt 1)  then return
        index = (*sPtr).bkgdIndex
        dataPtr = (*(*sPtr).bDataPtr)[index]
        updatePlotSlider, sPtr,'bkgd',dataPtr,err=err,emsg=emsg
    end

    fViewOn: begin
        if ((*sPtr).fastCount lt 1)  then return
        fans_plotData, sPtr,'fast',err=err,emsg=emsg ; can't view det cuts
    end 

    oViewOn: begin
        if ((*sPtr).sampCount lt 1)  then return
        fans_outpCalc, sPtr, err=err, emsg=emsg
        if (err ne 0) then begin
            printerror,err,emsg
            return
        endif
        fans_plotData, sPtr,'outp',err=err,emsg=emsg ; can't view det cuts
    end 

endcase
        
end 
;-------------------------------------------------------------------------------



;===============================================================================
; fans_ApplySettings
;
; PURPOSE:
;   Event handler. Deal with changes in ouput options (mask, bkgd sub
;   and fast bkgd sub settings).
;
; PARAMETERS:
;   event [in] - event structure
;
pro fans_outpOpts, event

widget_control, event.top, get_uvalue=sPtr

; store the state of the non-exclusive button in its uvalue
widget_control, event.id, set_uvalue=event.select

; reset status of output calculation flag
(*sPtr).doOutpCalc = 1

;; if the display option is 'Output', then update output plot
widget_control, (*sPtr).oViewID, get_uvalue=on
if (on) then begin
    ;; Calculate output data
    fans_outpCalc, sPtr, err=err, emsg=emsg
    if (err ne 0) then begin
        printerror,err,emsg
        return
    endif
    fans_plotData, sPtr, 'outp', err=err,emsg=emsg
endif

end 
;-------------------------------------------------------------------------------



;===============================================================================
; fans_OutpCalc
;
; PURPOSE:
;   Calculate the final reduced data according to the current dialog
;   settings.
;   Does the following:
;   1) Remove masked detectors if requested
;   2) Remove fast background (using fit parameters)
;   3) Remove sample independent background (empty can)
;   4) For fast bkgd subtraction, it is first normalized to counts
;      per minute and then scaled to the sample (or bkgd) scan.
;   5) When subtracting empty can bkgd, time or monitor normalization
;      can be performed. Again the bkgd is _scaled_to_ the sample.
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_OutpCalc, sPtr, err=err, emsg=emsg
err=0  &  emsg=''


; Only proceed if the doOutpCalc flag is set and 
; sample data exist
if ( (*sPtr).doOutpCalc eq 0 or (*sPtr).sampCount lt 1) then begin
    if ((*sPtr).outpCount eq 0) then begin
        err=1
        emsg='Output Calculation not performed! - No output data available'
    endif
    return
endif

; delete any existing output data
if (ptr_valid((*sPtr).oDataPtr)) then $
  fans_cleanDataPtrs, sPtr, 'outp'

; get ptr to sample data
sampDataPtr = (*(*sPtr).sDataPtr)[(*sPtr).sampIndex] 

; are masked detectors to be removed
widget_control, (*sPtr).mskOnID, get_uvalue=mskon
if (mskon) then begin 
    goodDets = fans_maskDets(sPtr)
    if ( (ptr_valid((*sPtr).mskListPtr) eq 0) or (goodDets[0] eq -1)) then mskon=0
endif

; should a fast bkgd be subtracted
widget_control, (*sPtr).subFastID, get_uvalue=subfast
if (subfast) then begin
    ;; if masked dets have changed since last fast fit then warn user
    ;; and exit
    if (mskon and (*sPtr).updFast) then begin
        message = ['Masked detectors have altered since last fit to fast bkgd!', $
                   'Please redo fast fit before displaying output']
        void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/information)
        err=1
        emsg = 'In "fans_OutpCalc" - fast fit is not in sync with mask dets'
        return 
    endif      
    ;; make sure fit parameters exist
    if (not (*sPtr).fitPresent) then begin 
        message = ['Fast bkgd fit has not been performed!', $
                   'Please do that before proceeding.']
        void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/error)
        err = 2
        emsg = 'In "fans_OutpCalc" - No fast fit parameters present'
        return 
    endif 
endif


; should bkgd be subtracted
nspec1 = (*sampDataPtr).nspec
widget_control, (*sPtr).subBkgdID, get_uvalue=subbkgd
if (subbkgd and (*sPtr).bkgdCount ge 1) then begin
    if (subbkgd) then $
      bkgdDataPtr = (*(*sPtr).bDataPtr)[(*sPtr).bkgdIndex] ; get ptr bkgd data
    nspec2 = (*bkgdDataPtr).nspec
    if (nspec1 ne nspec2) then begin
        message = ['Sample and bkgd files are dissimilar!', $
                   'They have different energy scans (range)', $
                   'Please examine your current file selections.']
        void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/error)
        err = 3
        emsg = 'In "fans_OutpCalc" - Samp & bkgd have different energies!' 
        return
    endif
endif


; now calculate output
;
; 1) make copies of samp and bkgd data
outpPtr = fans_initDataPtr(sampDataPtr)
(*(*outpPtr).historyPtr) = ['DAVE - FANS: Session log. Initialised on '+systime(), $
                            '----------------------------------------------------------', $
                         ' ', $
                         (*(*outpPtr).historyPtr)]
if (subbkgd) then begin 
    bkgdPtr = fans_initDataPtr(bkgdDataPtr)
    (*(*outpPtr).historyPtr) = [(*(*outpPtr).historyPtr), $
                                ' ', $
                                '____________________________________________________', $
                                'Sample independent background subtraction made', $
                                (*(*bkgdPtr).historyPtr)] 
endif 

; 2) apply mask to samp and bkgd by recalculating counts at each energy
if (mskon) then begin
    for i = 0,nspec1-1 do begin
        ;; in exp below, gooddets refer to the dets so need to sub 1
        ;; to get indices
        (*(*outpPtr).totCntsPtr)[i] = total((*(*outpPtr).detCntsPtr)[goodDets-1,i])
        (*(*outpPtr).totErrsPtr)[i] = sqrt(total(((*(*outpPtr).detErrsPtr)[goodDets-1,i])^2))
        if (subbkgd) then begin
            (*(*bkgdPtr).totCntsPtr)[i] = total((*(*bkgdPtr).detCntsPtr)[goodDets-1,i])
            (*(*outpPtr).totErrsPtr)[i] = sqrt(total(((*(*bkgdPtr).detErrsPtr)[goodDets-1,i])^2))
        endif
    endfor
endif

; 3) subtract fast bkgd from samp and bkgd - fast bkgd is scaled up before subtraction
if (subfast) then begin
    omega = (*(*outpPtr).energyPtr)
    ;fpar[5] = ((*(*sPtr).mfoPtr)[3] eq 0)? 1.0 : (*sPtr).mev2wnos ; convert factor b/n meV <--> wavenumbers 
    ;fastfit = (fans_fitfunc(omega,fpar))[*,0] ;; b/c fuction returns values + partial derivatives wrt fpar
    fpar = (*(*sPtr).fitParsPtr)
    parstr = strtrim(string(fpar,format='(G10.4)'),2)
    case (*sPtr).fitFunc of
        0: begin
            funcName = 'Arcsin'
            fastfit = (fans_Arcsin(omega,fpar))[*,0]
        end
        1: begin
            funcName = 'Linear'
            fastfit = (fans_Linear(omega,fpar))[*,0]
        end
        2: begin
            funcName = 'Quadratic'
            fastfit = (fans_Quad(omega,fpar))[*,0]
        end
        3: begin
            funcName = 'Cubic'
            fastfit = (fans_Cubic(omega,fpar))[*,0]
        end
    endcase
    if ((*sPtr).fastCount ge 1) then begin
        ;; get fast bkgd counting time - usually fast bkgd is
        ;; collected in fixed time mode so just use the fist value
        fastPtr = (*(*sPtr).fDataPtr)[(*sPtr).fastIndex]
        ftime = (*(*fastPtr).timePtr)[0]
        ;; update historyPtr
        (*(*outpPtr).historyPtr) = [(*(*outpPtr).historyPtr), $
                                    ' ', $
                                    '____________________________________________________', $
                                    (*(*fastPtr).historyPtr), $
                                    'Fast bkgd correction performed ... ', $
                                    'Model function: '+funcName, $
                                    '--- fit parameters used: '+parstr] 
    endif else begin
        ftime = 1.0             ; no fast data so assume fit is in per 1 minute counts
        ;fpar[4] = 3.35
        (*(*outpPtr).historyPtr) = [(*(*outpPtr).historyPtr), $
                                    ' ', $
                                    '____________________________________________________', $
                                    'Fast bkgd correction performed ... ', $
                                    'Model function: '+funcName, $
                                    '--- fit parameters used: '+parstr] 
    endelse 
    (*(*outpPtr).totCntsPtr) = (*(*outpPtr).totCntsPtr) - fastfit/ftime * (*(*outpPtr).timePtr) 
    if (subbkgd) then $
      (*(*bkgdPtr).totCntsPtr) = (*(*bkgdPtr).totCntsPtr) - fastfit/ftime * (*(*bkgdPtr).timePtr)
endif

; 4) Subtract bkgd if necessary - bkgd is (monitor) scaled to samp before subtraction
if (subbkgd) then begin
    smon = (*outpPtr).monitor * (*outpPtr).prefactor
    bmon = (*bkgdPtr).monitor * (*bkgdPtr).prefactor
    scalefactor = float(smon) / float(bmon)
    bkgdCnts = (*(*bkgdPtr).totCntsPtr) * scalefactor
    bkgdCErr = (*(*bkgdPtr).totErrsPtr) * scalefactor
    (*(*outpPtr).totCntsPtr) = (*(*outpPtr).totCntsPtr) - bkgdCnts
    (*(*outpPtr).totErrsPtr) = sqrt((*(*outpPtr).totErrsPtr)^2 + bkgdCerr^2)
endif

; 5) Move results to sPtr
tmparr = ptrarr(1)
tmparr[0] = outpPtr
(*sPtr).outpCount = 1
(*sPtr).outpIndex = 0
(*sPtr).oDataPtr = ptr_new(tmparr,/no_copy)
(*sPtr).doOutpCalc = 0

; remove any temperal heap data that was created
cleanDataPtr, bkgdPtr

;; update historyPtr with detector info
mTxt = strtrim((*outpPtr).monitor * (*outpPtr).prefactor,2)
text = ['Monitor counts = '+mTxt, $
        'Number of active detectors = '+strtrim((*outpPtr).ndets,2), $
        '--- the first two are monitors']
if (mskon) then begin
    widget_control, (*sPtr).mskListID, get_value=lstring
    text = [text, $
            'Number of masked detectors = '+strtrim(n_elements(*(*sPtr).mskListPtr),2), $
            'Masked detector list = '+lstring]
endif 
(*(*outpPtr).historyPtr) = [(*(*outpPtr).historyPtr), $
                         ' ', $
                         text]

; no errors
err=0  &  emsg=''

end
;-------------------------------------------------------------------------------



;===============================================================================
; fans_2Dave
;
; PURPOSE:
;   Save reduced data (integrated counts vs energy transfer) to a DAVE format file
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_2Dave, sPtr, err=err, emsg=emsg
err=0  &  emsg=''

; nothing to do if no sample data present
if ((*sPtr).sampCount lt 1) then begin
    err=1
    emsg='No sample data available - cannot determine output!'
    return
endif 

; Is output data up to date? - if not, evaluate it
if ( (*sPtr).doOutpCalc) then begin
    fans_OutpCalc, sPtr, err=err, emsg=emsg
    if (err ne 0) then return
endif

; copy data into the davePtr
; NOTE: the davePtr is passed down from DAVE.pro! This is not a good
; thing b/c there is no way of knowing what type of data it contains
; sice each instrument has a different specificPtr. Problems arise
; when the specificPtr contains ptrs, it will not be possible to
; clean the ptrs before rewriting the specificPtr. Thus there are
; memory leaks! 
;
; ########### Need a better solution!!!!!!!!!!!
;
davePtr = (*sPtr).davePtr
outpPtr = (*(*sPtr).oDataPtr)[(*sPtr).outpIndex]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 1) deal with the davePtr.commonStr

; first, free-up any pointers within commonStr
free_Ptr,(*(*davePtr).dataStrPtr).commonStr.histPtr
free_Ptr,(*(*davePtr).dataStrPtr).commonStr.treatmentPtr


qty = (*(*outpPtr).totCntsPtr)
qerr = (*(*outpPtr).totErrsPtr)
;Data is already normalized to monitor of first dataset 
;If using a custom monitor to normalize then need to rescale
monConstant = (*outpPtr).monitor * (*outpPtr).prefactor
wID = widget_info((*sPtr).tlbID,find_by_uname='USEMON')
if (widget_info(wID,/valid)) then begin
   widget_control, wID, get_uvalue=isSet
   if (isSet) then begin
      wID = widget_info((*sPtr).tlbID,find_by_uname='MONVAL')
      if (widget_info(wID,/valid)) then begin
         widget_control, wID, get_value=monValue
         if (monValue gt 0.0) then begin
             normFac = float(monValue)/monConstant
             qerr = qerr * normFac
             qty = qty * normFac
             monConstant=monValue
         endif
      endif
   endif
endif

;; title
switch (*outpPtr).base of
    'TIME': begin
        ftime = fix((*(*outpPtr).timePtr)[0] * 60)
        histlabel = 'Counts per '+strtrim(string(ftime),2)+' secs'
        break
    end
    
    'NEUT':
    else: begin
        monCnts = (*outpPtr).monitor * (*outpPtr).prefactor
        histlabel = 'Counts per '+strtrim(string(monConstant),2)+' mon'
        break
    end
endswitch
        
; fill it up
;dummy=strarr(1)

(*(*davePtr).dataStrPtr).commonStr.instrument = 'FANS'
(*(*davePtr).dataStrPtr).commonStr.histPtr = $
  ptr_new({qty:qty, $ ; *** may need to redefine as 2D
           err:qerr, $ ; *** may need to redefine as 2D
           x:(*(*outpPtr).energyPtr), $
           y:0.0 $              ; *** not clear what to do here - FANS data is 1D!
          })
(*(*davePtr).dataStrPtr).commonStr.xlabel = 'Energy Transfer'
(*(*davePtr).dataStrPtr).commonStr.xunits = ((*(*sPtr).mfoPtr)[3] eq 0)? 'meV' : '!ncm!u-1!n'
(*(*davePtr).dataStrPtr).commonStr.xtype = 'points'
(*(*davePtr).dataStrPtr).commonStr.ylabel = ''
(*(*davePtr).dataStrPtr).commonStr.yunits = ''
(*(*davePtr).dataStrPtr).commonStr.ytype = 'points'
monCnts = (*outpPtr).monitor * (*outpPtr).prefactor

(*(*davePtr).dataStrPtr).commonStr.histlabel = histlabel
(*(*davePtr).dataStrPtr).commonStr.histunits = ''
(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = ptr_new((*(*outpPtr).historyPtr))

;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 2) deal with davePtr.specificPtr - instrument specific attributes

; brutally clean specificPtr - have to use heap_free in order to
;                              recursively deal with inner ptrs
heap_free,(*(*davePtr).dataStrPtr).specificPtr

specificStr = { $
            filename:(*outpPtr).filename, $
            date:(*outpPtr).date, $
            scan:(*outpPtr).scan, $
            monitor:(*outpPtr).monitor, $
            prefactor:(*outpPtr).prefactor, $
            base:(*outpPtr).base, $
            type:(*outpPtr).type, $
            title:(*outpPtr).title, $
            collimation:(*outpPtr).collimation, $
            mosaic:(*outpPtr).mosaic, $
            orientation:(*outpPtr).orientation, $
            latticeA:(*outpPtr).latticeA, $
            latticeB:(*outpPtr).latticeB, $
            latticeC:(*outpPtr).latticeC, $
            latticeAl:(*outpPtr).latticeAl, $
            latticeBe:(*outpPtr).latticeBe, $
            latticeGa:(*outpPtr).latticeGa, $
            eStart:(*outpPtr).eStart, $
            eFinal:(*outpPtr).eFinal, $
            eDelta:(*outpPtr).eDelta, $
            monDSpace:(*outpPtr).monDSpace, $
            anaDSpace:(*outpPtr).anaDSpace, $
            tStart:(*outpPtr).tStart, $
            tInc:(*outpPtr).tInc, $
            qhCenter:(*outpPtr).qhCenter, $
            qkCenter:(*outpPtr).qkCenter, $
            qlCenter:(*outpPtr).qlCenter, $
            qhDelta:(*outpPtr).qhDelta, $
            qkDelta:(*outpPtr).qkDelta, $
            qlDelta:(*outpPtr).qlDelta, $
            hfield:(*outpPtr).hfield, $
            tflag:(*outpPtr).tflag, $
            ndets:(*outpPtr).ndets, $
            nspec:(*outpPtr).nspec $
           }

(*(*davePtr).dataStrPtr).specificPtr = ptr_new(specificStr)


; now save to disk
index = (*sPtr).sampIndex
if ( index eq (*sPtr).sampCount) then $
  index = 0

def_ofile = (*(*sPtr).outpListPtr)[index]+'.dave'

repeat begin
    okay = 1
    ofile = dialog_pickfile(dialog_parent=(*sPtr).tlbID, $
                            title="Select output filename", $
                            path=(*sPtr).workDir, $
                            file=def_ofile, $
                            get_path=newPath, $
                            filter='*.dave' $
                           )

    if (strtrim(ofile,2) eq '') then return ; dialog_pickfile was cancelled so do nothing further

    if ( strmid(ofile,strlen(ofile)-1,1) eq !delimiter) then $ ; ie only a directory was return...
      ofile = ofile+def_ofile   ; then append def_ofile to it
    
    ;; if file exists then make sure user wishes to overide it
    void = findfile(ofile,count=cnt)
    if (cnt) then begin
        message = [ofile+' exists','Do you want to overwrite it?']
        ans = dialog_message(message,dialog_parent=(*sPtr).tlbid,/question)
        okay = (strtrim(ans,2) eq 'Yes')? 1 : 0 
    endif
endrep until okay

if (strlen(newPath) gt 1) then $
  (*sPtr).workDir = newPath     ; update current working directory

; write to disk and inform user of ouput destination
save, davePtr, filename=ofile
message=['Output saved in:',ofile]
void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/information)

end
;-------------------------------------------------------------------------------



;===============================================================================
; fans_2DDave
;
; PURPOSE:
;   Save reduced data (individual detector counts) to a DAVE
;   format file. The save file is a 2D file
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_2DDave, sPtr, err=err, emsg=emsg
err=0  &  emsg=''

; nothing to do if no sample data present
if ((*sPtr).sampCount lt 1) then begin
    err=1
    emsg='No sample data available - cannot determine output!'
    return
endif 

; Is output data up to date? - if not, evaluate it
if ( (*sPtr).doOutpCalc) then begin
    fans_OutpCalc, sPtr, err=err, emsg=emsg
    if (err ne 0) then return
endif

; copy data into the davePtr
; NOTE: the davePtr is passed down from DAVE.pro! This is not a good
; thing b/c there is no way of knowing what type of data it contains
; sice each instrument has a different specificPtr. Problems arise
; when the specificPtr contains ptrs, it will not be possible to
; clean the ptrs before rewriting the specificPtr. Thus there are
; memory leaks! 
;
; ########### Need a better solution!!!!!!!!!!!
;
davePtr = (*sPtr).davePtr
outpPtr = (*(*sPtr).oDataPtr)[(*sPtr).outpIndex]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 1) deal with the davePtr.commonStr

; first, free-up any pointers within commonStr
free_Ptr,(*(*davePtr).dataStrPtr).commonStr.histPtr
free_Ptr,(*(*davePtr).dataStrPtr).commonStr.treatmentPtr

;; title
switch (*outpPtr).base of
    'TIME': begin
        ftime = fix((*(*outpPtr).timePtr)[0] * 60)
        histlabel = 'Detector Cnts per '+strtrim(string(ftime),2)+' secs'
        break
    end
    
    'NEUT':
    else: begin
        monCnts = (*outpPtr).monitor * (*outpPtr).prefactor
        histlabel = 'Detector Cnts per '+strtrim(string(monCnts),2)+' mon'
        break
    end
endswitch     

; fill it up
ndets = (*outpPtr).ndets
(*(*davePtr).dataStrPtr).commonStr.instrument = 'FANS'
(*(*davePtr).dataStrPtr).commonStr.histPtr = $
  ptr_new({qty:transpose((*(*outpPtr).detCntsPtr)[2:ndets-1,*]), $ 
           err:transpose((*(*outpPtr).detErrsPtr)[2:ndets-1,*]), $ 
           x:(*(*outpPtr).energyPtr), $
           y:indgen(ndets-2)+3 $ 
          })
(*(*davePtr).dataStrPtr).commonStr.xlabel = 'Energy Transfer'
(*(*davePtr).dataStrPtr).commonStr.xunits = ((*(*sPtr).mfoPtr)[3] eq 0)? 'meV' : '!ncm!u-1!n'
(*(*davePtr).dataStrPtr).commonStr.xtype = 'points'
(*(*davePtr).dataStrPtr).commonStr.ylabel = 'Detector Nos'
(*(*davePtr).dataStrPtr).commonStr.yunits = ''
(*(*davePtr).dataStrPtr).commonStr.ytype = 'points'
monCnts = (*outpPtr).monitor * (*outpPtr).prefactor
(*(*davePtr).dataStrPtr).commonStr.histlabel = histlabel
(*(*davePtr).dataStrPtr).commonStr.histunits = ''
(*(*davePtr).dataStrPtr).commonStr.treatmentPtr = ptr_new((*(*outpPtr).historyPtr))
(*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr) = $
  [(*(*(*davePtr).dataStrPtr).commonStr.treatmentPtr), $
   '=================================================', $
   'The data in the 2D FANS files have not been ', $
   'corrected for instrument fast background. ', $
   'The data are essentially the raw detector ', $
   'counts for the sample file(s) listed above. ', $
   'However, the following are performed: ', $
   '    * if present, subtraction of container scatt.', $
   '    * proper monitor norm. if more than one file.' ]  

;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 2) deal with davePtr.specificPtr - instrument specific attributes

; brutally clean specificPtr - have to use heap_free in order to
;                              recursively deal with inner ptrs
heap_free,(*(*davePtr).dataStrPtr).specificPtr

specificStr = { $
            filename:(*outpPtr).filename, $
            date:(*outpPtr).date, $
            scan:(*outpPtr).scan, $
            monitor:(*outpPtr).monitor, $
            prefactor:(*outpPtr).prefactor, $
            base:(*outpPtr).base, $
            type:(*outpPtr).type, $
            title:(*outpPtr).title, $
            collimation:(*outpPtr).collimation, $
            mosaic:(*outpPtr).mosaic, $
            orientation:(*outpPtr).orientation, $
            latticeA:(*outpPtr).latticeA, $
            latticeB:(*outpPtr).latticeB, $
            latticeC:(*outpPtr).latticeC, $
            latticeAl:(*outpPtr).latticeAl, $
            latticeBe:(*outpPtr).latticeBe, $
            latticeGa:(*outpPtr).latticeGa, $
            eStart:(*outpPtr).eStart, $
            eFinal:(*outpPtr).eFinal, $
            eDelta:(*outpPtr).eDelta, $
            monDSpace:(*outpPtr).monDSpace, $
            anaDSpace:(*outpPtr).anaDSpace, $
            tStart:(*outpPtr).tStart, $
            tInc:(*outpPtr).tInc, $
            qhCenter:(*outpPtr).qhCenter, $
            qkCenter:(*outpPtr).qkCenter, $
            qlCenter:(*outpPtr).qlCenter, $
            qhDelta:(*outpPtr).qhDelta, $
            qkDelta:(*outpPtr).qkDelta, $
            qlDelta:(*outpPtr).qlDelta, $
            hfield:(*outpPtr).hfield, $
            tflag:(*outpPtr).tflag, $
            ndets:(*outpPtr).ndets, $
            nspec:(*outpPtr).nspec $
           }

(*(*davePtr).dataStrPtr).specificPtr = ptr_new(specificStr)


; now save to disk
index = (*sPtr).sampIndex
if ( index eq (*sPtr).sampCount) then $
  index = 0

def_ofile = (*(*sPtr).outpListPtr)[index]+'_2D.dave'

repeat begin
    okay = 1
    ofile = dialog_pickfile(dialog_parent=(*sPtr).tlbID, $
                            title="Select output filename", $
                            path=(*sPtr).workDir, $
                            file=def_ofile, $
                            get_path=newPath, $
                            filter='*.dave' $
                           )

    if (strtrim(ofile,2) eq '') then return ; dialog_pickfile was cancelled so do nothing further

    if ( strmid(ofile,strlen(ofile)-1,1) eq !delimiter) then $ ; ie only a directory was return...
      ofile = ofile+def_ofile   ; then append def_ofile to it

    ;; if file exists then make sure user wishes to overide it
    void = findfile(ofile,count=cnt)
    if (cnt) then begin
        message = [ofile+' exists','Do you want to overwrite it?']
        ans = dialog_message(message,dialog_parent=(*sPtr).tlbid,/question)
        okay = (strtrim(ans,2) eq 'Yes')? 1 : 0 
    endif
endrep until okay

if (strlen(newPath) gt 1) then $
  (*sPtr).workDir = newPath     ; update current working directory

; write to disk and inform user of ouput destination
save, davePtr, filename=ofile
message=['Output saved in:',ofile]
void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/information)

end
;-------------------------------------------------------------------------------



;===============================================================================
; fans_2Ascii
;
; PURPOSE:
;   Save reduced data to an ASCII file
;
; PARAMETERS:
;   sPtr [in] - pointer to state structure
;
; KEYWORDS:
;   err [out] - err number, if an error occured
;
;   emsg [out] - explanatory error message, if an error occured
;
pro fans_2Ascii, sPtr, err=err, emsg=emsg
err=0  &  emsg=''

; nothing to do if no sample data present
if ((*sPtr).sampCount lt 1) then begin
    err=1
    emsg='No sample data available - cannot determine output!'
    return
endif 

; Is output data up to date? - if not, evaluate it
if ( (*sPtr).doOutpCalc) then begin
    fans_OutpCalc, sPtr, err=err, emsg=emsg
    if (err ne 0) then return
endif

; get the data for output
outpPtr = (*(*sPtr).oDataPtr)[(*sPtr).outpIndex]

; now save to disk
index = (*sPtr).sampIndex
if ( index eq (*sPtr).sampCount) then $
  index = 0

def_ofile = (*(*sPtr).outpListPtr)[index]+'.dat'

repeat begin
    okay = 1
    ofile = dialog_pickfile(dialog_parent=(*sPtr).tlbID, $
                            title="Select output filename", $
                            path=(*sPtr).workDir, $
                            file=def_ofile, $
                            get_path=newPath, $
                            filter='*.dat' $
                           )
    if (strtrim(ofile,2) eq '') then return ; dialog_pickfile was cancelled so do nothing further
    
    if ( strmid(ofile,strlen(ofile)-1,1) eq !delimiter) then $ ; ie only a directory was returned...
      ofile = ofile+def_ofile   ; then append def_ofile to it

    ;; if file exists then make sure user wishes to overide it
    void = findfile(ofile,count=cnt)
    if (cnt) then begin
        message = [ofile+' exists','Do you want to overwrite it?']
        ans = dialog_message(message,dialog_parent=(*sPtr).tlbid,/question)
        okay = (strtrim(ans,2) eq 'Yes')? 1 : 0 
    endif
endrep until okay

if (strlen(newPath) gt 1) then $
  (*sPtr).workDir = newPath     ; update current working directory
;(*sPtr).workDir = newPath

; open input file
openw, lun, ofile, /get_lun, error=err
if (err ne 0) then begin        ; if error, get the error msg and exit
    emsg = 'In "fans_readData" - '+!error_state.msg    
    return
endif

; write to disk and inform user of ouput destination
omega = (*(*outpPtr).energyPtr)
cnts = (*(*outpPtr).totCntsPtr)
cnts_err = (*(*outpPtr).totErrsPtr)
time = (*(*outpPtr).timePtr)
temp = (*(*outpPtr).tempPtr)

;Data is already normalized to monitor of first dataset 
;If using a custom monitor to normalize then need to rescale
monConstant = (*outpPtr).monitor * (*outpPtr).prefactor
wID = widget_info((*sPtr).tlbID,find_by_uname='USEMON')
if (widget_info(wID,/valid)) then begin
   widget_control, wID, get_uvalue=isSet
   if (isSet) then begin
      wID = widget_info((*sPtr).tlbID,find_by_uname='MONVAL')
      if (widget_info(wID,/valid)) then begin
         widget_control, wID, get_value=monValue
         if (monValue gt 0.0) then begin
             normFac = float(monValue)/monConstant
             cnts_Err = cnts_Err * normFac
             Cnts = Cnts * normFac
             monConstant=monValue
         endif
      endif
   endif
endif

nd = n_elements(omega)
for i = 0,nd-1 do printf, lun, format='(5(F11.4,2X))',omega[i],cnts[i],cnts_err[i],time[i],temp[i]

free_lun, lun, /force

message=['Output saved in:',ofile]
void = dialog_message(message,dialog_parent=(*sPtr).tlbID,/information)

end
