# %%
import numpy as np
import matplotlib.pyplot as plt

# %%
print(60 * "-")
print(__file__)
print("Aufgabe 2. Interpolationspolynome gemäss NEWTON-Schema berechnen")
# punkte [[x0, x1, x2, xn], [y0, y1, y2, yn]]
punkte = [
    [np.array([4, 8]), np.array([1, -1])],
    [np.array([-1, 1, 2]), np.array([15, 5, 9])],
    [np.array([-1, 0, 1, 2]), np.array([-5, -1, -1, 1])],
]

for punkt in punkte:

    # Parameter
    x_data = punkt[0]
    y_data = punkt[1]
    x_0 = x_data[0]
    x_E = x_data[-1]
    N = 201
    lw = 3
    fig = 1

    # Berechnung
    n = np.size(y_data)
    TAB = np.block([[y_data], [np.zeros((n - 1, n))]])
    c_data = np.zeros(n)
    c_data[0] = y_data[0]

    for i in range(1, n):
        for j in range(1, n):
            TAB[i][j] = (TAB[i - 1][j] - TAB[i - 1][j - 1]) / (
                x_data[j] - x_data[j - i]
            )
        c_data[i] = TAB[i][i]

    # Funktionen:
    def p(x):
        d = 1
        y = c_data[0]
        for k in range(1, n):
            d = d * (x - x_data[k - 1])
            y = y + c_data[k] * d
        return y

    # Daten
    u_data = np.linspace(x_0, x_E, N)
    v_data = p(u_data)

    fh = plt.figure(fig)
    plt.plot(u_data, v_data, linewidth=lw)
    plt.plot(x_data, y_data, "o", linewidth=lw)
    plt.xlabel(r"$x$")
    plt.ylabel(r"$y$")
    plt.grid(visible=True)
    plt.axis("image")

print(60 * "-")