/*      FILE: /mxtools/src/strawman/sshift_ansi.c
 *   PURPOSE: 
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001OCT26
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 *
 * COMMENT: 
 *
 *  This ANSI C function version of the K&R C function "sshift.c" from 
 *  the ZODIAC C library by 
 *
 *    Marc W. Buie (buie@lowell.edu)
 *    Lowell Observatory
 *    1400 West Mars Hill Road
 *    Flagstaff, AZ  86001
 *
 * CHANGES:
 *
 *   2001OCT26  Renamed pi macro to PI.
 *              Undefined all local macros at the end of the function.
 *
 *   2001OCT25  Converted function from K&R C to ANSI C.
 *
 */

#include <stdio.h>
#include <math.h>

void sshift_ansi(
  float x[],   /* Input data array. */
  int	n,     /* Number of points to shift. */
  float shift, /* Amount of the shift. */
  float hole,  /* value to assign to out of range data. */
  float	xp[]   /* Destination array for the shifted data. */
){

/*   This procedure 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].
 *
 *   Note: n is taken to be an absolute limit of the sizes of x and xp so
 * data will fall off one end or another 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.
 *
 */

#define DAMPFAC 3.25
#define	EPS	1.0e-5
#define PI	3.1415926535897932

   float  sinc[21]; /* This is the array containing the convolution points. */
   int    point;    /* This is a counter variable. */
   float  y;        /* Another counter variable. */
   int    ishift;   /* Integer part of shift. */
   float  fshift;   /* Fractional part of shift. */
   int    lobe,npix,i;     /* Counters */

/* First split the desired shift into a fractional and integer part. */
   ishift = shift;    
   fshift = shift - ishift;

/* Do the fractional shift first (if necessary). */
   if ( ( fabs(fshift) > EPS ) && ( fabs(fshift) < 1.0-EPS) ) {

/* Initialize the sinc array. */
      for ( y=fshift+10, point=0 ; point < 21; --y, point++ )
         sinc[point] = exp(-y*y/DAMPFAC/DAMPFAC)*sin(PI*y)/(PI*y);

/* Convolve the sinc array with the input data.  This is the shift. */
      for ( point=0; point < n; point++) 
         for( lobe=0, xp[point]=0; lobe < 21; lobe++) {
            npix = point - (lobe-10);
            if( (npix >= 0) && (npix < n) ) 
               xp[point] += sinc[lobe]*x[npix];
            else 
               xp[point] += sinc[lobe]*hole;
            }
      }
   else for( i=0; i<n; i++) xp[i]=x[i];

/* Now perform the integer shift. */
   if (ishift < 0) /* shift is negative */
      for ( i=0,point = -ishift; i < n; i++,point++) {
         if( point < 0 ) xp[i]=hole;
         else if( point >= n ) xp[i]=hole;
         else xp[i]=xp[point];
         }
   else /* shift is positive */
      for ( i=n-1,point = n-1-ishift; i >= 0; --i,--point) {
         if( point < 0 ) xp[i]=hole;
         else if( point >= n ) xp[i]=hole;
         else xp[i]=xp[point];
         }

   return;

}
#undef DAMPFAC
#undef EPS	
#undef PI	
/* end-of-file */