include "register.h"

#    plyvcw -- Calculates intensity at given (floating) position in vector
# 
#    Description:
#    ------------
#    This subroutine will, given a position in the input (vector) file,
#    calculate the interpolated intensity at that location using a
#    polynomial interpolation algorithm as given in Bevington (see the
#    description of the module without data mask checking for reference).
#    This version checks the data window mask and uses the general
#    algorithm where uneven spacing is assumed.
#
#        Version    Date        Author     Description
#                28-JAN-1993  Ray Williamson Conert 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 plyvcw (input, dim, posit, deg, wndo, badpix, numerr, out)

int dim	    	    	    	    	#  The size of the input array
int deg	    	    	    	    	#  The degree of polynomial (0-12)
int numerr  	    	    	    	#  Number of failed poly interps.
int x[13]   	    	    	    	#  Locations of good pixels
int i,j,k   	    	    	    	#  Local variables
int status  	    	    	    	#  The internal status, normally 'OK'
real input[dim]	    	    	    	#  The actual input array
real wndo[dim]	    	    	    	#  Window mask array
real posit  	    	    	    	#  Position of point to be interpolated
real out    	    	    	    	#  Returned value at that point
real badpix 	    	    	    	#  Fill value in data mask
real y[13]   	    	    	    	#  Values of good pixels
real a[13]  	    	    	    	#  Polynomial parameters
real prod, delta, dl[13]    	    	#  Local variables
real denom, sum	    	    	    	#  More of the same

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)) {

#  We interpolate, not extrapolate!
          out = badpix

        } else {

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

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

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

#       This routine assumes the presence of a mask so noneven spacing assumed

#  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.
          }

#  Denominator for delta calculations
          denom = x[2] - x[1]
#  This is the main or first delta
          delta = (posit - x[1])/denom
#  Now calculate the rest of the deltas
          do i = 2, deg + 1 {
#  and store them for future use.
            dl[i] = (x[i] - x[1])/denom
          }

#         Now calculate the a[j] polynomial coefficients

#  Set the first one
          a[1] = y[1]
#  For each additional coefficient
          do j = 2, deg + 1 {
#  Initialize the product of delta diffs
            prod = 1.
#  and the sum which will be the a[j]
            sum = 0.
#  The algorithm used here is adapted
            do i = 1, j - 1 {
#  from Bevington, p. 267 (prog. 13-2)
              k = j - i
              prod = prod * (dl[j] - dl[k])
              sum = sum - a[k]/prod
            }
            a[j] = sum + y[j]/prod
          }

#         We now have all the A coefficients calculated so let us interpolate

          out = a[1]
          do j = 2, deg + 1 {
            prod = 1.
            do i = 1, j - 1 {
              prod = prod * (delta - dl[i])
            }
            out = out + a[j] * prod
          }
        }

#       Finally, if we could not do polynomial, try nearest neighbor

10	if (status != OK) {

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

        }


end