'''Programma per interpolare con fft  e calcolare la derivata'''


import numpy as np
import scipy as sp
import matplotlib.pyplot as plt




def main():
    
    ifunc=input("Quale funzione vuoi: 1-triangolo 2-gradino 3:sin:\t")
    ifunc=int(ifunc)

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

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

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


    h=1.0/ninte

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

    if(ifunc==1):
        for i in range (ninte):
            if(x[i]<0.5):
                f[i]=x[i]
            else:
                f[i]=1.0-x[i]
    elif(ifunc==2):
        for i in range (ninte):
            if(x[i]>=0.25 and x[i]<=0.75):
                f[i]=1.0
            else:
                f[i]=0.0
    else:
        for i in range (ninte):
            f[i]=np.sin(2*np.pi*x[i])
    """Trasformata di Fourier"""

    g=sp.fft.fft(f)

    print(f)
    print(g)
    """Adesso riempi il grafico"""

    xg=np.zeros((npunti),dtype=np.float64)#da 0 a Npunti
    fg=np.zeros((npunti),dtype=np.float64)#da 0 a Npunti
    fftg=np.zeros((npunti),dtype=np.float64)#da 0 a Npunti
    derg=np.zeros((npunti),dtype=np.float64)#da 0 a Npunti

    hg=1.0/npunti

    for i in range (npunti):
        xg[i]=hg*i

    if(ifunc==1):
        for i in range (npunti):
            if(xg[i]<0.5):
                fg[i]=xg[i]
            else:
                fg[i]=1.0-xg[i]
    elif(ifunc==2):
        for i in range (npunti):
            if(xg[i]>=0.25 and xg[i]<=0.75):
                fg[i]=1.0
            else:
                fg[i]=0.0
    else:
        for i in range (npunti):
            fg[i]=np.sin(2*np.pi*xg[i])



    for i in range (npunti):
        for j in range(int(ninte/2)): 
            fftg[i]+=np.real(g[j]*np.exp(2j*np.pi*j*xg[i]))
        for j in range(1,int(ninte/2)): 
            fftg[i]+=np.real(g[ninte-j]*np.exp(-2j*np.pi*j*xg[i]))
    fftg*=1.0/ninte

    
    """E ora il grafico"""

    fig, ax = plt.subplots()

   
    plt.plot(xg,fg)
    plt.plot(xg,fftg)
    
    plt.show(block=True)

    
    """Adesso la derivata"""

    for i in range (npunti):
        for j in range(int(ninte/2)): 
            fftg[i]+=np.real(g[j]*(j*1j)*np.exp(2j*np.pi*j*xg[i]))
        for j in range(1,int(ninte/2)+1): 
            fftg[i]+=np.real(g[ninte-j]*(-j*1j)*np.exp(-2j*np.pi*j*xg[i]))

    fig, ax = plt.subplots()

    fftg*=1.0/ninte
    plt.plot(xg,fftg)
    
    plt.show(block=True)


if __name__ == "__main__":
    main()