Initial commit

1 files changed, 372 insertions, 0 deletions

diff --git a/math/iminterp/ii_1dinteg.x b/math/iminterp/ii_1dinteg.x
new file mode 100644
index 00000000..05b80542
--- /dev/null
+++ b/math/iminterp/ii_1dinteg.x
@@ -0,0 +1,372 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "im1interpdef.h"
+include <math/iminterp.h>
+
+# II_1DINTEG -- Find the integral of the interpolant from a to b be assuming
+# that both a and b land in the array. This routine is not used directly
+# in the 1D interpolation package but is actually called repeatedly from the
+# 2D interpolation package. Therefore the SINC function interpolator has
+# not been implemented, although it has been blocked in.
+
+real procedure ii_1dinteg (coeff, len_coeff, a, b, interp_type, nsinc, dx,
+	pixfrac)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	len_coeff	# length of coefficient array (used in sinc only)
+real	a		# lower limit for integral
+real	b		# upper limit for integral
+int	interp_type	# type of 1D interpolant
+int	nsinc		# width of sinc function
+real	dx		# sinc precision
+real	pixfrac		# drizzle pixel fraction
+
+int	neara, nearb, i, j, nterms
+real	deltaxa, deltaxb, accum, xa, xb, pcoeff[MAX_NDERIVS]
+
+begin
+	# Flip order and sign at end.
+	xa = a
+	xb = b
+	if (a > b) {
+	    xa = b
+	    xb = a
+	}
+
+	# Initialize.
+	neara = xa
+	nearb = xb
+	accum = 0.
+
+	switch (interp_type) {
+	case II_NEAREST:
+	    nterms = 0
+	case II_LINEAR:
+	    nterms = 1
+	case II_DRIZZLE:
+	    nterms = 0
+	case II_POLY3:
+	    nterms = 4
+	case II_POLY5:
+	    nterms = 6
+	case II_SPLINE3:
+	    nterms = 4
+	case II_SINC, II_LSINC:
+	    nterms = 0
+	}
+
+	switch (interp_type) {
+	# NEAREST
+	case II_NEAREST:
+
+	    # Reset segment to center values.
+	    neara = xa + 0.5
+	    nearb = xb + 0.5
+
+	    # Set up for first segment.
+	    deltaxa = xa - neara
+
+	    # For clarity one segment case is handled separately.
+
+	    # Only one segment involved.
+	    if (nearb == neara) {
+
+		deltaxb = xb - nearb
+		accum = accum + (deltaxb - deltaxa) * coeff[neara]
+
+	    # More than one segment.
+	    } else {
+
+		# First segment.
+		accum = accum + (0.5 - deltaxa) * coeff[neara]
+
+		# Middle segment.
+		do j = neara + 1, nearb - 1 {
+		    accum = accum + coeff[j]
+		}
+
+		# Last segment.
+		deltaxb = xb - nearb
+		accum = accum + (deltaxb + 0.5) * coeff[nearb]
+	    }
+
+	# LINEAR
+	case II_LINEAR:
+
+	    # Set up for first segment.
+	    deltaxa = xa - neara
+
+	    # For clarity one segment case is handled separately.
+
+	    # Only one segment involved.
+	    if (nearb == neara) {
+
+		deltaxb = xb - nearb
+		accum = accum + (deltaxb - deltaxa) * coeff[neara] +
+		     0.5 * (coeff[neara+1] - coeff[neara]) *
+		     (deltaxb * deltaxb - deltaxa * deltaxa)
+
+	    # More than one segment.
+	    } else {
+
+		# First segment.
+		accum = accum + (1. - deltaxa) * coeff[neara] +
+		     0.5 * (coeff[neara+1] - coeff[neara]) *
+		     (1. - deltaxa * deltaxa)
+
+		# Middle segment.
+		do j = neara + 1, nearb - 1 {
+		    accum = accum + 0.5 * (coeff[j+1] + coeff[j])
+		}
+
+		# Last segment.
+		deltaxb = xb - nearb
+		accum = accum + coeff[nearb] * deltaxb + 0.5 *
+			(coeff[nearb+1] - coeff[nearb]) * deltaxb * deltaxb
+	    }
+
+	# DRIZZLE-- Note that to get get pixfrac an interface change was
+	# required.
+	case II_DRIZZLE:
+	    if (pixfrac >= 1.0)
+		call ii_dzigrl1 (a, b, accum, coeff)
+	    else
+		call ii_dzigrl (a, b, accum, coeff, pixfrac)
+
+	# SINC -- Note that to get ncoeff, nsinc, and dx,  an interface change
+	# was required. 
+	case II_SINC, II_LSINC:
+	    call ii_sincigrl (xa, xb, accum, coeff, len_coeff, nsinc, dx)
+
+	# A higher order interpolant.
+	default:
+
+	    # Set up for first segment.
+	    deltaxa = xa - neara
+
+	    # For clarity one segment case is handled separately.
+
+	    # Only one segment involved.
+	    if (nearb == neara) {
+
+		deltaxb = xb - nearb
+		call ii_getpcoeff (coeff, neara, pcoeff, interp_type)
+		do i = 1, nterms
+		    accum = accum + (1./i) * pcoeff[i] *
+		    	    (deltaxb ** i - deltaxa ** i)
+
+	    # More than one segment.
+	    } else {
+
+		# First segment.
+		call ii_getpcoeff (coeff, neara, pcoeff, interp_type)
+		do i = 1, nterms
+		    accum = accum + (1./i) * pcoeff[i] * (1. - deltaxa ** i)
+
+		# Middle segment.
+		do j = neara + 1, nearb - 1 {
+		    call ii_getpcoeff (coeff, j, pcoeff, interp_type)
+
+		    do i = 1, nterms
+			accum = accum + (1./i) * pcoeff[i]
+		}
+
+		# Last segment.
+		deltaxb = xb - nearb
+		call ii_getpcoeff (coeff, nearb, pcoeff, interp_type)
+		do i = 1, nterms
+		    accum = accum + (1./i) * pcoeff[i] * deltaxb ** i
+	    }
+	}
+
+	if (a < b)
+	    return (accum)
+	else
+	    return (-accum)
+end
+
+	
+# II_GETPCOEFF -- Generates polynomial coefficients if the interpolant is
+# SPLINE3, POLY3 or POLY5.
+
+procedure ii_getpcoeff (coeff, index, pcoeff, interp_type)
+
+real	coeff[ARB]	# coefficient array
+int	index		# coefficients wanted for index < x < index + 1
+real	pcoeff[ARB]	# polynomial coefficients
+int	interp_type	# type of interpolant
+
+int	i, k, nterms
+real	diff[MAX_NDERIVS]
+
+begin
+	# generate polynomial coefficients, first for spline
+
+	if (interp_type == II_SPLINE3) {
+
+	    pcoeff[1] = coeff[index-1] + 4. * coeff[index] + coeff[index+1]
+	    pcoeff[2] = 3. * (coeff[index+1] - coeff[index-1])
+	    pcoeff[3] = 3. * (coeff[index-1] - 2. * coeff[index] +
+	    		coeff[index+1])
+	    pcoeff[4] = -coeff[index-1] + 3. * coeff[index] -
+	    		3. * coeff[index+1] + coeff[index+2]
+	} else {
+
+	    if (interp_type == II_POLY5)
+		nterms = 6
+
+	    # must be POLY3
+	    else
+		nterms = 4
+
+	    # Newton's form written in line to get polynomial from data
+
+	    # load data
+	    do i = 1, nterms
+		diff[i] = coeff[index - nterms/2 + i]
+
+	    # generate difference table
+	    do k = 1, nterms - 1
+		do i = 1, nterms - k
+		    diff[i] = (diff[i+1] - diff[i]) / k
+
+	    # shift to generate polynomial coefficients of (x - index)
+	    do k = nterms, 2, -1
+		do i = 2, k
+		    diff[i] = diff[i] + diff[i-1] * (k - i - nterms/2)
+
+	    do i = 1, nterms
+		pcoeff[i] = diff[nterms + 1 - i]
+	}
+end
+
+
+# II_SINCIGRL -- Evaluate integral of sinc interpolator.
+# The integral is computed by dividing interval into a number of equal
+# size subintervals which are at most one pixel wide.  The integral
+# of each subinterval is the central value times the interval width.
+
+procedure ii_sincigrl (a, b, sum, data, npix, nsinc, mindx)
+
+real	a, b			# integral limits
+real	sum			# output integral value
+real	data[npix]		# input data array
+int	npix			# number of pixels
+int	nsinc			# sinc truncation length
+real	mindx			# interpolation minimum
+
+int	n
+real	x, y, dx, x1, x2
+
+begin
+	x1 = min (a, b)
+	x2 = max (a, b)
+	n = max (1, nint (x2 - x1))
+	dx = (x2 - x1) / n
+
+	sum = 0.
+	for (x = x1 + dx / 2; x < x2; x = x + dx) {
+	    call ii_sinc (x, y, 1, data, npix, nsinc, mindx)
+	    sum = sum + y * dx
+	}
+end
+
+
+# II_DZIGRL -- Procedure to integrate the drizzle interpolant.
+
+procedure ii_dzigrl (a, b, sum, data, pixfrac)
+
+real    a, b         	# x start and stop values, must be within [1,npts]
+real    sum		# integgral value returned to the user
+real    data[ARB]       # data to be interpolated
+real    pixfrac         # the drizzle pixel fraction
+
+int     j, neara, nearb
+real    hpixfrac, xa, xb, dx, accum
+
+begin
+        hpixfrac = pixfrac / 2.0
+
+        # Define the interval of integration.
+        xa = min (a, b)
+        xb = max (a, b)
+        neara = xa + 0.5
+        nearb = xb + 0.5
+
+        # Initialize the integration
+        accum = 0.0
+        if (neara == nearb) {
+
+            dx = min (xb, nearb + hpixfrac) - max (xa, neara - hpixfrac)
+            if (dx > 0.0)
+                accum = accum + dx * data[neara]
+
+        } else {
+
+            # first segement
+            dx = neara + hpixfrac - max (xa, neara - hpixfrac)
+            if (dx > 0.0)
+                accum = accum + dx * data[neara]
+
+            # interior segments.
+            do j = neara + 1, nearb - 1
+                accum = accum + pixfrac * data[j]
+
+            # last segment
+            dx = min (xb, nearb + hpixfrac) - (nearb - hpixfrac)
+            if (dx > 0.0)
+                accum = accum + dx * data[nearb]
+        }
+
+        if (a > b)
+            sum = -accum
+        else
+            sum = accum 
+end
+
+
+# II_DZIGRL1 -- Procedure to integrate the drizzle interpolant in the case
+# where pixfrac = 1.0.
+
+procedure ii_dzigrl1 (a, b, sum, data)
+
+real    a, b            # x start and stop values, must be within [1,npts]
+real    sum		# integgral value returned to the user
+real    data[ARB]       # data to be interpolated
+
+int     j, neara, nearb
+real    xa, xb, deltaxa, deltaxb, accum
+
+begin
+        # Define the interval of integration.
+        xa = min (a, b)
+        xb = max (a, b)
+        neara = xa + 0.5
+        nearb = xb + 0.5
+        deltaxa = xa - neara
+        deltaxb = xb - nearb
+
+        # Only one segment involved.
+        accum = 0.0
+        if (neara == nearb) {
+
+            accum = accum + (deltaxb - deltaxa) * data[neara]
+
+        } else {
+
+            # First segment.
+            accum = accum + (0.5 - deltaxa) * data[neara]
+
+            # Middle segment.
+            do j = neara + 1, nearb - 1
+                accum = accum + data[j]
+
+            # Last segment.
+            accum = accum + (deltaxb + 0.5) * data[nearb]
+        }
+
+        if (a > b)
+           sum = -accum 
+        else
+           sum = accum
+end