interp hw

import numpy
import pylab

data = numpy.loadtxt('data_interp.dat')

dx = data[:,0]
dy = data[:,1]

sort_index = numpy.argsort(dx)

dx = dx[sort_index]
dy = dy[sort_index]

x = numpy.arange(0,4.60,.01)

for i in range (0,x.size):
while x > dx[i]:
	i += 1

	dddy = dy[i+1]

	ddy = dy[i]

	ddyy = dy[i-1]

	ddyyy = dy[i-2]

	dddx = dx[i+1]

	ddx = dx[i]

	ddxx = dx[i-1]

	ddxxx = dx[i-2]

	y_pw =


            xa= x[i-1]
            ya= y[i-1]
            xb= x[i]
            yb= x[i]
            yy= (yb-ya)*((xx-xa)/(xb-xa))+ ya

	y = (ddy-ddyy)*(x-ddxx)/(ddx-ddxx) + ddyy

	mA = ((ddy-ddyy)*(ddxx-ddxxx)**2 + (ddyy-ddyyy)*(ddx-ddxx)**2)/((ddx-ddxx)*(ddxx-ddxxx)*(ddx-ddxxx))

	mB = ((dddy-ddy)*(ddx-ddxx)**2 + (ddy-ddyy)*(dddx-ddx)**2)/((dddx-ddx)*(ddx-ddxx)*(dddx-ddxx))

	x_s = (x - ddxx)/(ddx - ddxx)

	a0 = ddyy
	a1 = (ddx - ddxx)*mA
	a2 = 3*(ddy - ddyy) - 2*(ddx-ddxx)*mA - (ddx - ddxx)*mB
	a3 = -2*(ddy - ddyy) + (ddx-ddxx)*mA + (ddx - ddxx)*mB

y2 = a0 + a1*x_s + a2*x_s**2 + a3*x_s**3

the_function = numpy.cos(numpy.arange(0,4.60,0.01))*numpy.tanh(numpy.arange(0,4.60,0.01))

pylab.figure(1)
pylab.plot(dx,dy)
pylab.plot(x,y_pw,'o')
pylab.plot(x,y_cs,'^')
pylab.plot(x,the_function,'--')
pylab.show()