import matplotlib.pyplot as plt

from itertools import zip_longest
from math import sin, cos
from math import pi

#kod dla metody Hermite'a, parametry działają analogicznie.

x0 = 1/3
xn = 3.0
n = 30
czebyszew = 0

def getOrigValues(x):
    return x * sin(3 * pi / x)

def getDerivValues(x):
    return sin( 3 * pi / x) - 3 * pi * cos(3 * pi / x) / x

def calculateBasicPoints(n, czebyszew):
    points = []
    interval = (xn - x0) / (n - 1)
    for x in range(n):
        if czebyszew:
            currX = 1/2 * (xn + x0) + 1/2 * (xn - x0) * cos((2 * (x + 1) - 1) / 2 / n * pi)
        else:
            currX = x0 + x * interval
        points.append([currX, getOrigValues(currX)])
    return points

def maxDev(l1, l2):
    temp = 0
    for x,y in zip(l1, l2):
        if abs(x - y) > temp: temp = abs(x - y)
    return temp

def makePyramid(xVals, yVals, derVals):
    n = len(xVals)

    pyramid = []
    for x,y in zip(xVals, yVals):
        pyramid.append([x, y])
    for x in range(int(n/2)):
        pyramid[x * 2 + 1].append(derVals[x])
    for x in range(int(n/2) - 1):
        pyramid[x * 2 + 2].append( (pyramid[x * 2 + 2][1] - pyramid[x * 2 + 1][1]) / (xVals[x * 2 + 2] - xVals[x * 2 + 1]) )

    for x in range(2, n):
        for y in range(x - 1):
            pyramid[x].append( (pyramid[x][2 + y] - pyramid[x - 1][2 + y]) / ( xVals[x] - xVals[x - (2 + y)] ))

    coeffs = []
    for x in range(n):
        coeffs.append(pyramid[x][x + 1])
    return coeffs

def getHermiteVal(coefficients, x, xVals):
    points = list(reversed(xVals))
    coeffs = list(reversed(coefficients))
    total = coeffs[0]
    for i in range(1, len(coeffs)):
        total *= (x - points[i])
        total += coeffs[i]
    return total


def do(n):
    points = calculateBasicPoints(n, czebyszew)
    xVal = [x for x, y in points]
    yVal = [y for x, y in points]
    xVals = []
    yVals = []
    derVals = [getDerivValues(x) for x in xVal]
    for x in xVal:
        xVals.append(x)
        xVals.append(x)
    for y in yVal:
        yVals.append(y)
        yVals.append(y)

    coefficients = makePyramid(xVals, yVals, derVals)

    point_count = 500
    pointX = []
    pointY = []
    pointYorig = []
    interval = (xn - x0) / point_count
    for x in range(point_count):
        currX = x0 + x * interval
        pointX.append(currX)
        pointY.append(getHermiteVal(coefficients, currX, xVals))
        pointYorig.append(getOrigValues(currX))

    print(n, maxDev(pointY, pointYorig))


    plt.plot(pointX, pointY, label="wielomian interpolujący")


    plt.plot(pointX, pointYorig, label="oryginalna funkcja")
    for point in calculateBasicPoints(n, czebyszew):
        plt.scatter(point[0], point[1], s=50)
    plt.legend(loc=1)
    plt.show()

do(n)