'''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(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

    

    for i in range(1,ninte-1):
        alfa[i]=p2[i-1+1]/(6*h)
        beta[i]=-p2[i-1]/(6*h)
        gamma[i]=-p2[i-1+1]*h/6+f[i+1]/h
        eta[i]=p2[i-1]*h/6-f[i]/h

    alfa[0]=p2[0]/(6*h)
    beta[0]=0.0
    gamma[0]=-p2[0]*h/6+f[1]/h
    eta[0]=-f[0]/h


    alfa[ninte-1]=0
    beta[ninte-1]=-p2[ninte-2]/(6*h)
    gamma[ninte-1]=+f[ninte]/h
    eta[ninte-1]=p2[ninte-2]*h/6-f[ninte-1]/h


    """E ora il grafico"""

    xgraf=np.zeros((npunti),dtype=np.float64)
    fgraf=np.zeros((npunti),dtype=np.float64)
    fspline=np.zeros((npunti),dtype=np.float64)

    hgraf=1.0/npunti

    for i in range(npunti):
        xgraf[i]=hgraf*i
        m=int(xgraf[i]//h)
        fspline[i]=alfa[m]*(xgraf[i]-x[m])**3+beta[m]*(xgraf[i]-x[m+1])**3+gamma[m]*(xgraf[i]-x[m])+eta[m]*(xgraf[i]-x[m+1])

    if(ifunc==1):
        for i in range (npunti):
            fgraf[i]=np.sin(xgraf[i])
    elif(ifunc==2):
        for i in range (npunti):
            fgraf[i]=np.exp(xgraf[i])     
    else:
        for i in range (npunti):
            fgraf[i]=xgraf[i]*np.sin(20*xgraf[i])


    fig, ax = plt.subplots()

   
    plt.plot(xgraf,fgraf)
    plt.plot(xgraf,fspline)
    plt.plot(x, f, marker='o', markersize=6, markerfacecolor='red',
         markeredgecolor='black', markeredgewidth=2, linestyle='')
    
    plt.show(block=True)
    

if __name__ == "__main__":
    main()