1 files changed, 96 insertions, 0 deletions
diff --git a/math/iminterp/msisqgrl.x b/math/iminterp/msisqgrl.x
new file mode 100644
index 00000000..82a3122b
--- /dev/null
+++ b/math/iminterp/msisqgrl.x
@@ -0,0 +1,96 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <mach.h>
+include "im2interpdef.h"
+include <math/iminterp.h>
+
+# MSISQGRL -- Procedure to integrate the 2D interpolant over a rectangular
+# region.
+
+real procedure msisqgrl (msi, x1, x2, y1, y2)
+
+pointer	msi		# pointer to the interpolant descriptor structure
+real	x1, x2		# x integration limits
+real	y1, y2		# y integration limits
+
+int	i, interp_type, nylmin, nylmax, offset 
+pointer	xintegrl, ptr
+real	xmin, xmax, ymin, ymax, accum
+real	ii_1dinteg()
+
+begin
+	# set up 1D interpolant type to match 2-D interpolant
+	switch (MSI_TYPE(msi)) {
+	case II_BINEAREST:
+	    interp_type = II_NEAREST
+	case II_BILINEAR:
+	    interp_type = II_LINEAR
+	case II_BIDRIZZLE:
+	    interp_type = II_DRIZZLE
+	case II_BIPOLY3:
+	    interp_type = II_POLY3
+	case II_BIPOLY5:
+	    interp_type = II_POLY5
+	case II_BISPLINE3:
+	    interp_type = II_SPLINE3
+	case II_BISINC:
+	    interp_type = II_SINC
+	case II_BILSINC:
+	    interp_type = II_LSINC
+	}
+
+	# set up temporary storage for x integrals
+	call calloc (xintegrl, MSI_NYCOEFF(msi), TY_REAL)
+
+	# switch order of x integration at the end
+	xmin = x1
+	xmax = x2
+	if (x2 < x1) {
+	    xmax = x2
+	    xmin = x1
+	}
+
+	# switch order of y integration at end
+	ymin = y1
+	ymax = y2
+	if (y2 < y1) {
+	    ymax = y2
+	    ymin = y1
+	}
+
+	# find the appropriate range in y in the coeff array
+	offset = mod (MSI_FSTPNT(msi), MSI_NXCOEFF(msi)) 
+	nylmin = max (1, min (MSI_NYCOEFF(msi), int (ymin + 0.5))) 
+	nylmax = min (MSI_NYCOEFF(msi), max (1, int (ymax + 0.5)))
+	nylmin = nylmin + offset
+	nylmax = nylmax + offset
+
+	# integrate in x
+	ptr = MSI_COEFF(msi) + offset + (nylmin - 1) * MSI_NXCOEFF(msi)
+	do i = nylmin, nylmax {
+	    Memr[xintegrl+i-1] = ii_1dinteg (COEFF(ptr), MSI_NXCOEFF(msi),
+	        xmin, xmax, interp_type, MSI_NSINC(msi), DX, MSI_XPIXFRAC(msi))
+	    if (x2 < x1)
+		Memr[xintegrl+i-1] = - Memr[xintegrl+i-1]
+	    ptr = ptr + MSI_NXCOEFF(msi)
+	}
+
+	# integrate in y
+	if (interp_type == II_SPLINE3) {
+	    call amulkr (Memr[xintegrl], 6.0, Memr[xintegrl],
+	        MSI_NYCOEFF(msi))
+	    accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi),
+	        ymin, ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
+	} else
+	    accum = ii_1dinteg (Memr[xintegrl+offset], MSI_NYCOEFF(msi),
+	        ymin, ymax, II_NEAREST, MSI_NSINC(msi), DY, MSI_YPIXFRAC(msi))
+
+	# free space
+	call mfree (xintegrl, TY_REAL)
+
+	# correct for integration error.
+	if (y2 < y1)
+	    return (-accum)
+	else
+	    return (accum)
+end