;+
; NAME:
; sshift
; PURPOSE: (one line)
; Shift data using a damped sinc function for fractional part.
; DESCRIPTION:
;
; This function will shift an array of data pointed to by x and
; extending for n points. The amount of the shift is given by shift.
; The result of the operation is placed at xp. A shift that is within
; 0.0001 of a whole number is treated to be that of the whole number. If
; the shift is by an integral number of pixels then the shift involves
; reindexing the data, no interpolation is done. If the shift is some
; non-integral amount then the data is resampled using a damped sinc
; function.
;
; The sense of the shift is as follows: think of the array plotted on a
; fixed scale. A shift of 1 corresponds to shifting the data by one data
; point to the right relative to the fixed scale, ie. x[3]=xp[4].
;
; The data will fall off one end or another of the output vector as a
; result of the shift. However, this is not as significant as the edge
; effect, the convolution is not complete for any data point within 10
; points of the edge, so those points cannot be trusted. The missing
; points in the convolution are assumed to be equal to the end points.
;
; CATEGORY:
; Numerical
; CALLING SEQUENCE:
; xp = sshift(x,shift)
; INPUTS:
; x - Input data array to be shifted.
; shift - Amount to shift (negative to left).
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
; OUTPUTS:
; Return value is the shifted array.
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; The input and output arrays cannot be the same.
; PROCEDURE:
; MODIFICATION HISTORY:
; Adapted from Zodiac routine: shiftc/sshift
; Marc W. Buie, Lowell Observatory, 1992 October 2
;-
function sshift,x,shift
EPS = 1.0e-5 ; Smallest fractional shift allowed.
DAMPFAC = 3.25 ; Damping factor for gaussian.
NS = 21 ; Number of points in the sinc convolution kernal.
npts = n_elements(x)
xp = fltarr(npts)
; First, split the desired shift into a fractional and integer part.
ishift = fix(shift)
fshift = shift - ishift
; Do the fractional shift first (if necessary).
if ( ( abs(fshift) gt EPS ) and ( abs(fshift) lt 1.0-EPS) ) then begin
; Initialize the sinc array.
y = fshift - (findgen(NS)-10)
py = !pi * y
sinc = exp( -y^2/DAMPFAC^2 ) * sin(py)/py
sinc = rotate(sinc,2)
; Convolve the sinc array with the input data. This is the shift.
for point=0L,n_elements(x)-1 do begin
lobe = (indgen(NS) - 10) + point
vals = fltarr(NS)
z = where(lobe lt 0,count)
if (count ne 0) then vals[z] = x[0]
z = where(lobe ge 0 and lobe lt npts,count)
if (count ne 0) then vals[z] = x[lobe[z]]
z = where(lobe ge npts,count)
if (count ne 0) then vals[z] = x[npts-1]
xp[point] = total(sinc*vals)
endfor
; No fractional shift, just copy the data.
endif else begin
xp = x
endelse
; Now perform the integer shift.
left = -ishift
right = npts-1 - ishift
if (left lt 0) then begin
xp[ishift:npts-1] = xp[0:npts-1-ishift]
xp[0:ishift-1] = xp[ishift]
endif else if (left gt 0) then begin
xp[0:npts-1-left] = xp[left:npts-1]
xp[npts-left:npts-1] = xp[npts-1-left]
endif
return,xp
end