interp hw

import numpy
import pylab

def lpw(sizer,dx,dy):
        lpw_array = numpy.zeros(dx.size)
        qq = 0
        for sort_index in range(1, sizer):
            while sizer[sort_index] > x[qq]:
                qq += 1
                xa= dx[qq-1]
                ya= dy[qq-1]
                xb= dx[qq]
                yb= dx[qq]
                yy = (yb-ya)*((xx-xa)/(xb-xa))+ ya
                lpw_array[sort_index] = yy
        return lpw_array

def csp(sizer,dx,dy):
        csp_array = numpy.zeros(dx.size)
        for sort_index in range(1, sizer):
            while sizer[sort_index] > dx[sort_index]:
                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 = (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
            csp_array = a0 + a1*x_s + a2*x_s**2 + a3*x_s**3
        return csp_array

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

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

sizer = numpy.arange(0,4.60,0.01)
sort_index = numpy.argsort(dx)
dx = dx[sort_index]
dy = dy[sort_index]

the_function = numpy.cos(sizer)*numpy.tanh(sizer)

lpw_y = lpw(sizer,dx,dy)
csp_y = csp(sizer,dx,dy)

pylab.figure(1)
pylab.plot(dx,dy)
pylab.plot(sizer,lpw_y,'o')
pylab.plot(sizer,csp_y,'^')
pylab.plot(sizer,the_function,'--')
pylab.show()