from numpy import *
import pylab
from matplotlib import mlab
from scipy.linalg import solve
from scipy import array
import matplotlib.pyplot as plt 
from pylab import *
f1=array([
         [4.0,102.56],
         [4.2,104.18],
         [4.5,100.11]])
def sweep(f,d2fa,d2fb):
    k=2
    A = zeros(len(f))
    B = zeros(len(f))
    C = zeros(len(f))
    F = zeros(len(f))
    res = zeros(len(f))
    alpha = zeros(len(f))
    beta = zeros(len(f))
    while (k<len(f)):
        A[k]=h(f,k-1)
        B[k]=h(f,k)
        C[k]=2*(h(f,k-1)+h(f,k))
        F[k]=3*(func(f,k)-func(f,k-1))
        k+=1
    alpha[2]=-B[2]/C[2]
    beta[2]=F[2]/C[2]
    k=2
    while (k<len(f)-1):
        alpha[k+1]=-B[k]/(A[k]*alpha[k]+C[k])
        beta[k+1]=(F[k]-A[k]*beta[k])/(A[k]*alpha[k]+C[k])                
        k+=1
    k=len(f)-2
    res[0]=d2fa
    res[len(f)-1]=d2fb
    while(k>0):
        res[k]=alpha[k+1]*res[k+1]+beta[k+1]    
        k-=1
    return res        
             
def h(f,k):
    return f[k,0]-f[k-1,0]
def func(f,k):
    return (f[k,1]-f[k-1,1])/h(f, k)             
             
def sweep2(f,d2fa,d2fb):
    k=2
    sigma = zeros(10)
    lamda = zeros(10)
    C=zeros(10)
    C[0]=d2fa
    C[9]=d2fb
    sigma[0]=-h(f, 2)/(2*h(f,1)+2*h(f,2))
    lamda[0]=3*(func(f,2)-func(f,1))/(2*(h(f,1)+h(f,2)))
    while k<len(f):
        sigma[k-1]=h(f,k)/(2*h(f,k-1)+2*h(f,k)+h(f,k-1)*sigma[k-2])
        lamda[k-1]=(3*(func(f,k)-func(f,k-1))-h(f,k-1)*lamda[k-2])/(2*h(f,k-1)+2*h(f,k)+h(f,k-1)*sigma[k-2])
        k+=1
    k=9
    while k>1:
        C[k-1]=sigma[k-1]*C[k]+lamda[k-1]
        k-=1
    return C
    
def ispline(f,x,C):
    if x<f[0,0]:
        print "viberite bol'she x"
    elif x>f[len(f)-1,0]:
        print "viberite men'she x"
    else: 
        k=0
        count=0
        A = f[:,1]
        B = zeros(len(f))
        D = zeros(len(f))
        while (k<len(f)):
            if (x-f[k,0]<=0):
                count=k
                k=10
            else:
                k+=1
        k=1
        while (k<len(f)):
            B[k]= func(f,k)+2/3*h(f,k)*C[k]+1/3*h(f,k)*C[k-1]
            D[k]= (C[k]-C[k-1])/(3*h(f,k))
            k+=1
        g=A[count]+B[count]*(x-f[count,0])+C[count]*(x-f[count,0])**2+D[count]*(x-f[count,0])**3
        return g
def graf(f):   
    xmin = 4.0
    xmax = 4.5
    dx = 0.01
    xlist = mlab.frange (xmin, xmax, dx)
    ylist = [ispline(f,x,sweep2(f,0,0)) for x in xlist]
    pylab.plot (xlist, ylist)
    y2list = [ispline(f,x,sweep(f,0,0)) for x in xlist]
    pylab.plot (xlist, y2list)
    x = f[:,0]
    y = f[:,1]
    area = 100
    scatter(x,y,s=area, marker='^', c='r')
    pylab.show()
    
graf(f1)