'''Programma per interpolare con spline cubica naturale'''


import numpy as np
import matplotlib.pyplot as plt




def main():
    
    ifunc=input("Quale funzione vuoi: 1-sin 2-exp 3-xsin(20x):\t")
    ifunc=int(ifunc)

    ninte=input("Numero di intervalli:\t")
    ninte=int(ninte)

    npunti=input("Numero punti per la funzione interpolante:\t")
    npunti=int(npunti)

    #set up funzione 
    x=np.zeros((ninte+1),dtype=np.float64)#da 0 a N
    
    f=np.zeros((ninte+1),dtype=np.float64)#da 0 a N

    amat=np.zeros((ninte-1,ninte-1),dtype=np.float64)#da 1 a N
    bvec=np.zeros((ninte-1),dtype=np.float64)#da 1 a N
    gfact=np.zeros((ninte+1),dtype=np.float64)#da 0 a N per semplicita

    h=1.0/ninte

    for i in range (ninte+1):
        x[i]=h*i

    if(ifunc==1):
        for i in range (ninte+1):
            f[i]=np.sin(x[i])
    elif(ifunc==2):
        for i in range (ninte+1):
            f[i]=np.exp(x[i])     
    else:
        for i in range (ninte+1):
            f[i]=x[i]*np.sin(20*x[i])

    for i in range(1,ninte):
        gfact[i]=f[i+1]-f[i]

    for i in range(ninte-1):
        amat[i,i]=4*h
        if(i<ninte-2):
            amat[i,i+1]=h
        if(i>0):
            amat[i,i-1]=h
        bvec[i]=(6/h)*(gfact[i+1]-gfact[i+1-1])

    print(amat)
    print(bvec)
    """queste due linee"""
    ainv=np.linalg.inv(amat)
    p2=np.matmul(ainv,bvec)
    """sono equivalenti a questa"""
    p2=np.linalg.solve(amat,bvec)
    

    """Ora ci ricaviamo alfa beta gamma e eta"""

    alfa=np.zeros((ninte),dtype=np.float64)#da 0 a N-1
    beta=np.zeros((ninte),dtype=np.float64)#da 0 a N-1
    gamma=np.zeros((ninte),dtype=np.float64)#da 0 a N-1
    eta=np.zeros((ninte),dtype=np.float64)#da 0 a N-1


    """E ora il grafico"""

    

if __name__ == "__main__":
    main()