import matplotlib.pyplot as plt

x = [15.5, 12.5, 8, 10, 7, 4, 8, 10, 9.5, 14, 18, 17, 22, 25, 19, 24.5, 23, 17, 16, 12.5, 16.5, 21, 17, 11, 5.5, 7.5, 10, 12]
y = [32.5, 28.5, 29, 33, 33, 37, 39.5, 38.5, 42, 43.5, 42, 40, 41.5, 37, 35, 33.5, 29.5, 30.5, 32, 19.5, 24.5, 22, 15, 10.5, 2.5, 8, 14.5, 20]
t = [k/27 for k in range(len(x))]


def gauss(t_, x_):
    n = len(t_)
    h = [0 for k in range(n-1)]
    b = [0 for k in range(n-1)]
    for i in range(n-1):
        h[i] = t_[i+1] - t_[i]
        b[i] = 6*(x_[i+1]-x_[i])/h[i]
    u = [0 for k in range(n-1)]
    v = [0 for k in range(n-1)]
    u[0] = 2*(h[0] + h[1])
    v[0] = b[1] - b[0]
    for i in range(1,n-1):
        u[i] = 2*(h[i-1] + h[i]) - h[i-1]*h[i-1]/u[i-1]
        v[i] = b[i] - b[i-1] - h[i-1]*v[i-1]/u[i-1]
    z = [0 for k in range(n)]
    for i in range(n-2,0,-1):
        z[i] = (v[i] - h[i]*z[i+1])/u[i]
    return z


z_x = gauss(t,x)
z_y = gauss(t,y)


def s(x_,y_,z_):
    i = 0
    n = len(t)
    for k in range(n-1,0,-1):
        if x_-k >= 0:
            i = k
    h_i = t[i+1]-t[i]
    s_i = z_[i]/(6*h_i)*(t[i+1]-x_)**3 + z_[i+1]/(6*h_i)*(x_-t[i])**3 \
          + (y_[i+1]/h_i - z_[i+1]*h_i/6)*(x_- t[i])+ (y_[i]/h_i - z_[i]*h_i/6)*(t[i+1] - x_)
    return s_i

M=1000
s_u_x = [s(k/M,x,z_x) for k in range(M)]
s_u_y = [s(k/M,y,z_y) for k in range(M)]
plt.plot(s_u_x, s_u_y)
plt.show()