diff options
Diffstat (limited to 'math/iminterp/ii_pc2deval.x')
| -rw-r--r-- | math/iminterp/ii_pc2deval.x | 444 |
1 files changed, 444 insertions, 0 deletions
diff --git a/math/iminterp/ii_pc2deval.x b/math/iminterp/ii_pc2deval.x new file mode 100644 index 00000000..c26d2095 --- /dev/null +++ b/math/iminterp/ii_pc2deval.x @@ -0,0 +1,444 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math.h> + +# II_PCPOLY3 -- Procedure to evaluate the polynomial coefficients +# of third order in x and y using Everetts formuala. + +procedure ii_pcpoly3 (coeff, index, len_coeff, pcoeff, len_pcoeff) + +real coeff[ARB] # 1D array of interpolant coeffcients +int index # pointer into coeff array +int len_coeff # row length of coeffcients +real pcoeff[len_pcoeff,ARB] # polynomial coefficients +int len_pcoeff # row length of pcoeff + +int tptr +int i, j +real cd20, cd21, temp[4] + +begin + # determine polynomial coefficients in x + tptr = index + do i = 1, 4 { + + # calculate the central differences + cd20 = 1./6. * (coeff[tptr+1] - 2. * coeff[tptr] + coeff[tptr-1]) + cd21 = 1./6. * (coeff[tptr+2] - 2. * coeff[tptr+1] + coeff[tptr]) + + # calculate the polynomial coefficients in x at each y + pcoeff[1,i] = coeff[tptr] + pcoeff[2,i] = coeff[tptr+1] - coeff[tptr] - 2. * cd20 - cd21 + pcoeff[3,i] = 3. * cd20 + pcoeff[4,i] = cd21 - cd20 + + tptr = tptr + len_coeff + } + + # calculate polynomial coefficients in y + do j = 1, 4 { + + # calculate the central differences + cd20 = 1./6. * (pcoeff[j,3] - 2. * pcoeff[j,2] + pcoeff[j,1]) + cd21 = 1./6. * (pcoeff[j,4] - 2. * pcoeff[j,3] + pcoeff[j,2]) + + # calculate the final coefficients + temp[1] = pcoeff[j,2] + temp[2] = pcoeff[j,3] - pcoeff[j,2] - 2. * cd20 - cd21 + temp[3] = 3. * cd20 + temp[4] = cd21 - cd20 + + do i = 1, 4 + pcoeff[j,i] = temp[i] + } +end + + +# II_PCPOLY5 -- Procedure to evaluate the polynomial coefficients +# of fifth order in x and y using Everetts formuala. + +procedure ii_pcpoly5 (coeff, index, len_coeff, pcoeff, len_pcoeff) + +real coeff[ARB] # 1D array of interpolant coeffcients +int index # pointer into coeff array +int len_coeff # row length of coeffcients +real pcoeff[len_pcoeff,ARB] # polynomial coefficients +int len_pcoeff # row length of pcoeff array + +int tptr +int i, j +real cd20, cd21, cd40, cd41, temp[6] + +begin + # determine polynomial coefficients in x + tptr = index + do i = 1, 6 { + + # calculate the central differences + cd20 = 1./6. * (coeff[tptr+1] - 2. * coeff[tptr] + coeff[tptr-1]) + cd21 = 1./6. * (coeff[tptr+2] - 2. * coeff[tptr+1] + coeff[tptr]) + cd40 = 1./120. * (coeff[tptr-2] - 4. * coeff[tptr-1] + + 6. * coeff[tptr] - 4. * coeff[tptr+1] + + coeff[tptr+2]) + cd41 = 1./120. * (coeff[tptr-1] - 4. * coeff[tptr] + + 6. * coeff[tptr+1] - 4. * coeff[tptr+2] + + coeff[tptr+3]) + + # calculate coefficients in x for each y + pcoeff[1,i] = coeff[tptr] + pcoeff[2,i] = coeff[tptr+1] - coeff[tptr] - 2. * cd20 - cd21 + + 6. * cd40 + 4. * cd41 + pcoeff[3,i] = 3. * cd20 - 5. * cd40 + pcoeff[4,i] = cd21 - cd20 - 5. * (cd40 + cd41) + pcoeff[5,i] = 5. * cd40 + pcoeff[6,i] = cd41 - cd40 + + tptr = tptr + len_coeff + } + + # calculate polynomial coefficients in y + do j = 1, 6 { + + # calculate the central differences + cd20 = 1./6. * (pcoeff[j,4] - 2. * pcoeff[j,3] + pcoeff[j,2]) + cd21 = 1./6. * (pcoeff[j,5] - 2. * pcoeff[j,4] + pcoeff[j,3]) + cd40 = 1./120. * (pcoeff[j,1] - 4. * pcoeff[j,2] + + 6. * pcoeff[j,3] - 4. * pcoeff[j,4] + pcoeff[j,5]) + cd41 = 1./120. * (pcoeff[j,2] - 4. * pcoeff[j,3] + + 6. * pcoeff[j,4] - 4. * pcoeff[j,5] + pcoeff[j,6]) + + # calculate the final coefficients + temp[1] = pcoeff[j,3] + temp[2] = pcoeff[j,4] - pcoeff[j,3] - 2. * cd20 - cd21 + + 6. * cd40 + 4. * cd41 + temp[3] = 3. * cd20 - 5. * cd40 + temp[4] = cd21 - cd20 - 5. * (cd40 + cd41) + temp[5] = 5. * cd40 + temp[6] = cd41 - cd40 + + do i = 1, 6 + pcoeff[j,i] = temp[i] + + } + +end + + +# II_PCSPLINE3 -- Procedure to evaluate the polynomial coefficients +# of bicubic spline. + +procedure ii_pcspline3 (coeff, index, len_coeff, pcoeff, len_pcoeff) + +real coeff[ARB] # 1D array of interpolant coeffcients +int index # pointer into coeff array +int len_coeff # row length of coeffcients +real pcoeff[len_pcoeff,ARB] # polynomial coefficients +int len_pcoeff # row length of pcoeff + +int tptr +int i, j +real temp[4] + +begin + # determine polynomial coefficients in x + tptr = index + do i = 1, 4 { + + pcoeff[1,i] = coeff[tptr+1] + 4. * coeff[tptr] + coeff[tptr-1] + pcoeff[2,i] = 3. * (coeff[tptr+1] - coeff[tptr-1]) + pcoeff[3,i] = 3. * (coeff[tptr-1] - 2. * coeff[tptr] + + coeff[tptr+1]) + pcoeff[4,i] = -coeff[tptr-1] + 3. * coeff[tptr] - + 3. * coeff[tptr+1] + coeff[tptr+2] + + tptr = tptr + len_coeff + } + + # calculate polynomial coefficients in y + do j = 1, 4 { + + temp[1] = pcoeff[j,3] + 4. * pcoeff[j,2] + pcoeff[j,1] + temp[2] = 3. * (pcoeff[j,3] - pcoeff[j,1]) + temp[3] = 3. * (pcoeff[j,1] - 2. * pcoeff[j,2] + pcoeff[j,3]) + temp[4] = -pcoeff[j,1] + 3. * pcoeff[j,2] - 3. * pcoeff[j,3] + + pcoeff[j,4] + + do i = 1, 4 + pcoeff[j,i] = temp[i] + } +end + +# II_BISINCDER -- Evaluate the derivatives of the 2D sinc interpolator. If the +# function value only is needed call ii_bisinc. This routine computes only +# the first 2 derivatives in x and y. The second derivative is computed +# even if only the first derivative is needed. The sinc truncation length +# is nsinc. The taper is a cosbell approximated by a quartic polynomial. +# The data value if returned if x is closer to x[i] than mindx and y is +# closer to y[i] than mindy. + +procedure ii_bisincder (x, y, der, nxder, nyder, len_der, coeff, first_point, + nxpix, nypix, nsinc, mindx, mindy) + +real x, y # the input x and y values +real der[len_der,ARB] # the output derivatives array +int nxder, nyder # the number of derivatives to compute +int len_der # the width of the derivatives array +real coeff[ARB] # the coefficient array +int first_point # offset of first data point into the array +int nxpix, nypix # size of the coefficient array +int nsinc # the sinc truncation length +real mindx, mindy # the precision of the sinc interpolant + +double sumx, normx[3], normy[3], norm[3,3], sum[3,3] +int i, j, k, jj, kk, xc, yc, nconv, index +int minj, maxj, offj, mink, maxk, offk, last_point +pointer sp, ac, ar +real sconst, a2, a4, dx, dy, dxn, dyn, dx2, taper, sdx, ax, ay, ctanx, ctany +real zx, zy +real px[3], py[3] + +begin + # Return if no derivatives ar to be computed. + if (nxder == 0 || nyder == 0) + return + + # Initialize the derivatives to zero. + do jj = 1, nyder { + do kk = 1, nxder + der[kk,jj] = 0.0 + } + + # Return if the data is outside range. + xc = nint (x) + yc = nint (y) + if (xc < 1 || xc > nxpix || yc < 1 || yc > nypix) + return + + # Call ii_bsinc if only the function value is requested. + if (nxder == 1 && nyder == 1) { + call ii_bisinc (coeff, first_point, nxpix, nypix, x, y, der[1,1], + 1, nsinc, mindx, mindy) + return + } + + # Compute the constants for the cosine bell taper approximation. + sconst = (HALFPI / nsinc) ** 2 + a2 = -0.49670 + a4 = 0.03705 + + # Allocate some working space. + nconv = 2 * nsinc + 1 + call smark (sp) + call salloc (ac, 3 * nconv, TY_REAL) + call salloc (ar, 3 * nconv, TY_REAL) + call aclrr (Memr[ac], 3 * nconv) + call aclrr (Memr[ar], 3 * nconv) + + # Initialize. + dx = x - xc + dy = y - yc + if (dx == 0.0) + ctanx = 0.0 + else + ctanx = 1.0 / tan (PI * dx) + if (dy == 0.0) + ctany = 0.0 + else + ctany = 1.0 / tan (PI * dy) + index = - 1 - nsinc + dxn = -1 - nsinc - dx + dyn = -1 - nsinc - dy + if (mod (nsinc, 2) == 0) + sdx = 1.0 + else + sdx = -1.0 + do jj = 1, 3 { + normy[jj] = 0.0d0 + normx[jj] = 0.0d0 + } + + do i = 1, nconv { + dx2 = sconst * (i + index) ** 2 + taper = sdx * (1.0 + a2 * dx2 + a4 * dx2 * dx2) ** 2 + #ax = dxn + i + #ay = dyn + i + ax = -dxn - i + ay = -dyn - i + if (ax == 0.0) { + px[1] = 1.0 + px[2] = 0.0 + px[3] = - 1.0 / 3.0 + } else if (dx == 0.0) { + px[1] = 0.0 + px[2] = 0.0 + px[3] = 0.0 + } else { + zx = 1.0 / ax + px[1] = taper * zx + px[2] = px[1] * (ctanx - zx) + px[3] = -px[1] * (1.0 + 2.0 * zx * (ctanx - zx)) + } + if (ay == 0.0) { + py[1] = 1.0 + py[2] = 0.0 + py[3] = - 1.0 / 3.0 + } else if (dy == 0.0) { + py[1] = 0.0 + py[2] = 0.0 + py[3] = 0.0 + } else { + zy = 1.0 / ay + py[1] = taper * zy + py[2] = py[1] * (ctany - zy) + py[3] = -py[1] * (1.0 + 2.0 * zy * (ctany - zy)) + } + + Memr[ac+i-1] = px[1] + Memr[ac+nconv+i-1] = px[2] + Memr[ac+2*nconv+i-1] = px[3] + Memr[ar+i-1] = py[1] + Memr[ar+nconv+i-1] = py[2] + Memr[ar+2*nconv+i-1] = py[3] + + do jj = 1, 3 { + normx[jj] = normx[jj] + px[jj] + normy[jj] = normy[jj] + py[jj] + } + + sdx = -sdx + } + + # Normalize. + do jj = 1, 3 { + do kk = 1, 3 + norm[kk,jj] = normx[kk] * normy[jj] + } + + + # Do the convolution. + minj = max (1, yc - nsinc) + maxj = min (nypix, yc + nsinc) + mink = max (1, xc - nsinc) + maxk = min (nxpix, xc + nsinc) + do jj = 1, nyder { + offj = ar + (jj - 1) * nconv - yc + nsinc + do kk = 1, nxder { + offk = ac + (kk - 1) * nconv - xc + nsinc + sum[kk,jj] = 0.0d0 + + # Do the convolutions. + do j = yc - nsinc, minj - 1 { + sumx = 0.0d0 + do k = xc - nsinc, mink - 1 + sumx = sumx + Memr[k+offk] * coeff[first_point+1] + do k = mink, maxk + sumx = sumx + Memr[k+offk] * coeff[first_point+k] + do k = maxk + 1, xc + nsinc + sumx = sumx + Memr[k+offk] * coeff[first_point+nxpix] + + sum[kk,jj] = sum[kk,jj] + Memr[j+offj] * sumx + } + + + do j = minj, maxj { + index = first_point + (j - 1) * nxpix + sumx = 0.0d0 + do k = xc - nsinc, mink - 1 + sumx = sumx + Memr[k+offk] * coeff[index+1] + do k = mink, maxk + sumx = sumx + Memr[k+offk] * coeff[index+k] + do k = maxk + 1, xc + nsinc + sumx = sumx + Memr[k+offk] * coeff[index+nxpix] + + sum[kk,jj] = sum[kk,jj] + Memr[j+offj] * sumx + } + + do j = maxj + 1, yc + nsinc { + last_point = first_point + (nypix - 1) * nxpix + sumx = 0.0d0 + do k = xc - nsinc, mink - 1 + sumx = sumx + Memr[k+offk] * coeff[last_point+1] + do k = mink, maxk + sumx = sumx + Memr[k+offk] * coeff[last_point+k] + do k = maxk + 1, xc + nsinc + sumx = sumx + Memr[k+offk] * coeff[last_point+nxpix] + + sum[kk,jj] = sum[kk,jj] + Memr[j+offj] * sumx + } + + } + } + + # Build the derivatives. + der[1,1] = sum[1,1] / norm[1,1] + if (nxder > 1) + der[2,1] = sum[2,1] / norm[1,1] - (sum[1,1] * norm[2,1]) / + norm[1,1] ** 2 + if (nxder > 2) + der[3,1] = sum[3,1] / norm[1,1] - (norm[3,1] * sum[1,1] + + 2.0d0 * sum[2,1] * norm[2,1]) / norm[1,1] ** 2 + + 2.0d0 * sum[1,1] * norm[2,1] * norm[2,1] / norm[1,1] ** 3 + if (nyder > 1) { + der[1,2] = sum[1,2] / norm[1,1] - (sum[1,1] * norm[1,2]) / + norm[1,1] ** 2 + if (nxder > 1) + der[2,2] = sum[2,2] / norm[1,1] - (sum[2,1] * norm[1,2] + + sum[1,2] * norm[2,1] + norm[2,2] * sum[1,1]) / + norm[1,1] ** 2 + (2.0d0 * sum[1,1] * norm[2,1] * + norm[1,2]) / norm[1,1] ** 3 + if (nxder > 2) + der[3,2] = sum[3,2] / norm[1,1] - (sum[3,1] * norm[1,2] + + 2.0 * norm[2,2] * sum[2,1] + 2.0 * sum[2,2] * + norm[2,1] + norm[3,1] * sum[1,2] + norm[3,2] * + sum[1,1]) / norm[1,1] ** 2 + (4.0 * norm[2,1] * + sum[2,1] * norm[1,2] + 2.0 * norm[2,1] * sum[1,2] * + norm[2,1] + 4.0 * norm[2,1] * norm[2,2] * sum[1,1] + + 2.0 * norm[3,1] * norm[1,2] * sum[1,1]) / + norm[1,1] ** 3 - 6.0 * norm[2,1] * norm[2,1] * + norm[1,2] * sum[1,1] / norm[1,1] ** 4 + + } + if (nyder > 2) { + der[1,3] = sum[1,3] / norm[1,1] - (norm[1,3] * sum[1,1] + + 2.0d0 * sum[1,2] * norm[1,2]) / norm[1,1] ** 2 + + 2.0d0 * sum[1,1] * norm[1,2] * norm[1,2] / norm[1,1] ** 3 + if (nxder > 1) + der[2,3] = sum[2,3] / norm[1,1] - (sum[1,3] * norm[2,1] + + 2.0 * norm[2,2] * sum[1,2] + 2.0 * sum[2,2] * + norm[1,2] + norm[1,3] * sum[2,1] + norm[2,3] * + sum[1,1]) / norm[1,1] ** 2 + (4.0 * norm[1,2] * + sum[1,2] * norm[2,1] + 2.0 * norm[1,2] * sum[2,1] * + norm[1,2] + 4.0 * norm[1,2] * norm[2,2] * sum[1,1] + + 2.0 * norm[1,3] * norm[2,1] * sum[1,1]) / + norm[1,1] ** 3 - 6.0 * norm[1,2] * norm[1,2] * + norm[2,1] * sum[1,1] / norm[1,1] ** 4 + if (nxder > 2) + der[3,3] = sum[3,3] / norm[1,1] - (2.0 * sum[2,3] * norm[2,1] + + norm[3,1] * sum[1,3] + 2.0 * norm[3,2] * sum[1,2] + + 4.0 * sum[2,2] * norm[2,2] + 2.0 * sum[3,2] * + norm[1,2] + 2.0 * norm[2,3] * sum[2,1] + sum[3,1] * + norm[1,3] + norm[3,3] * sum[1,1]) / norm[1,1] ** 2 + + (2.0 * norm[2,1] * norm[2,1] * sum[1,3] + 8.0 * + norm[2,1] * sum[1,2] * norm[2,2] + 8.0 * norm[2,1] * + norm[1,2] * sum[2,2] + 4.0 * norm[2,1] * sum[2,1] * + norm[1,3] + 4.0 * norm[2,1] * norm[2,3] * sum[1,1] + + 4.0 * norm[1,2] * sum[1,2] * norm[3,1] + 8.0 * + norm[2,2] * sum[2,1] * norm[1,2] + 2.0 * norm[1,2] * + norm[1,2] * sum[3,1] + 4.0 * norm[2,2] * norm[2,2] * + sum[1,1] + 4.0 * norm[1,2] * norm[3,2] * sum[1,1] + + 2.0 * norm[1,3] * norm[3,1] * sum[1,1]) / + norm[1,1] ** 3 - (12.0 * norm[2,1] * norm[2,1] * + norm[1,2] * sum[1,2] + 12.0 * norm[2,1] * norm[1,2] * + norm[1,2] * sum[2,1] + 24.0 * norm[2,1] * norm[1,2] * + norm[2,2] * sum[1,1] + 6.0 * norm[2,1] * norm[2,1] * + norm[1,3] * sum[1,1] + 6.0 * norm[1,2] * norm[1,2] * + norm[3,1] * sum[1,1]) / norm[1,1] ** 4 + ( 24.0 * + norm[1,2] * norm[1,2] * norm[2,1] * norm[2,1] * + sum[1,1]) / norm[1,1] ** 5 + } + + + + + call sfree (sp) +end |