Notes on the Interpolator Package for IRAF Version 2

	Some parts of the image interpolator package are available.
This is a reminder of changes to the document TRP prepared by
Doug Tody.  Also it is a place to keep track of why things were
done the way they were.  The individual procedures (subroutines)
will be discussed; making it easy to add, subtract and change things.

**************************************************************************
**									**
**									**

		GENERAL NOTES FOR 1-D PACKAGE:

1. Handling of errors.

	Some routines call the error routines of the IRAF package.
	Logically it is reasonable to trap silly errors at certain points.
	For instance if the interpolator type is set to a nonexistant type
	an error call is generated.  

	To the programmer, the two most important error features are:

	1. All pixels are assumed to be good.

	2. All x references are assumed to land in the closed interval
	   1 <= x <= NPTS.
	
	Point 2 probably needs some explanation, because at first hand
	it seems rather easy to test for x < 1 || x > NPTS and take
	appropriate action.  There are two objections to including this
	test.  First it is often possible to generate x values in such
	a way that they are guaranteed to lie in the given interval.
	Testing again within the interpolator package is a waste of
	computer time.  Second the interpolator package becomes too large
	if code is included to handle out of bounds references and such
	code duplicates other parts of the IRAF package that handles
	out of bound values.  The alternative is to return INDEF when
	x is out of bounds.  This is not very appealing for many
	routine applications such as picture shifting.

2. Size of "coeff" array:

	coeff is a real 1-d array of size  -- 2 * NPTS + SZ_ASI

    where SZ_ASI is currently 30.  It, like the other defines needed for
    the interp package, is available in the file -- terpdef.h .

       This coeff array serves as work area for asifit, a structure to save
    interpolator type, flags etc., an array to save data points or b-spline
    coefficients as needed.

   The polynomial interpolators are saved as the data points and the
interpolation is done "on the fly".  Rather than compute polynomial
coefficients before-hand, the interpolated value is calculated from
the data values at time of request.  (The actual calculation is done
by adding and multiplying things in an order suggested by Everett's
interpolation formula rather than the polynomial coefficient form.
The resulting values are the same up to differences due to numerical
roundoff.  As a practical matter, all forms work except polynomials
shifted far from the point of interest.  Polynomials shifted to a
center far away suffer from bad roundoff errors.)

   The splines must be calculated at a single pass because (formally)
they are global (i.e use all of the data points).  The b-splines (b for
basis) are just certain piecewise polynomials for solving the interpo-
lation problem, it is convenient to calculate the coefficients of these
b-splines and to save these in the array.

3. Method of solution:

**************************************************************************
**									**
**									**
**									**
			PROCEDURES:

***************************************************************************

	*****        asiset(coeff,interptype) 		       *****

The interpolator type can be set to:

	IT_NEAREST		nearest neighbor
	IT_LINEAR		linear interpolation
	IT_POLY3		interior polynomial 3rd order
	IT_POLY5		interior polynomial 5th order
	IT_SPLINE3		cubic natural spline

Subroutines needed:

	NONE



*************************************************************************

	*******		asifit (data, npts, coeff)            *****

   This fits the interpolator to the data.  Takes care of bad points by
replacing them with interpolated values to generate a
uniformly spaced array.  For polynomial interpolators, this array is
sent to the evaluation routines.  For the spline interpolator, the uniform
array is fitted and the basis-spline coefficients are saved.

Interior polynomial near boundary:
   Array is extended according to choice of boundary conditions.
Out of bounds values and values near the edge depend on choice of
boundary conditions.

Spline near boundary:
   Natural (second derivative equals zero at edge) spline is fitted to
data only.  The resulting function is reflected, wrapped etc. to generate
out of bounds values.  Out of bounds values depend on choice of boundary
conditions but values within the array near the edges do not.
   For the wrap boundary condition, the one segment connecting point
number n to point 1 is not generated by the spline fit.  Instead the
polynomial that is continuous and has a continuous first and second
derivative n is inserted.  In this one case, the first and second derivative
are not continuous at 1.

Replacing bad points:
   The same kind of interpolator is used to replace bad points as that
used to generate interpolated values.  The boundary conditions are not
handled in the same manner.  Bad points near the edge are replaced by
using lower order polynomial interpolators if there are not enough
points on both sides of the bad point.  Bad points that are not
straddled by good points are assigned the value of the nearest good
point. This is done in the temporary array before the uniformly spaced
interpolator is applied.

subroutines needed:

	real iipint(x,y,nterm,x0)
		Returns interpolated value at x0 using polynomial with
			<nterm> terms to connect the <nterm> data pairs
			in x,y.

	cubspl(tau,c,n,ibcbeg,ibcend)
		Unmodified cubic spline interpolator for non-uniformly
			spaced data taken from C. de Boor's book:
			A Practical Guide to Splines, (1978
			Springer-Verlag).

	iifits(b,d,n)
		Fits cubic spline to uniformly spaced array of y values
			using a basis-spline representation for
			the spline.



*****************************************************************************

	*****		y = asival (x, coeff)                  *****

   Subroutine to return interpolated value after fitting.  

Subroutines needed:

	real iivalu (x, coeff)
		Returns interpolated value assuming x is in array.  Also
			has no code for bad point propagation.




***************************************************************************

	*****		asieva (x, y, n, coeff)                  *****

   Given an ordered array of x values returns interpolated values in y.
This is needed for speed.  Reduces number of subroutine calls per point,
reduces the number of tests for out of bounds, and bad pixel
propagation can be made more efficient.

Subroutines needed:

	real iivalu (x, coeff)



**************************************************************************

	*****		asider (x, derivs, nderiv, coeff)      *****

   This evaluates derivatives at point x.  Returns the first
nderiv - 1 derivatives.  derivs[1] is function value, derivs[2] is
first derivative, derivs[3] is second derivative etc.
   There is no check that nderiv is a reasonable value.  Polynomials
have a limited number of non-zero derivatives, if more are asked for
they come back as zero.

Subroutines needed:

	iivalu (x, coeff)

	iiderv (x, derivs, nderiv, coeff)
		Assumes x lands in array returns derivatives.  No handling
			of bad points.

****************************************************************************

	******		integral = asigrl (a, b, coeff)           *******

    This integrates the array from to a to b.  The boundary conditions are
incorporated into the code so if function values can be returned from all
points between a and b, the integral will be defined.  Thus for WRAP,
REFLECT and PROJECT boundary conditions, the distance from a to b can be
as large as 3 times the array length if a and b are chosen appropriately.
For these boundary conditions, attempts to extend beyond +/- one cycle
produce INDEF.

Subroutines needed:

	iiigrl (a, b, coeff)
		returns integral when a and b land in array.