include "register.h"

#    plyvec -- Calculates intensity at given (floating) position in vector
#    This subroutine will, given a position in the input (vector) file,
#    calculate the interpolated intensity at that location using a
#    polynomial interpolation algorithm based on that given in
#    "Data Reduction and Error Analysis for the Physical Sciences"
#    by P.R. Bevington (pp 265 ff). In this case it is assumed that
#    the data are evenly spaced, as this version assumes there is no
#    data window mask.
#
#        Version    Date        Author     Description
#                28-JAN-1993  Ray Williamson Convert to SPP
#           3    22-AUG-1983  PAT MURPHY   Generalize to allow window mask
#           2     8-AUG-1983  PAT MURPHY   Add Fortran Code
#           1     8-AUG-1983  PAT MURPHY   Prolog Generation and PDL generation

procedure plyvec (input, dim, posit, deg, badpix, numerr, out)

int dim	    	    	    	    	#  (I) The size of the input array
int deg	    	    	    	    	#  (I) The degree of polynomial (0-12)
int numerr  	    	    	    	#  (O) Running error count for NN
int x[13]   	    	    	    	#  (L) Locations of good pixels
int i, j, aj  	    	    	    	#  (L) Misc. Local variables
int status  	    	    	    	#  Status Flag
real input[dim]	    	    	    	#  (I) The actual input array
real posit  	    	    	    	#  (I) Position of point to be interp'ed
real badpix 	    	    	    	#  (I) Fill value in data mask
real out    	    	    	    	#  (O) Returned value at that point
real y[13]  	    	    	    	#  (L) Values of good pixels
real prod, delta    	    	    	#  (L) Misc. Local variables
int isgn ()	    	    	    	#  External sign function
double factor ()    	    	    	#  Factorial function

begin

#  If there aren't enough points for fit
#  or if factorial will overflow, then quit

        if ((deg+1) > dim || deg    > 12) {
        call error(1,
	    "Not enough points for fit or possible factorial OVERFLOW")

        } else if (posit < 1. || posit > float(dim)) {

          out = badpix

        } else {

#       We're ok so extract the x and y vectors for the interpolation

          call expixv (input, dim, posit, deg+1, x, y, status)

#  If problems, try Nearest N.
          if (status != OK ) goto 10

#       We assume that there is no mask, so we have even spacing

#  Initialize output value = A[1]
          out = y[1]

#  If the user tries to do zeroth order
          if (deg <= 0) {
#  polynomial, use Nearest Neighbor
            status = INFO
#  instead and give warning message when
            goto 10
#  we're through with entire vector.
          }

#  Get difference; Note X[2]-X[1] is 1
          delta = posit - x[1]
#  For each additional 'A' term
          do j = 2, deg + 1 {
#  Initialize the product
            prod = 1.
#  and the sum
            aj = 0.
#  This is the A[j] calculation loop
            do i = 1, j {
#Note Product loop only goes to J-1
              if (i < j) prod = prod * (delta - i + 1)

#   ((-1)**(j-i))*y[i]
              aj = aj + isgn (j - i) * y[i] /
#	    	(i-1)!*(j-i)!
	    	(int(factor (i - 1)) * int(factor (j - i)))
            }
            out = out + aj * prod
          }
        }

#       If we got this far with Status = ok, we're through.
#       iF NOT, WE WANT TO TRY NEAREST NEIGHBOR.
#       Use NUMERR to keep track of # of Nearest Neighbor interpolations done. 
10	if (status != OK) {

            call nnbvec (input, dim, posit, badpix, out)
            numerr = numerr + 1

        }


end