From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- math/iminterp/ii_1dinteg.x | 372 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 372 insertions(+) create mode 100644 math/iminterp/ii_1dinteg.x (limited to 'math/iminterp/ii_1dinteg.x') 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 + +# 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 -- cgit