# -*- coding: utf-8 -*-
"""
Created on Tue May 21 12:21:11 2019

@author: HAL900
"""

import numpy as np
import matplotlib.pyplot as plt

def baryGew(x):
    n = len(x)
    lambdas = np.ones_like(x)
    for i in range(n):
        for k in range(n):
            if k == i:
                continue
            lambdas[i] *=  1 / (x[i]-x[k])
    return lambdas

def bary(x, f, lam, X):
    n = len(x)
    Y = np.zeros_like(X)

    for j in range(len(X)):
        sum1 = 0
        sum2 = 0
        zero = False
        for i in range(n):
            if X[j] == x[i]:
                Y[j] = f[i]
                zero = True
            else:
                sum1 += lam[i]*f[i] / (X[j]-x[i])
                sum2 += lam[i] / (X[j]-x[i])

        if not zero:
            Y[j] = sum1 / sum2

    return Y

def f(x):
    return 1 / (1 + 25*np.power(x,2))

def test():
    xn = np.linspace(-1,1,401)

    xi_10 = [2*i / 10 - 1 for i in range(11)]
    xi_100 = [2*i / 100 - 1 for i in range(101)]
    xi2_10 = [np.cos( (2*i+1)*np.pi / (2*10+2) ) for i in range(11)]
    xi2_100 = [np.cos( (2*i+1)*np.pi / (2*100+2) ) for i in range(101)]

    b10 = bary(xi_10, f(xi_10), baryGew(xi_10), xn)
    b100 = bary(xi_100, f(xi_100), baryGew(xi_100), xn)
    plt.figure(1)
    plt.plot(xn, f(xn), label="f")
    plt.plot(xn, b10, label="p,10")
    plt.plot(xn, b100, label="p,100")
    plt.ylim(-1,1)
    plt.legend()

    b10 = bary(xi2_10, f(xi2_10), baryGew(xi2_10), xn)
    b100 = bary(xi2_100, f(xi2_100), baryGew(xi2_100), xn)
    plt.figure(2)
    plt.plot(xn, f(xn), label="f")
    plt.plot(xn, b10, label="p,10")
    plt.plot(xn, b100, label="p,100")
    plt.ylim(-1,1)
    plt.legend()


test()