Repatch (from linux) of OSX IRAF

1 files changed, 525 insertions, 0 deletions

diff --git a/math/iminterp/doc/im1dinterp.spc b/math/iminterp/doc/im1dinterp.spc
new file mode 100644
index 00000000..ce5b8680
--- /dev/null
+++ b/math/iminterp/doc/im1dinterp.spc
@@ -0,0 +1,525 @@
+.help iminterp Jul84 "Math Package"
+
+.ce
+Specifications for the Image Interpolator Package
+.ce
+Lindsey Davis
+.ce
+Vesa Junkkarinen
+.ce
+August 1984
+
+.sh
+1. Introduction
+
+    One of the most common operations in image processing is
+interpolation in a data array. Due to the large amount of data involved,
+efficiency is highly important. The advantage of having locally written
+interpolators, includes the ability to optimize for uniformly spaced data
+and the possibility of adding features that are useful to the final
+application.
+
+.sh
+2. Requirements
+
+.ls (1)
+The package shall take as input a one-dimensional array containing image
+data. The pixels are assumed to be equally spaced along a line.
+The coordinates of a pixel are assumed to be
+the same as the subscript of the pixel in the data array.
+The coordinate of the first pixel in the array and the spacing between pixels
+is assumed to be 1.0. All pixels are assumed to be good.
+Checking for INDEF valued and out of bounds pixels is the responsibility of the
+user. A routine to remove INDEF valued pixels from a data array shall be
+included in the package.
+.le
+.ls (2)
+The package is divided into array sequential interpolators and array
+random interpolators. The sequential interpolators have been optimized
+for returning many values as is the case when an array is shifted, or
+oversampled at many points in order to produce a
+smooth plot.
+The random interpolators allow the evaluation of a few interpolated
+points without the computing time and storage overhead required for
+setting up the sequential version.
+.le
+.ls (3)
+The quality of the interpolant will be set at run time. The options are:
+
+.nf
+    II_NEAREST		- nearest neighbour
+    II_LINEAR		- linear interpolation
+    II_POLY3		- 3rd order divided differences
+    II_POLY5		- 5th order divided differences
+    II_SPLINE3		- cubic spline
+.fi
+
+The calling sequences shall be invariant to the interpolant selected.
+Routines should be designed so that new interpolants can be added
+with minimal changes to the code and no change to the calling sequences.
+.le
+.ls (4)
+The interpolant parameters and the arrays necessary to store the coefficients
+are stored in a structure referenced by a pointer. The pointer is returned
+to the user program by the initial call to ASIINIT or ASIRESTORE and freed
+by a call to ASIFREE (see section 3.1).
+.le
+.ls (5)
+The package routines shall be able to:
+.ls o
+Calculate the coefficients of the interpolant and store these coefficients in
+the appropriate part of the interpolant descriptor structure.
+.le
+.ls o
+Evaluate the interplant at a given x(s) coordinate(s).
+.le
+.ls o
+Calculate the derivatives of the interpolant at a given value of x.
+.le
+.ls o
+Integrate the interpolant over a specified x interval.
+.le
+.ls o
+Store the interpolant in a user supplied array. Restore the saved interpolant
+to the interpolant descriptor structure for later use by ASIEVAL, ASIVECTOR,
+ASIDER and ASIGRL.
+.le
+
+.sh
+3. Specifications
+
+.sh
+3.1. The Array Sequential Interpolator Routines
+
+    The package prefix is asi and the package routines are:
+
+.nf
+        asiinit	(asi, interp_type)
+         asifit	(asi, datain, npix)
+    y = asieval	(asi, x)
+      asivector	(asi, x, yfit, npix)
+         asider	(asi, x, der, nder)
+     v = asigrl	(asi, a, b)
+        asisave	(asi, interpolant)
+     asirestore	(asi, interpolant)
+        asifree	(asi)
+.fi
+
+.sh
+3.2. The Array Random Interpolator Routines
+
+    The package prefix is ari and the package routines are:
+
+.nf
+    y = arieval	(x, datain, npix, interp_type)
+	 arider	(x, datain, npix, der, nder, interp_type)
+.fi
+
+.sh
+3.3. Miscellaneous
+
+    A routine has been included in the package to remove INDEF valued
+pixels from an array.
+
+.nf
+    arbpix (datain, dataout, npix, interp_type, boundary_type)
+.fi
+
+.sh
+3.4. Algorithms
+
+.sh
+3.4.1. Coefficients
+
+    The coefficient array used by the asi routines is calculated by ASIFIT.
+ASIFIT accepts an array of data, checks that the number
+of data points is appropriate for the interpolant selected, allocates
+space for the interpolant, and calculates the coefficients.
+Boundary coefficient values are calculated
+using boundary projection.  With the exception of the cubic spline interpolant,
+the coefficients are stored as the data points.
+The B-spline coefficients are 
+calculated using natural end conditions (Prenter 1975).
+After a call to ASIFIT the coefficient array contains the following.
+
+.nf
+    case II_NEAREST:
+
+        # no boundary conditions necessary
+	coeff[1] = datain[1]
+	.
+	.
+	.
+	coeff[npts] = datain[npix]
+
+    case II_LINEAR:
+
+        # coeff[npxix+1] required if x = npix
+	coeff[1] = datain[1]
+	.
+	.
+	.
+	coeff[npix] = datain[npix]
+	coeff[npix+1] = 2. * datain[npix] - datain[npix-1]
+
+    case II_POLY3:
+
+        # coeff[0] required if x = 1
+	# coeff[npix+1], coeff[npix+2] required if x = npix
+	coeff[0] = 2. * datain[1] - datain[2]
+	coeff[1] = datain[1]
+	.
+	.
+	.
+	coeff[npix] = datain[npix]
+	coeff[npix+1] = 2. * datain[npix] - datain[npix-1]
+	coeff[npix+2] = 2. * datain[npix] - datain[npix-2]
+
+    case II_POLY5:
+
+        # coeff[1], coeff[0] reqired if x = 1
+	# coeff[npix+1], coeff[npix+2], coeff[npix=3]
+	# required if x = npix
+
+	coeff[-1] = 2. * datain[1] - datain[3]
+	coeff[0] = 2. * datain[1] - datain[2]
+	coeff[1] = datain[1]
+	.
+	.
+	.
+	coeff[npix] = datain[npix]
+	coeff[npix+1] = 2. * datain[npix] - datain[npix-1]
+	coeff[npix+2] = 2. * datain[npix] - datain[npix-2]
+	coeff[npix+3] = 2. * datain[npix] - datain[npix-3]
+
+    case SPLINE3:
+
+        # coeff[0] = 2nd der at x = 1, coeff[0] = 0.
+	# coeff[npix+1] = 2nd der at x = npts, coeff[npix+1] = 0.
+	# coeff[npix+2] = 0., required if x = npix
+	coeff[0] = b[1]
+	coeff[1] = b[2]
+	.
+	.
+	.
+	coeff[npix] = b[npix+1]
+	coeff[npix+1] = b[npix+2]
+	coeff[npix+2] = 0.
+.fi
+
+.sh
+3.4.2. Evaluation
+
+    The ASIEVAL and ASIVECTOR routines have been optimized to be as efficient
+as possible. The values of the II_NEAREST and II_LINEAR interpolants
+are calculated directly. The II_SPLINE3 interpolant is evaluated using
+polynomial coefficients calculated directly from the B-spline coefficients
+(de Boor 1978). Values of the higher order polynomial interpolants
+are calculated using central differences. The equations for each case are
+listed below.
+
+.nf
+case II_NEAREST:
+    
+    y = coeff[int (x + 0.5)]
+
+case II_LINEAR:
+
+    nx = x
+    y = (x - nx) * coeff[nx+1] + (nx + 1 - x) * coeff[nx]
+
+case II_POLY3:
+
+    nx = x
+    s = x - nx
+    t = 1. - s
+
+    # second central differences
+    cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1])
+    cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx])
+
+    y = s * (coeff[nx+1] + (s * s - 1.) * cd21) + t * (coeff[nx] +
+	(t * t - 1.) * cd20)
+
+case II_POLY5:
+
+    nx = x
+    s = x - nx
+    t = 1. - s
+
+    # second central differences
+    cd20 = 1./6. * (coeff[nx+1] - 2. * coeff[nx] + coeff[nx-1])
+    cd21 = 1./6. * (coeff[nx+2] - 2. * coeff[nx+1] + coeff[nx])
+
+    # fourth central diffreences
+    cd40 = 1./120. * (coeff[nx-2] - 4. * coeff[nx-1] + 6. * coeff[nx] - 4. *
+	   coeff[nx+1] + a[nx+2])
+    cd41 = 1./120. * (coeff[nx-1] - 4. * coeff[nx] + 6. * coeff[nx+1] - 4. *
+	   coeff[nx+2] + coeff[nx+3]
+
+    y = s * (coeff[nx+1] + (s * s - 1.) * (cd21 + (s * s - 4.) * cd41)) +
+	t * (coeff[nx] + (t * t - 1.) * (cd20 + (t * t - 4.) * cd40))
+
+case II_SPLINE3:
+
+    nx = x
+    s = x - nx
+
+    pc[1] = coeff[nx-1] + 4. * coeff[nx] + coeff[nx+1]
+    pc[2] = 3. * (coeff[nx+1] - coeff[nx-1])
+    pc[3] = 3. * (coeff[nx-1] - 2. * coeff[nx] + coeff[nx+1])
+    pc[4] = -coeff[nx-1] + 3. * coeff[nx] - 3. * coeff[nx+1] + coeff[nx+2]
+
+    y = pc[1] + s * (pc[2] + s * (pc[3] + s * pc[4]))
+.fi
+
+
+    The ARIEVAL routine uses the expressions above to evaluate the
+interpolant. However unlike ASIEVAL, ARIEVAL does not use a previously
+calculated coefficient array. Instead ARIEVAL selects the appropriate
+portion of the data array, calculates the coefficients and boundary
+coefficients if necessary, and evaluates the interpolant at the time it
+is called. The cubic spline interpolant uses at most SPLTS (currently 16)
+data points to calculate the B-spline coefficients.
+
+.sh
+3.4.3. Derivatives
+
+    Derivatives of the interpolant are calculated by evaluating the
+derivatives of the interpolating polynomial. For all interpolants der[1]
+equals the value of the interpolant at x.
+For the sake of efficiency the derivatives
+of the II_NEAREST and II_LINEAR interpolants are calculated directly.
+
+.nf
+    case II_NEAREST:
+
+	der[1] = coeff[int (x+0.5)]
+
+    case II_LINEAR:
+
+	der[1] = (x - nx) * coeff [nx+1] + (nx + 1 - x) * coeff[nx]
+	der[2] = coeff[nx+1] - coeff[nx]
+.fi
+
+    In order to calculate the derivatives of the cubic spline and
+polynomial interpolants
+the coefficients of the interpolating polynomial must be calculated.
+The polynomial
+coefficients for the cubic spline interpolant are computed directly from the
+B-spline coefficients (see 3.4.2.). The higher order polynomial
+interpolant coefficients are calculated as follows.
+
+First the appropriate portion of the coefficient array is loaded.
+
+.nf
+	do i = 1, nterms
+	    d[i] = coeff[nx - nterms/2 + i]
+.fi
+
+Next the divided differences are calculated (Conte and de Boor 1972).
+
+.nf
+	do k = 1, nterms - 1
+	    do i = 1, nterms - k
+		d[i] = (d[i+1] - d[i]) / k
+.fi
+
+The d[i] are the coefficients of an interpolating polynomial of the
+following form. The x[i] are the nterms data points surrounding the
+point of interest.
+
+.nf
+	p(x) = d[1] * (x-x[1]) * ... * (x-x[nterms-1) +
+	       d[2] * (x-x[2]) * ... * (x-x[nterms-1]) + ... + d[nterms]
+.fi
+
+Next a new set of polynomial coefficients are calculated
+(Conte and de Boor 1972).
+
+.nf
+	do k = nterms, 2, -1
+	    do i = 2, k
+		d[i] = d[i] + d[i-1] * (k - i - nterms/2)
+.fi
+
+The new d[i] are the coefficients of the follwoing polynomial.
+
+.nf
+	nx = x
+	p(x) = d[1] * (x-nx) ** (nterms-1) + d[2] * (x-nx) ** (nterms-2) + ...
+	       d[nterms]
+.fi
+
+The d[i] array is flipped. The value and derivatives
+of the interpolant are then calculated using the d[i] array and
+nested multiplication.
+
+.nf
+	s = x - nx
+
+	do k = 1, nder {
+
+	    accum = d[nterms-k+1]
+
+	    do j = nterms - k, 1, -1
+		accum = d[j] + s * accum
+
+	    der[k] = accum
+
+	    # differnetiate
+	    do j = 1, nterms - k
+		d[j] = j * d[j + 1]
+	}
+.fi
+
+    ARIDER calculates the derivatives of the interpolant using the same
+technique ASIDER. However ARIDER does not use a previously calculated
+coefficient array like ASIDER. Instead ARIDER selects the appropriate portion
+of the data array, calculates the coefficients and boundary coefficients,
+and computes the derivatives at the time it is called.
+
+.sh
+3.4.5. Integration
+
+    ASIGRL calculates the integral of the interpolant between fixed limits
+by integrating the interpolating polynomial. The coefficients of the
+interpolating polynomial are calculated as discussed in section 3.4.4.
+
+.sh 
+4. Usage
+
+.sh
+4.1. User Notes
+
+The following series of steps illustrates the use of the package.
+
+.ls 4
+.ls (1)
+Insert an include <iminterp.h> statement in the calling program to make
+the IINTERP definitions available to the user program.
+.le
+.ls (2)
+Remove INDEF valued pixels from the data using ARBPIX.
+.le
+.ls (3)
+Call ASIINIT to initialize the interpolant parameters.
+.le
+.ls (4)
+Call ASIFIT to calculate the coefficients of the interpolant.
+.le
+.ls (5)
+Evaluate the interpolant at a given value of x(s) using ASIEVAL or
+ASIVECTOR.
+.le
+.ls (6)
+Calculate the derivatives and integral or the interpolant using
+ASIDER and ASIGRL.
+.le
+.ls (7)
+Free the interpolator structure by calling ASIFREE.
+.le
+.le
+
+    The interpolant can be saved and restored using ASISAVE and ASIRESTORE.
+If the values and derivatives of only a few points in an array are desired
+ARIEVAL and ARIDER can be called.
+
+.sh
+4.2. Examples
+
+.nf
+Example 1: Shift a data array by a constant amount
+
+    include <iminterp.h>
+    ...
+    call asiinit (asi, II_POLY5)
+    call asifit (asi, inrow, npix)
+
+    do i = 1, npix
+	outrow[i] = asieval (asi, i + shift)
+
+    call asifree (asi)
+    ...
+
+Example 2: Calculate the integral under the data array
+
+    include <iminterp.h>
+    ...
+    call asiinit (asi, II_POLY5)
+    call asifit (asi, datain, npix)
+
+    integral =  asigrl (asi, 1. real (npix))
+
+    call asifree (asi)
+    ...
+
+Example 2: Store interpolant for later use by ASIEVAL
+	   LEN_INTERP must be at least npix + 8 units long where npix is
+	   is defined in the call to ASIFIT.
+
+    include <iminterp.h>
+
+    real   interpolant[LEN_INTERP]
+    ...
+    call asiinit (asi, II_POLY3)
+    call asifit (asi, datain, npix)
+    call asisave (asi, interpolant)
+    call asifree (asi)
+    ...
+    call asirestore (asi, interpolant)
+    do i = 1, npts
+	yfit[i] = asieval (asi, x[i])
+    call asifree (asi)
+    ...
+.fi
+.sh
+5. Detailed Design
+
+.sh
+5.1. Interpolator Descriptor
+
+    The interpolant parameters and coefficients are stored in a
+structure listed below.
+
+.nf
+    define	LEN_ASISTRUCT	4		# Length in structure units of
+						# interpolant descriptor
+
+    define	ASI_TYPE	Memi[$1]	# Interpolant type
+    define	ASI_NCOEFF	Memi[$1+1]	# No. of coefficients
+    define	ASI_OFFSET	Memi[$1+2]	# First "data" point in
+						# coefficient array
+    define	ASI_COEFF	Memi[$1+3]	# Pointer to coefficient array
+.fi
+
+.sh
+5.2. Storage Requirements
+
+    The interpolant descriptor requires LEN_ASISTRUCT storage units. The
+coefficient array requires ASI_NCOEFF(asi) real storage elements, where
+ASI_NCOEFF(asi) is defined as follows.
+
+.nf
+    ASI_NCOEFF(asi) = npix	# II_NEAREST
+    ASI_NCOEFF(asi) = npix+1	# II_LINEAR
+    ASI_NCOEFF(asi) = npix+3	# II_POLY3
+    ASI_NCOEFF(asi) = npix+5	# II_POLY5
+    ASI_NCOEFF(asi) = npix+3	# II_SPLINE3
+.fi
+
+.sh
+6. References
+
+.ls (1)
+Carl de Boor, "A Practical Guide to Splines", 1978, Springer-Verlag New
+York Inc.
+.le
+.ls (2)
+S.D. Conte and C. de Boor, "Elementary Numerical Analysis", 1972, McGraw-Hill,
+Inc.
+.le
+.ls (3)
+P.M. Prenter, "Splines and Variational Methods", 1975, John Wiley and Sons Inc.
+.le
+.endhelp