#    pt4img -- Performs 3, 4, or 6 point interpolation on an image
# 
#    Description:
#    ------------
#    This module will perform four point interpolation at a given point
#    in an image.
#
#     Version    Date        Author     Description
#             28-JAN-1993  RAY WILLIAMSON Port to STSDAS
#        6    12-APR-1985  D.L.BALL     Improve efficiency of algorithm
#        5    10-FEB-1984  Pat Murphy   Add more fortran
#        3    19-JAN-1984  Pat Murphy   Initial FORTRAN coding
#        2    19-JAN-1984  Pat Murphy   Add PDL
#        1    19-JAN-1984  PAT MURPHY   Prolog Generation

procedure pt4img (x, y, input, dim1,dim2, badpix, output)

int dim1, dim2	    	    	    	    #  Size of input array
int ix, iy, ixplus1, iyplus1	    	    #  Integer pixel locations
real input[dim1,dim2]	    	    	    #  Input array
real x, y   	    	    	    	    #  co-ords of point to interp at
real badpix 	    	    	    	    #  Fill value if we can't interp
real output 	    	    	    	    #  result
real p, q, oneminusp, oneminusq	    	    #  temporary variables

begin
        if (x < 1 || x > dim1 || y < 1 || y > dim2) {

            output = badpix

        } else {

#  Get integer location of lower left
#  pixel nearest given co-ord.
            ix = int[x]
            iy = int[y]

# -------------- Added 840620 PPM: catch exact correspondance -----------------

            if (x == float(dim1)) ix = ix - 1
            if (y == float(dim2)) iy = iy - 1

# ------------------------ End of 840620 changes ------------------------------

#  Store 1 + these for efficiency
            ixplus1 = ix + 1
            iyplus1 = iy + 1
#  Get fractional part of co-ords
            p = x - ix
#  in original image
            q = y - iy

#  store these too as they are used
            oneminusp = 1.0 - p
#  more than once in the calculation
            oneminusq = 1.0 - q

            output = oneminusq *(oneminusp * input [ix,      iy]+
                               p         * input [ixplus1, iy])+
                   q         *(oneminusp * input [ix,      iyplus1]+
                               p         * input [ixplus1, iyplus1])

        }

end