;+
; Shift the data in the data container by the given number of channels 
; (which may be a floating point number containing fractional
; channels) such that data in channel i before the shift is found in 
; channel i+offset after the shift.  The data's reference channel is
; also shifted so that features should remain stationary when
; displayed except when the x-axis is channels.  Shifted data is
; replaced with zeros unless the /wrap keyword is set in which case
; the data wraps around as necessary during the shift. Data are first
; shifted by the nearest integer number of channels and then the
; result is shifted by the remaining fractional channel if that
; fraction channel is more than ftol.  
;
; <p>Currently, fractional shifting is
; done using an FFT windowed using a Welch function to reduce ringing
; as a result of the FFT.  Set ftol to a number >= 1 to turn off all
; fractional shifting.  Eventually other interpolation methods will be
; available.
;
; <p> Use one of xshift, vshift, or fshift to calculate the offset to
; align this data with data in an ongoing accumulation.
;
; <p> Windowing using the Welsh function reduces ringing at the
; expense of the high-order terms in the FFT.  This very slightly
; broadens features.  This can be turned off using the /nowelsh
; keyword.
;
; <p> The data are padded with 0s to the next higher power of two as
; an intermediate step (the final result will have the original size
; of the data).  This is done to reduce ringing due to any
; discontinuities between the two ends of the data.  This can be
; turned off using the /nopad keyword.
;
; @param dc {in}{out}{required}{type=spectrum} The data container to
; shift.  The shift is done in place.
;
; @param offset {in}{required}{type=floating point} The number of
; channels to shift the data (positive shifts things towards higher
; channels, negative shifts things towards lower channels).
;
; @keyword wrap {in}{optional}{type=boolean} Data shifted off one end
; of the array appears on the other end of the array (it wraps around
; as a result of the shift) when this is set.  Otherwise, as data is
; shifted it is replaced by 0s and data shifted off the end is lost.
; 
; @keyword ftol {in}{optional}{type=floating point}{default=0.01}
; Fractional shifts (the non-integer portion of offset) are only done
; when they are larger than ftol.  Set this value to >= 1.0 to turn
; off all fractional shifts.
; 
; @keyword nowelsh {in}{optional}{type=boolean} When set, the shifted
; data is NOT windowed using the Welsh function.
;
; @keyword nopad {in}{optional}{type=boolean} When set, the data is
; NOT padded with 0s to the next higher power of 2 prior to the FFT
; and shift.
;
; @uses <a href="data_valid.html">data_valid</a>
;
; @version $Id: dcshift.pro,v 1.2 2005/05/18 19:10:49 bgarwood Exp $
;-
pro dcshift, dc, offset, wrap=wrap, ftol=ftol, nowelsh=nowelsh, nopad=nopad
    compile_opt idl2

    if n_params() ne 2 then begin
        message,'Usage: dcshift, dc, offset [,wrap=wrap, ftol=ftol]',/info
        return
    endif

    nch = data_valid(dc)
    
    if nch le 0 then begin
        message,'No data in data container',/info
        return
    endif

    if abs(offset) ge nch then begin
        message,'offset is more than the number of channels, can not shift',/info
        return
    endif

    if not keyword_set(nowelsh) then nowelsh = 0
    if not keyword_set(nopad) then nopad = 0

    ishift = round(offset)
    fshift = offset-ishift

    if abs(ishift) gt 0 then begin
        ; use builtin shift - which wraps
        *dc.data_ptr = shift(*dc.data_ptr,ishift)

        if not keyword_set(wrap) then begin
            ; zero out the wrapped stuff
            if ishift gt 0 then begin
                ibeg = 0
                iend = ishift-1
            endif else begin
                iend = nch -1
                ibeg = iend + ishift - 1
            endelse
            (*dc.data_ptr)[ibeg:iend] = 0.0
        endif
    endif

    ; do fractional shift here
    if n_elements(ftol) eq 0 then ftol = 0.01

    if (abs(fshift) gt ftol) then begin
        if not nopad then begin
            ; expand the data to next power of 2 and pad with zeros to 
            ; prevent aliasing
            pow2 = round(alog10(nch)/alog10(2))
            pow2 += 1
            newsize = 2^pow2
        endif else begin
            newsize = nch
        endelse
        apad = dblarr(newsize)
        npad = newsize - nch
        nskip = round(npad/2.0)
        apad[nskip:(nskip+nch-1)] = *dc.data_ptr

        ; fft the data
        fta = fft(apad,/inverse,/overwrite)

        ; add in the phase shift and do welch windowing
        phshift = 2.0d * !dpi * fshift/double(newsize)

        amp = sqrt(real_part(fta)^2 + imaginary(fta)^2)
        phs = atan(imaginary(fta),real_part(fta))

        ; do each half separately
        half = (newsize/2)
        farr = findgen(half)
        ; this multiplies by the Welch function
        if not nowelsh then amp[0:(half-1)] *= (1.0 - (farr/float(newsize/2))^2)
        ; and this does the shift
        phs[0:(half-1)] += phshift*farr

        farr += float(half - newsize)
;        ; this multiplies by the Welch function
        if not nowelsh then amp[half:(newsize-1)] *= (1.0 - (farr/float(newsize/2))^2)
        ; and this does this shift
        phs[half:(newsize-1)] += phshift*farr

        ; return to complex form
        fta = complex(amp*cos(phs),amp*sin(phs))

        ; and back to frequency domain
        apad = fft(fta,/overwrite)

        ; and back to a
        *dc.data_ptr = real_part(apad[nskip:(nskip+nch-1)])
    endif else begin
        fshift = 0.0
    endelse
    
    ; and reset the reference pixel
    trueShift = double(ishift) + fshift
    dc.reference_channel += trueShift

end