From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001
From: Joseph Hunkeler <jhunkeler@gmail.com>
Date: Wed, 8 Jul 2015 20:46:52 -0400
Subject: Initial commit

---
 math/iminterp/asifit.x | 146 +++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 146 insertions(+)
 create mode 100644 math/iminterp/asifit.x

(limited to 'math/iminterp/asifit.x')

diff --git a/math/iminterp/asifit.x b/math/iminterp/asifit.x
new file mode 100644
index 00000000..3ea33041
--- /dev/null
+++ b/math/iminterp/asifit.x
@@ -0,0 +1,146 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/iminterp.h>
+include "im1interpdef.h"
+
+define	TEMP	Memr[P2P($1)]
+
+# ASIFIT -- Fit the interpolant to the data.
+
+procedure asifit (asi, datain, npix)
+
+pointer	asi		# interpolant descriptor
+real	datain[ARB]	# data array
+int	npix		# nunber of data points
+
+int	i
+pointer	c0ptr, cdataptr, cnptr, temp
+
+begin
+	# Check the data array for size and allocate space for the coefficient
+	# array.
+
+	switch (ASI_TYPE(asi)) {
+
+	case II_SPLINE3:
+	    if (npix < 4)
+		call error (0, "ASIFIT: too few points for SPLINE3")
+	    else {
+		ASI_NCOEFF(asi) = npix + 3
+		ASI_OFFSET(asi) = 1
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+		call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+		call malloc (temp, ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	case II_POLY5:
+	    if (npix < 6)
+		call error (0,"ASIFIT: too few points for POLY5")
+	    else {
+		ASI_NCOEFF(asi) = npix + 5
+		ASI_OFFSET(asi) = 2
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+		call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	case II_POLY3:
+	    if (npix < 4)
+		call error (0, "ASIFIT: too few points for POLY3")
+	    else {
+		ASI_NCOEFF(asi) = npix + 3
+		ASI_OFFSET(asi) = 1
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+		call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	case II_DRIZZLE, II_LINEAR:
+	    if (npix < 2)
+		call error (0, "ASIFIT: too few points for LINEAR")
+	    else {
+		ASI_NCOEFF(asi) = npix + 1
+		ASI_OFFSET(asi) = 0
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+	        call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	case II_SINC, II_LSINC:
+	    if (npix < 1)
+		call error (0, "ASIFIT: too few points for SINC")
+	    else {
+		ASI_NCOEFF(asi) = npix
+		ASI_OFFSET(asi) = 0
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+	        call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	default:
+	    if (npix < 1)
+		call error (0," ASIFIT: too few points for NEAREST")
+	    else {
+		ASI_NCOEFF(asi) = npix
+		ASI_OFFSET(asi) = 0
+		if (ASI_COEFF(asi) != NULL)
+		    call mfree (ASI_COEFF(asi), TY_REAL)
+		call malloc (ASI_COEFF(asi), ASI_NCOEFF(asi), TY_REAL)
+	    }
+
+	}
+
+
+	# Define the pointers.
+	# (c0ptr + 1) points to first element in the coefficient array.
+	# (cdataptr + 1) points to first data element in the coefficient array.
+	# (cnptr + 1) points to the first element after the last data point in
+	# coefficient array.
+
+	c0ptr = ASI_COEFF(asi) - 1
+	cdataptr = ASI_COEFF(asi) - 1 + ASI_OFFSET(asi)
+	cnptr = cdataptr + npix
+
+	# Put data into the interpolant structure.
+	do i = 1, npix 
+	    COEFF(cdataptr + i) = datain[i]
+
+	# Specify the end conditions.
+	switch (ASI_TYPE(asi)) {
+
+	case II_SPLINE3:
+	    # Natural end conditions - second deriv. zero
+	    COEFF(c0ptr + 1) = 0.
+	    COEFF(cnptr + 1) = 0.
+	    COEFF(cnptr + 2) = 0.	# if x = npts
+
+	    # Fit spline - generate b-spline coefficients.
+	    call ii_spline (COEFF(ASI_COEFF(asi)), TEMP(temp), npix)
+	    call mfree (temp, TY_REAL)
+
+	case II_NEAREST, II_SINC, II_LSINC:
+	    # No end conditions required.
+
+	case II_LINEAR, II_DRIZZLE:
+	    COEFF(cnptr + 1) = 2. * COEFF(cdataptr + npix) -	# if x = npts
+	    		       COEFF(cdataptr + npix - 1)
+
+	case II_POLY3:
+	    COEFF(c0ptr + 1) = 2. * COEFF(cdataptr + 1) - COEFF(cdataptr + 2) 
+	    COEFF(cnptr + 1) = 2. * COEFF(cdataptr + npix) -
+	    		       COEFF(cdataptr + npix - 1)
+	    COEFF(cnptr + 2) = 2. * COEFF(cdataptr + npix) -
+	    		       COEFF(cdataptr + npix - 2) 
+
+	case II_POLY5:
+	    COEFF(c0ptr + 1) = 2. * COEFF(cdataptr + 1) - COEFF(cdataptr + 3)
+	    COEFF(c0ptr + 2) = 2. * COEFF(cdataptr + 1) - COEFF(cdataptr + 2)
+	    COEFF(cnptr + 1) = 2. * COEFF(cdataptr + npix) -
+	    		       COEFF(cdataptr + npix - 1)
+	    COEFF(cnptr + 2) = 2. * COEFF(cdataptr + npix) -
+	    		       COEFF(cdataptr + npix - 2)
+	    COEFF(cnptr + 3) = 2. * COEFF(cdataptr + npix) -
+	    		       COEFF(cdataptr + npix - 3)
+	}
+end
-- 
cgit