import numpy as np
import matplotlib.pyplot as plt
# b
x = np.array([1850, 1875, 1900, 1925, 1950, 1975, 2000])
y = np.array([285.2, 288.6, 295.7, 305.3, 311.3, 331.36,
              369.64])
h = 25
n = len(x) - 1

z = np.zeros(n-1)
for i in range(n-1):
    z[i] = (y[i] - 2 * y[i+1] + y[i+2]) * 6 / (h**2)

M = np.zeros(n+1)
a = np.zeros(n)
b = np.zeros(n)
c = np.zeros(n)
d = np.zeros(n)

z_matrix = z

matrix = np.array([
    [5, 1, 0, 0, 0],
    [1, 4, 1, 0, 0],
    [0, 1, 4, 1, 0],
    [0, 0, 1, 4, 1],
    [0, 0, 0, 1, 5]
])

inverse_matrix = np.linalg.inv(matrix)
M_matrix = np.dot(inverse_matrix, z_matrix)
for i in range(n-1):
    M[i+1] = M_matrix[i]
M[0] = M[1]
M[n] = M[n-1]

for i in range(n):
    a[i] = (M[i+1] - M[i]) / (6 * h)
    b[i] = M[i] / 2
    c[i] = (y[i+1] - y[i]) / h - h * (M[i+1] + 2 * M[i]) / 6
    d[i] = y[i]

for i in range(n):
    print(f"S_{i+1}(x) = {a[i]:.4e}(x - {x[i]})^3 "
          f"+ {b[i]:.4e}(x - {x[i]})^2 + {c[i]:.4e}(x - {x[i]}) "
          f"+ {d[i]:.4e}, for {x[i]} <= x < {x[i+1]}")

year_estimates = [1990, 2010]
spline_estimates = np.zeros(2)
for j in range(2):
    year = year_estimates[j]
    i = int((year - 1850) / h)
    if i >= n:
        i = n - 1
    dx = year - x[i]
    spline_estimates[j] = (a[i]*dx**3 + b[i]*dx**2
                           + c[i]*dx + d[i])

print(f"\nCO2 concentration estimate for 1990: "
      f"{spline_estimates[0]:.2f} ppm")
print(f"CO2 concentration estimate for 2010: "
      f"{spline_estimates[1]:.2f} ppm")

print("a coefficients=", a)
print("b coefficients=", b)
print("c coefficients=", c)
print("d coefficients=", d)
print("M values=", M)

xx = np.linspace(x[0], x[-1], 1000)
yy = np.zeros(len(xx))
for k in range(len(xx)):
    year = xx[k]
    i = int((year - 1850) / h)
    if i >= n:
        i = n - 1
    dx = year - x[i]
    yy[k] = a[i]*dx**3 + b[i]*dx**2 + c[i]*dx + d[i]

plt.figure(figsize=(12, 6))
plt.plot(x, y, 'o', label='Original data')
plt.plot(xx, yy, label='Parabolic Runout Cubic Spline')
plt.xlabel('Year')
plt.ylabel('CO2 Concentration (ppm)')
plt.title('Parabolic Runout Cubic Spline Interpolation of CO2 Concentration')
plt.legend()
plt.grid(True)
plt.show()