FUNCTION RINTER, P, X, Y, DFDX, DFDY, INITIALIZE = initialize
;+
; NAME:
;      RINTER
; PURPOSE:
;      Cubic interpolation of an image at a set of reference points.
; EXPLANATION:
;      This interpolation program is equivalent to using the intrinsic 
;      INTERPOLATE() function with CUBIC = -0.5.   However,
;      RINTER() has two advantages: (1) one can optionally obtain the 
;      X and Y derivatives at the reference points, and (2) if repeated
;      interpolation is to be applied to an array, then some values can
;      be pre-computed and stored in Common.   RINTER() was originally  
;      for use with the DAOPHOT procedures, but can also be used for 
;      general cubic interpolation.
;
; CALLING SEQUENCE:
;      Z = RINTER( P, X, Y, [ DFDX, DFDY ] )
;               or
;      Z = RINTER(P, /INIT)
;
; INPUTS:                 
;      P  - Two dimensional data array, 
;      X  - Either an N element vector or an N x M element array,
;               containing X subscripts where cubic interpolation is desired.
;      Y -  Either an N element vector or an N x M element array, 
;               containing Y subscripts where cubic interpolation is desired.
;
; OUTPUT:
;      Z -  Result = interpolated vector or array.  If X and Y are vectors,
;              then so is Z, but if X and Y are arrays then Z will be also.
;              If P is DOUBLE precision, then so is Z, otherwise Z is REAL.
;
; OPTIONAL OUTPUT:
;      DFDX - Vector or Array, (same size and type as Z), containing the 
;               derivatives with respect to X
;      DFDY - Array containing derivatives with respect to Y
;
; OPTIONAL KEYWORD INPUT:
;     /INIT - Perform computations associated only with the input array (i.e. 
;             not with X and Y) and store in common.    This can save time if
;             repeated calls to RINTER are made using the same array.  
;        
; EXAMPLE:
;      suppose P is a 256 x 256 element array and X = FINDGEN(50)/2. + 100.
;      and Y = X.  Then Z will be a 50 element array, containing the
;      cubic interpolated points.
;
; SIDE EFFECTS:
;      can be time consuming.
;
; RESTRICTION:
;      Interpolation is not possible at positions outside the range of 
;       the array (including all negative subscripts), or within 2 pixel
;       units of the edge.  No error message is given but values of the 
;       output array are meaningless at these positions.
;
; PROCEDURE:
;       invokes CUBIC interpolation algorithm to evaluate each element
;       in Z at virtual coordinates contained in X and Y with the data
;       in P.                                                          
;
; COMMON BLOCKS:
;       If repeated interpolation of the same array is to occur, then
;       one can save time by initializing the common block RINTER.    
;
; REVISION HISTORY:
;       March 1988 written W. Landsman STX Co.
;       Checked for IDL Version 2, J. Isensee, September, 1990
;       Corrected call to HISTOGRAM, W. Landsman November 1990
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Fix output derivatives for 2-d inputs, added /INIT W. Landsman May 2000
;       
;-
 On_error, 2     
 common rinter, c1, c2, c3, init         

 if (N_params() LT 3) and (NOT keyword_set(INIT)) then begin     
        print, 'Syntax:  Z = RINTER( P, X, Y, [ DFDX, DFDY] )  '
        print, '   or    Z = RINTER( P, /INIT)    to initialize common block
        print,'P -  Array to be interpolated' 
        print,'X -  Vector or array of X positions' 
        print,'Y -  Vector or array of Y Positions'
        print,'DFDX, DFDY - Optional output derivatives '
        return,0
 endif 

 c = size(p)
 if c[0] NE 2 then $    
        message,'Input array (first parameter) must be 2 dimensional'

 if keyword_set(initialize) then begin

; Don't use SHIFT function to avoid wraparound at the end points

      nx = c[1]
      p_1 = p & p1 = p & p2 = p
      p_1[1,0] = p[0:nx-2,*]
      p1[0,0] = p[1:*,*]
      p2[0,0] = p[2:*,*]
      c1 = 0.5*(p1 - p_1)
      c2 = 2.*p1 + p_1 - 0.5*(5.*p + p2)
      c3 = 0.5*(3.*(p-p1) + p2 - p_1)
      init = 1
      if N_params() LT 3 then return,0
 endif
 
 sx = size(x)
 npts = sx[sx[0]+2]
 c[3] = c[3] > 4                        ;Make sure output array at least REAL

 i = long( x[*] ) 
 j = long( y[*] )
 xdist = x[*] - i  
 ydist = y[*] - j
 x_1 = c[1]*(j-1) + i
 x0 = x_1 + c[1] 
 x1 = x0 + c[1] 
 x2 = x1 + c[1]

 if N_elements(init) EQ 0 then init = 0    ;Has COMMON block been initialized?

 if init EQ 0 then begin 

   xgood = [ x_1,x0,x1,x2 ]
   num = histogram( xgood, MIN=0)
   xgood = where( num GE 1 ) 
   p_1 = p[xgood-1] & p0 = p[xgood] & p1 = p[xgood+1] & p2 = p[xgood+2]
   c1 = p*0.  & c2 = c1 & c3 = c1
   c1[xgood] = 0.5*( p1 - p_1)
   c2[xgood] = 2.*p1 + p_1 - 0.5*(5.*p0 + p2)
   c3[xgood] = 0.5*(3.*(p0 - p1) + p2 - p_1)
 endif

 y_1 = xdist*( xdist*( xdist*c3[x_1] +c2[x_1]) + c1[x_1]) + p[x_1]
 y0 =  xdist*( xdist*( xdist*c3[x0] +c2[x0]) + c1[x0]) + p[x0]
 y1 =  xdist*( xdist*( xdist*c3[x1] +c2[x1]) + c1[x1]) + p[x1]
 y2 =  xdist*( xdist*( xdist*c3[x2] +c2[x2]) + c1[x2]) + p[x2]

 if N_params() GT 3 then begin
 
    dy_1 = xdist*(xdist*c3[x_1]*3. + 2.*c2[x_1]) + c1[x_1]
    dy0  = xdist*(xdist*c3[x0 ]*3. + 2.*c2[x0]) + c1[x0]
    dy1  = xdist*(xdist*c3[x1 ]*3. + 2.*c2[x1]) + c1[x1]
    dy2  = xdist*(xdist*c3[x2 ]*3. + 2.*c2[x2]) + c1[x2]
    d1 = 0.5*(dy1 - dy_1)
    d2 = 2.*dy1 + dy_1 - 0.5*(5.*dy0 +dy2)
    d3 = 0.5*( 3.*( dy0-dy1 ) + dy2 - dy_1)
    dfdx =  ydist*( ydist*( ydist*d3 + d2 ) + d1 ) + dy0

 endif

 d1 = 0.5*(y1 - y_1)
 d2 = 2.*y1 + y_1 - 0.5*(5.*y0 +y2)
 d3 = 0.5*(3.*(y0-y1) + y2 - y_1)
 z =  ydist*(ydist*(ydist*d3 + d2) + d1) + y0  
 if N_params() GT 3 then dfdy = ydist*(ydist*d3*3.+2*d2) + d1   
 
 if ( sx[0] EQ 2 ) then begin        ;Convert results to 2-D if desired

    z = reform(z,sx[1],sx[2] ) 
    if N_params() GT 3 then  begin      ;Create output derivative arrays?
          dfdx = reform(dfdx,sx[1],sx[2])
          dfdy = reform(dfdy,sx[1],sx[2])          
     endif

 endif

 return,z
 end