# -*- coding: utf-8 -*-
"""
Created on Mon May 16 17:11:30 2011

@author: FRANCES
"""
import numpy as np
import scipy.linalg as las
import matplotlib.pylab as plt


def lagrange(xnodi,fnodi,x):
     """
     Funzione che determina in un insieme di punti
     il valore del polinomio interpolante ottenuto
     dalla formula di Lagrange e la funzione di Lebesgue
 
     f,L = lagrange(xnodi,fnodi,x)
 
     Dati di input:
 
     xnodi vettore con i nodi dell'interpolazione
     fnodi vettore con i valori nei nodi
     x   vettore con i punti in cui si vuole
          calcolare il valore del polinomio
          interpolante
     Dati di output:
     f   vettore contenente i valori assunti
          dal polinomio interpolante
     L  vettore contenente i valori assunti dalla
         funzione di lebesgue in x
 
     """
     n=xnodi.shape[0]
     m=x.shape[0]

     l=np.zeros([n,m])
     ind=np.zeros([n,n-1],dtype=int)
     
     ind[0,:]=np.arange(1,n)
    
     for i in range(1,(n)):
        ind[i,0:(n-1)]=np.floor(np.concatenate((np.arange(0,(i)),np.arange(i+1,n))))
     ind[n-1,:]=np.arange(0,n-1)
    

     for i in range(0,n):
         den = np.prod( xnodi[i]-xnodi[(ind[i,:])])
         for j in range(0,m):
             l[i,j]=np.prod( x[j] - xnodi[ind[i,:]] )/den
    


     y = np.sum(np.dot(np.diag(fnodi),l),axis=0)
     

     L=np.sum(abs(l),axis=0)

     return  y,L

def potenze(xnodi,fnodi,xx):
     """
     Funzione che determina in un insieme di punti
     il valore del polinomio interpolante ottenuto
     dalla formula di Lagrange e la funzione di Lebesgue
 
     f = potenze(xnodi,fnodi,x)
 
     Dati di input:
 
     xnodi vettore con i nodi dell'interpolazione
     fnodi vettore con i valori nei nodi
     x   vettore con i punti in cui si vuole
          calcolare il valore del polinomio
          interpolante
     Dati di output:
     f   vettore contenente i valori assunti
          dal polinomio interpolante
 
     """
     n=xnodi.shape[0]
     A=np.zeros([n,n])
     for j in range(0,n):
        A[:,j]=xnodi**(j)

     p = las.solve(A,fnodi)
     f = np.polyval(p[np.arange(n-1,-1,-1)],xx)
     return f
     
def  cheby(a,b,n):
     """
     Nodi di Chebyshev
     """
     c =  (a + b + (b-a)*np.cos((2*(np.arange(0,n+1))+1)*np.pi/(2*(n+1))))/2
     return c
     
def runge(x):
    """
    Funzione di Runge
    """
    y=1/(1+25*x**2)
    return y
    
def plotinterp(f=runge,a=-1,b=1,n=8,tn=0):
    """
    plot del polinomio interpolante di grado n la funzione
    f in [a,b] usando n+1 nodi equidistanti se tn=0
             n+1 nodi di Chebyshev se tn=1
    """
    
    if (tn==0):
       xnodi = np.linspace(a,b,n+1)
    else:
       xnodi = cheby(a,b,n) 
       
    fnodi = f(xnodi)
    xx = np.linspace(a,b,100)

    yl, L = lagrange(xnodi,fnodi,xx)
    ye = f(xx)
    error = np.max(np.abs(yl-ye))
    plt.figure(1)
    plt.cla()
    plt.plot(xx,yl,xnodi,fnodi,'o',xx,ye,'--',markersize=12,linewidth=2) 
    plt.figure(2)
    plt.cla()
    plt.plot(xx,L,linewidth=2)
    plt.show()
    return  error
            
    
    
    
    