# dernum.py

def dernum (x,y):
  derivate=[]; n=len(x)-1
  L=hermite(zip(x[:3],y[:3])); derivate.append(L.der()(x[0]))
  L=hermite(zip(x[:4],y[:4])); derivate.append(L.der()(x[1]))
  for k in xrange(2,n-1):
    L=hermite(zip(x[k-2:k+3],y[k-2:k+3])); derivate.append(L.der()(x[k]))
  L=hermite(zip(x[n-3:],y[n-3:])); derivate.append(L.der()(x[n-1]))
  L=hermite(zip(x[n-2:],y[n-2:])); derivate.append(L.der()(x[n]))
  return derivate

def interspline (x,y):
  derivate=dernum(x,y); n=len(x)-1
  def f (z):
    for k in xrange(n):
      x1=x[k]; x2=x[k+1]; y1=y[k]; y2=y[k+1]
      if z==x1: return y1
      elif z<x2:
        yp1=derivate[k]; yp2=derivate[k+1]
        theta=1/(x2-x1); v=y2-y1; dx1=z-x1; dx2=z-x2
        return y1+yp1*dx1+(v*theta-yp1)*theta*dx1*dx1+\
            (yp2-2*v*theta+yp1)*theta*theta*dx1*dx1*dx2
    return y2 # Solo per z=x[n].
  return f