import numpy as np
import matplotlib.pyplot as plt
from cycler import cycler

from splineA import splineA


def fig5_18():
    a = -1
    b = 1

    # interpolation data
    nx = 13
    xd = np.linspace(a, b, nx)
    yd = np.zeros(nx)
    for iy in range(nx):
        yd[iy] = 1 / (1 + 25 * xd[iy] ** 2)

    # exact values
    n = 400
    x = np.linspace(a, b, n)
    ye = np.zeros(n)
    for ie in range(n):
        ye[ie] = 1 / (1 + 25 * x[ie] ** 2)

    #### evaluate and then plot spline function
    # ysM = MATLAB spline function
    # ys = natural spline (requires splineA.m file)

    ys = splineA(xd, yd, x)
    # ysM = spline(xd,yd1,x);

    # p_n(x)
    for ix in range(nx):
        cc = (2 * (ix + 1) - 1) * np.pi / (2 * nx)
        xd[ix] = 0.5 * (a + b + (b - a) * np.cos(cc))
        yd[ix] = f(xd[ix])
    p = np.zeros(n)
    for ii in range(n):
        for k in range(nx):
            p[ii] = p[ii] + yd[k] * ell(k, x[ii], xd)

    plt.close("all")
    # get(gcf)
    fig, ax = plt.subplots(figsize=(9.0, 3.2))
    fig.subplots_adjust(left=0.06, right=0.98, bottom=0.21, top=0.97)
    co = np.array([
        [0, 0, 1],
        [0, 0.5, 0],
        [1, 0, 0],
        [0, 0.75, 0.75],
        [0.75, 0, 0.75],
        [0.75, 0.75, 0],
        [0.25, 0.25, 0.25],
    ])
    ax.set_prop_cycle(cycler(color=co))

    ax.plot(x, ys, "r", linewidth=1.6)
    ax.plot(x, p, "b", linewidth=1.6)
    ax.plot(x, ye, "--k", linewidth=1.6)
    ax.set_yticks([0, 1])
    ax.set_xticks([-1, -0.5, 0, 0.5, 1])
    ax.set_xlim(-1, 1)
    ax.set_ylim(-0.02, 1.05)
    ax.grid(True)
    ax.set_xlabel("x-axis")
    ax.legend(
        [" Spline", " Cheb", " Runge"],
        loc="upper right",
        prop={"size": 14, "weight": "bold"},
    )
    ax.tick_params(labelsize=14)
    for label in ax.get_xticklabels() + ax.get_yticklabels():
        label.set_fontweight("bold")
    for spine in ax.spines.values():
        spine.set_visible(True)
    plt.show()

    # plt.savefig('/Users/mark/Desktop/splineR.eps')


def f(x):
    g = 1 / (1 + 25 * x ** 2)
    return g


# ell_i(x) function
def ell(i, x, xd):
    n2 = xd.size
    p = 1.0
    for j in range(n2):
        if j != i:
            p = p * (x - xd[j]) / (xd[i] - xd[j])
    return p


if __name__ == "__main__":
    fig5_18()
