def stineman_interp(xi,x,y,yp=None):
<EOL>x=np.asarray(x, np.float_)<EOL>y=np.asarray(y, np.float_)<EOL>assert x.shape == y.shape<EOL>N=len(y)<EOL>if yp is None:<EOL><INDENT>yp = slopes(x,y)<EOL><DEDENT>else:<EOL><INDENT>yp=np.asarray(yp, np.float_)<EOL><DEDENT>xi=np.asarray(xi, np.float_)<EOL>yi=np.zeros(xi.shape, np.float_)<EOL>dx = x[<NUM_LIT:1>:] - x[:-<NUM_LIT:1>]<EOL>dy = y[<NUM_LIT:1>:] - y[:-<NUM_LIT:1>]<EOL>s = dy/dx  <EOL>idx = np.searchsorted(x[<NUM_LIT:1>:-<NUM_LIT:1>], xi)<EOL>sidx = s.take(idx)<EOL>xidx = x.take(idx)<EOL>yidx = y.take(idx)<EOL>xidxp1 = x.take(idx+<NUM_LIT:1>)<EOL>yo = yidx + sidx * (xi - xidx)<EOL>dy1 = (yp.take(idx)- sidx) * (xi - xidx)       <EOL>dy2 = (yp.take(idx+<NUM_LIT:1>)-sidx) * (xi - xidxp1) <EOL>dy1dy2 = dy1*dy2<EOL>yi = yo + dy1dy2 * np.choose(np.array(np.sign(dy1dy2), np.int32)+<NUM_LIT:1>,<EOL>((<NUM_LIT:2>*xi-xidx-xidxp1)/((dy1-dy2)*(xidxp1-xidx)),<EOL><NUM_LIT:0.0>,<EOL><NUM_LIT:1>/(dy1+dy2),))<EOL>return yi<EOL>
Given data vectors *x* and *y*, the slope vector *yp* and a new
abscissa vector *xi*, the function :func:`stineman_interp` uses
Stineman interpolation to calculate a vector *yi* corresponding to
*xi*.

Here's an example that generates a coarse sine curve, then
interpolates over a finer abscissa::

  x = linspace(0,2*pi,20);  y = sin(x); yp = cos(x)
  xi = linspace(0,2*pi,40);
  yi = stineman_interp(xi,x,y,yp);
  plot(x,y,'o',xi,yi)

The interpolation method is described in the article A
CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell
W. Stineman. The article appeared in the July 1980 issue of
Creative Computing with a note from the editor stating that while
they were:

  not an academic journal but once in a while something serious
  and original comes in adding that this was
  "apparently a real solution" to a well known problem.

For *yp* = *None*, the routine automatically determines the slopes
using the :func:`slopes` routine.

*x* is assumed to be sorted in increasing order.

For values ``xi[j] < x[0]`` or ``xi[j] > x[-1]``, the routine
tries an extrapolation.  The relevance of the data obtained from
this, of course, is questionable...

Original implementation by Halldor Bjornsson, Icelandic
Meteorolocial Office, March 2006 halldor at vedur.is

Completely reworked and optimized for Python by Norbert Nemec,
Institute of Theoretical Physics, University or Regensburg, April
2006 Norbert.Nemec at physik.uni-regensburg.de