I'm very sorry, I previously sent the wrong patch. Here's the correct patch. (The other was a previous version of arrayfnsmodule.c, that added an interpf() command. This version adds an optional typecode to allow specifying single precision results from interp) -- Joe VanAndel National Center for Atmospheric Research http://www.atd.ucar.edu/~vanandel/ Internet: vanandel@ucar.edu Index: arrayfnsmodule.c =================================================================== RCS file: /cvsroot/numpy/Numerical/Src/arrayfnsmodule.c,v retrieving revision 1.1 diff -u -r1.1 arrayfnsmodule.c --- arrayfnsmodule.c 2000/01/20 18:18:28 1.1 +++ arrayfnsmodule.c 2000/03/28 17:55:04 @@ -575,6 +575,96 @@ return result ; } +static int +binary_searchf(float dval, float dlist [], int len) +{ + /* binary_search accepts three arguments: a numeric value and + * a numeric array and its length. It assumes that the array is sorted in + * increasing order. It returns the index of the array's + * largest element which is <= the value. It will return -1 if + * the value is less than the least element of the array. */ + /* self is not used */ + int bottom , top , middle, result; + + if (dval < dlist [0]) + result = -1 ; + else { + bottom = 0; + top = len - 1; + while (bottom < top) { + middle = (top + bottom) / 2 ; + if (dlist [middle] < dval) + bottom = middle + 1 ; + else if (dlist [middle] > dval) + top = middle - 1 ; + else + return middle ; + } + if (dlist [bottom] > dval) + result = bottom - 1 ; + else + result = bottom ; + } + + return result ; +} +/* return float, rather than double */ + +static PyObject * +arr_interpf(PyObject *self, PyObject *args) +{ + /* interp (y, x, z) treats (x, y) as a piecewise linear function + * whose value is y [0] for x < x [0] and y [len (y) -1] for x > + * x [len (y) -1]. An array of floats the same length as z is + * returned, whose values are ordinates for the corresponding z + * abscissae interpolated into the piecewise linear function. */ + /* self is not used */ + PyObject * oy, * ox, * oz ; + PyArrayObject * ay, * ax, * az , * _interp; + float * dy, * dx, * dz , * dres, * slopes; + int leny, lenz, i, left ; + + PyObject *tpo = Py_None; /* unused argument, we've already parsed it*/ + + Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz, &tpo)); + GET_ARR(ay,oy,PyArray_FLOAT,1); + GET_ARR(ax,ox,PyArray_FLOAT,1); + if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) { + SETERR ("interp: x and y are not the same length."); + Py_DECREF(ay); + Py_DECREF(ax); + return NULL ;} + GET_ARR2(az,oz,PyArray_FLOAT,1,MAX_INTERP_DIMS); + lenz = A_SIZE (az); + dy = (float *) A_DATA (ay); + dx = (float *) A_DATA (ax); + dz = (float *) A_DATA (az); + /* create output array with same size as 'Z' input array */ + Py_Try (_interp = (PyArrayObject *) PyArray_FromDims + (A_NDIM(az), az->dimensions, PyArray_FLOAT)); + dres = (float *) A_DATA (_interp) ; + slopes = (float *) malloc ( (leny - 1) * sizeof (float)) ; + for (i = 0 ; i < leny - 1; i++) { + slopes [i] = (dy [i + 1] - dy [i]) / (dx [i + 1] - dx [i]) ; + } + for ( i = 0 ; i < lenz ; i ++ ) + { + left = binary_searchf (dz [i], dx, leny) ; + if (left < 0) + dres [i] = dy [0] ; + else if (left >= leny - 1) + dres [i] = dy [leny - 1] ; + else + dres [i] = slopes [left] * (dz [i] - dx [left]) + dy [left]; + } + + free (slopes); + Py_DECREF(ay); + Py_DECREF(ax); + Py_DECREF(az); + return PyArray_Return (_interp); +} + static char arr_interp__doc__[] = "" ; @@ -592,8 +682,24 @@ PyArrayObject * ay, * ax, * az , * _interp; double * dy, * dx, * dz , * dres, * slopes; int leny, lenz, i, left ; + PyObject *tpo = Py_None; + char type_char = 'd'; + char *type = &type_char; - Py_Try(PyArg_ParseTuple(args, "OOO", &oy, &ox, &oz)); + Py_Try(PyArg_ParseTuple(args, "OOO|O", &oy, &ox, &oz,&tpo)); + if (tpo != Py_None) { + type = PyString_AsString(tpo); + if (!type) + return NULL; + if(!*type) + type = &type_char; + } + if (*type == 'f' ) { + return arr_interpf(self, args); + } else if (*type != 'd') { + SETERR ("interp: unimplemented typecode."); + return NULL; + } GET_ARR(ay,oy,PyArray_DOUBLE,1); GET_ARR(ax,ox,PyArray_DOUBLE,1); if ( (leny = A_SIZE (ay)) != A_SIZE (ax)) {