Repatch (from linux) of OSX IRAF

1 files changed, 294 insertions, 0 deletions

diff --git a/math/iminterp/msider.x b/math/iminterp/msider.x
new file mode 100644
index 00000000..e66d9119
--- /dev/null
+++ b/math/iminterp/msider.x
@@ -0,0 +1,294 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "im2interpdef.h"
+include <math/iminterp.h>
+
+# MSIDER -- Calculate the derivatives of the interpolant. The derivative
+# der[i,j] = d f(x,y) / dx (i-1) dy (j-1). Therefore der[1,1] contains
+# the value of the interpolant, der[2,1] the 1st derivative in x and
+# der[1,2] the first derivative in y.
+
+procedure msider (msi, x, y, der, nxder, nyder, len_der)
+
+pointer	msi			# pointer to interpolant descriptor structure
+real	x[ARB]			# x value
+real	y[ARB]			# y value
+real	der[len_der,ARB]	# derivative array
+int	nxder			# number of x derivatives
+int	nyder			# number of y derivatives
+int	len_der			# row length of der, len_der >= nxder
+
+int	first_point, len_coeff
+int	nxterms, nyterms, nx, ny, nyd, nxd
+int	i, j, ii, jj
+real	sx, sy, tx, ty, xmin, xmax, ymin, ymax
+real	pcoeff[MAX_NTERMS,MAX_NTERMS], pctemp[MAX_NTERMS,MAX_NTERMS]
+real	sum[MAX_NTERMS], accum, deltax, deltay, tmpx[4], tmpy[4]
+pointer	index, ptr
+
+begin
+	if (nxder < 1 || nyder < 1)
+	    return
+
+	# set up coefficient array parameters
+	len_coeff = MSI_NXCOEFF(msi)
+	index = MSI_COEFF(msi) + MSI_FSTPNT(msi) - 1
+
+	# zero the derivatives
+	do j = 1, nyder {
+	    do i = 1, nxder
+		der[i,j] = 0.
+	}
+
+	# calculate the appropriate number of terms of the polynomials in
+	# x and y
+
+	switch (MSI_TYPE(msi)) {
+
+	case II_BINEAREST:
+
+	    nx = x[1] + 0.5
+	    ny = y[1] + 0.5
+
+	    ptr = index + (ny - 1) * len_coeff + nx
+	    der[1,1] = COEFF(ptr) 
+
+	    return
+
+	case II_BISINC, II_BILSINC:
+
+	    call ii_bisincder (x[1], y[1], der, nxder, nyder, len_der,
+	        COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi), MSI_NXCOEFF(msi),
+		MSI_NYCOEFF(msi), MSI_NSINC(msi), DX, DY)
+
+	    return
+
+	case II_BILINEAR:
+
+	    nx = x[1]
+	    ny = y[1]
+	    sx = x[1] - nx
+	    sy = y[1] - ny
+	    tx = 1. - sx
+	    ty = 1. - sy
+
+	    ptr = index + (ny - 1) * len_coeff + nx
+	    der[1,1] = tx * ty * COEFF(ptr) + sx * ty * COEFF(ptr+1) +
+		       sy * tx * COEFF(ptr+len_coeff) +
+		       sx * sy * COEFF(ptr+len_coeff+1)
+
+	    if (nxder > 1)
+		der[2,1] = -ty * COEFF(ptr) + ty * COEFF(ptr+1) -
+			   sy * COEFF(ptr+len_coeff) +
+			   sy * COEFF(ptr+len_coeff+1)
+
+	    if (nyder > 1)
+		der[1,2] = -tx * COEFF(ptr) - sx * COEFF(ptr+1) +
+			   tx * COEFF(ptr+len_coeff) +
+			   sx * COEFF(ptr+len_coeff+1)
+
+	    if (nyder > 1 && nxder > 1)
+		der[2,2] = COEFF(ptr) - COEFF(ptr+1) - COEFF(ptr+len_coeff) +
+			   COEFF(ptr+len_coeff+1)
+
+	    return
+
+	case II_BIDRIZZLE:
+            if (MSI_XPIXFRAC(msi) >= 1.0 && MSI_YPIXFRAC(msi) >= 1.0)
+                call ii_bidriz1 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                    MSI_NXCOEFF(msi), x, y, der[1,1], 1, MSI_BADVAL(msi))
+            #else if (MSI_XPIXFRAC(msi) <= 0.0 && MSI_YPIXFRAC(msi) <= 0.0)
+                #call ii_bidriz0 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                    #MSI_NXCOEFF(msi), x, y, der[1,1], 1, MSI_BADVAL(msi))
+            else
+                call ii_bidriz (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                    MSI_NXCOEFF(msi), x, y, der[1,1], 1, MSI_XPIXFRAC(msi),
+                    MSI_YPIXFRAC(msi), MSI_BADVAL(msi))
+
+	    if (nxder > 1) {
+		xmax = max (x[1], x[2], x[3], x[4])
+		xmin = min (x[1], x[2], x[3], x[4])
+		ymax = max (y[1], y[2], y[3], y[4])
+		ymin = min (y[1], y[2], y[3], y[4])
+		deltax = xmax - xmin
+		if (deltax == 0.0)
+		    der[2,1] = 0.0
+		else {
+		    tmpx[1] = xmin; tmpy[1] = ymin
+		    tmpx[2] = (xmax - xmin) / 2.0; tmpy[2] = ymin
+		    tmpx[3] = (xmax - xmin) / 2.0; tmpy[3] = ymax
+		    tmpx[4] = xmin; tmpy[4] = ymax
+            	    if (MSI_XPIXFRAC(msi) >= 1.0 && MSI_YPIXFRAC(msi) >= 1.0)
+                        call ii_bidriz1 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, accum, 1,
+			    MSI_BADVAL(msi))
+            	    #else if (MSI_XPIXFRAC(msi) <= 0.0 &&
+		        #MSI_YPIXFRAC(msi) <= 0.0)
+                        #call ii_bidriz0 (COEFF(MSI_COEFF(msi)),
+			    #MSI_FSTPNT(msi), MSI_NXCOEFF(msi), tmpx, tmpy,
+			    #accum, 1, MSI_BADVAL(msi))
+            	    else
+                        call ii_bidriz (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, accum, 1,
+			    MSI_XPIXFRAC(msi), MSI_YPIXFRAC(msi),
+			    MSI_BADVAL(msi))
+		    tmpx[1] = (xmax - xmin) / 2.0; tmpy[1] = ymin
+		    tmpx[2] = xmax; tmpy[2] = ymin
+		    tmpx[3] = xmax; tmpy[3] = ymax
+		    tmpx[4] = (xmax - xmin) / 2.0; tmpy[4] = ymax
+            	    if (MSI_XPIXFRAC(msi) >= 1.0 && MSI_YPIXFRAC(msi) >= 1.0)
+                        call ii_bidriz1 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, der[2,1], 1,
+			    MSI_BADVAL(msi))
+            	    #else if (MSI_XPIXFRAC(msi) <= 0.0 &&
+		        #MSI_YPIXFRAC(msi) <= 0.0)
+                        #call ii_bidriz0 (COEFF(MSI_COEFF(msi)),
+			    #MSI_FSTPNT(msi), MSI_NXCOEFF(msi), tmpx, tmpy,
+			    #der[2,1], 1, MSI_BADVAL(msi))
+            	    else
+                        call ii_bidriz (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, der[2,1], 1,
+			    MSI_XPIXFRAC(msi), MSI_YPIXFRAC(msi),
+			    MSI_BADVAL(msi))
+		    der[2,1] = 2.0 * (der[2,1] - accum) / deltax
+		}
+	    }
+	    if (nyder > 1) {
+		deltay = ymax - ymin
+		if (deltay == 0.0)
+		    der[1,2] = 0.0
+		else {
+		    tmpx[1] = xmin; tmpy[1] =  ymin
+		    tmpx[2] = xmax; tmpy[2] = ymin
+		    tmpx[3] = xmax; tmpy[3] = (ymax - ymin) / 2.0
+		    tmpx[4] = xmin; tmpy[4] = (ymax - ymin) / 2.0
+            	    if (MSI_XPIXFRAC(msi) >= 1.0 && MSI_YPIXFRAC(msi) >= 1.0)
+                        call ii_bidriz1 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, accum, 1,
+			    MSI_BADVAL(msi))
+            	    #else if (MSI_XPIXFRAC(msi) <= 0.0 &&
+		        #MSI_YPIXFRAC(msi) <= 0.0)
+                        #call ii_bidriz0 (COEFF(MSI_COEFF(msi)),
+			    #MSI_FSTPNT(msi), MSI_NXCOEFF(msi), tmpx, tmpy,
+			    #accum, 1, MSI_BADVAL(msi))
+            	    else
+                        call ii_bidriz (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, accum, 1,
+			    MSI_XPIXFRAC(msi), MSI_YPIXFRAC(msi),
+			    MSI_BADVAL(msi))
+		    tmpx[1] = xmin; tmpy[1] = (ymax - ymin) /  2.0
+		    tmpx[2] = xmax; tmpy[2] = (ymax - ymin) / 2.0
+		    tmpx[3] = xmax; tmpy[3] = ymax
+		    tmpx[4] = xmin; tmpy[4] = ymax
+            	    if (MSI_XPIXFRAC(msi) >= 1.0 && MSI_YPIXFRAC(msi) >= 1.0)
+                        call ii_bidriz1 (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, der[1,2], 1,
+			    MSI_BADVAL(msi))
+            	    #else if (MSI_XPIXFRAC(msi) <= 0.0 &&
+		        #MSI_YPIXFRAC(msi) <= 0.0)
+                        #call ii_bidriz0 (COEFF(MSI_COEFF(msi)),
+			    #MSI_FSTPNT(msi), MSI_NXCOEFF(msi), tmpx, tmpy,
+			    #der[1,2], 1, MSI_BADVAL(msi))
+            	    else
+                        call ii_bidriz (COEFF(MSI_COEFF(msi)), MSI_FSTPNT(msi),
+                            MSI_NXCOEFF(msi), tmpx, tmpy, der[1,2], 1,
+			    MSI_XPIXFRAC(msi), MSI_YPIXFRAC(msi),
+			    MSI_BADVAL(msi))
+		    der[1,2] = 2.0 * (der[1,2] - accum) / deltay
+		}
+	    }
+
+	    return
+
+	case II_BIPOLY3:
+
+	    nxterms = 4
+	    nyterms = 4
+
+	    nxd = min (nxder, 4)
+	    nyd = min (nyder, 4) 
+
+	    nx = x[1]
+	    sx = x[1] - nx
+	    ny = y[1]
+	    sy = y[1] - ny
+
+	    first_point = MSI_FSTPNT(msi) + (ny - 2) * len_coeff + nx
+	    call ii_pcpoly3 (COEFF(MSI_COEFF(msi)), first_point, len_coeff,
+	    		     pcoeff, 6)
+
+	case II_BIPOLY5:
+
+	    nxterms = 6
+	    nyterms = 6
+
+	    nxd = min (nxder, 6)
+	    nyd = min (nyder, 6)
+
+	    nx = x[1]
+	    sx = x[1] - nx
+	    ny = y[1]
+	    sy = y[1] - ny
+
+	    first_point = MSI_FSTPNT(msi) + (ny - 3) * len_coeff + nx
+	    call ii_pcpoly5 (COEFF(MSI_COEFF(msi)), first_point, len_coeff,
+	    		    pcoeff, 6)
+
+	case II_BISPLINE3:
+
+	    nxterms = 4
+	    nyterms = 4
+
+	    nxd = min (nxder, 4)
+	    nyd = min (nyder, 4)
+
+	    nx = x[1]
+	    sx = x[1] - nx
+	    ny = y[1]
+	    sy = y[1] - ny
+
+	    first_point = MSI_FSTPNT(msi) + (ny - 2) * len_coeff + nx
+	    call ii_pcspline3 (COEFF(MSI_COEFF(msi)), first_point, len_coeff,
+	    			pcoeff, 6)
+	}
+
+	# evaluate the derivatives by nested multiplication
+	do j = 1, nyd {
+
+	    # set pctemp
+	    do jj =  nyterms, j, -1 {
+		do ii = 1, nxterms
+		    pctemp[ii,jj] = pcoeff[ii,jj]
+	    }
+
+	    do i = 1, nxd {
+
+		# accumulate the partial sums in x
+		do jj = nyterms, j, -1 {
+		    sum[jj] = pctemp[nxterms,jj]
+		    do ii = nxterms - 1, i, -1
+			sum[jj] = pctemp[ii,jj] + sum[jj] * sx
+		}
+
+		# accumulate the sum in y
+		accum = sum[nyterms]
+		do jj = nyterms - 1, j, -1
+		    accum = sum[jj] + accum * sy
+
+		# evaluate derivative
+		der[i,j] = accum
+
+		# differentiate in x
+		do jj = nyterms, j, -1 {
+		    do ii = nxterms, i + 1, -1
+		        pctemp[ii,jj] = (ii - i) * pctemp[ii,jj]
+		}
+	    }
+
+	    # differentiate in y
+	    do jj = 1, nxterms {
+		do ii = nyterms, j + 1, -1
+		    pcoeff[jj,ii] = (ii - j) * pcoeff[jj,ii]
+	    }
+	}
+end