# lagSymSine.py
# Andrew Davison, ad@coe.psu.ac.th, June 2026
'''
Plots f(x) = sin(x) and the Lagrange polynomial for that
function based on the interpolation nodes in xs. The 
polynomial is generated using Sympy's interpolating_poly()
and the result is printed in Algebraic form to stdout in
various forms.

The function and polynomial are plotted between x = 0 to \pi
in steps of 1/1000 to show how the choice of interpolation
nodes can affect its accuracy by causing it to oscillate.

The absolute error between the two plots is also shown as
a separate graph.
'''

import math

import sympy as sp
from sympy.polys.specialpolys import interpolating_poly
   # https://docs.sympy.org/latest/modules/polys/reference.html

import matplotlib.pyplot as plt


x = sp.symbols('x')

# Interpolation nodes: x_k = k*pi/3, k = 0..3
xs = [k*sp.pi/3 for k in range(4)]
ys = [sp.sin(v) for v in xs]

# Interpolating polynomial using Lagrange
ipoly = interpolating_poly(len(xs), x, X=xs, Y=ys)
print("Poly:")
print(ipoly); print()

p = sp.expand(ipoly)   # multiply out all the terms
print("Expanded:")
print(p); print()
print("Factored:")
print(sp.factor(p)); print()

# Convert to a Python function
pFunc = sp.lambdify(x, p, "math")

# Generate plot data
xVals = [(i * math.pi) / 1000 for i in range(1001)]
sinVals = [math.sin(v) for v in xVals]
polyVals = [pFunc(v) for v in xVals]
errVals = [abs(poly - sine)
                  for poly, sine in zip(polyVals, sinVals)]

# Plot curves and error
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6,7), sharex=True)

# Top plot: sin(x) and interpolation
ax1.plot(xVals, sinVals, label="sin(x)")
ax1.plot(xVals, polyVals, "--", label="Lagrange")

# Plot interpolation nodes
nodeXs = [float(v.evalf()) for v in xs]
nodeYs = [float(v.evalf()) for v in ys]
ax1.plot(nodeXs, nodeYs, "o", label="Data")

ax1.axhline(0, color="black")
ax1.axvline(0, color="black")
ax1.set_ylabel("y")
ax1.grid(True)
ax1.legend()
ax1.set_title("sin(x) vs Lagrange Approx")

# Bottom plot: absolute error
ax2.plot(xVals, errVals, label="|Lagrange - sin(x)|")

ax2.axhline(0, color="black")
ax2.axvline(0, color="black")
ax2.set_xlabel("x")
ax2.set_ylabel("abs error")
ax2.grid(True)
ax2.legend()
ax2.set_title("Absolute Error")

plt.tight_layout()
plt.show()