#!/usr/bin/python

# ----------------------------------------------------------------
# John Kerl
# kerl.john.r@gmail.com
# 2005-01-11
#
# This is an example of Lagrange interpolation, using degree 3.
#
# Given x_0 .. x_3, y_0 .. y_3, compute the polynomial
#
#              (x - x_1)(x - x_2)(x - x_3)
# y_0 ------------------------------------ +
#              (x_0-x_1)(x_0-x_2)(x_0-x_3)
#
#     (x - x_0)         (x - x_2)(x - x_3)
# y_1 ------------------------------------ +
#     (x_1-x_0)         (x_1-x_2)(x_1-x_3)
#
#     (x - x_0)(x - x_1)         (x - x_3)
# y_2 ------------------------------------ +
#     (x_2-x_0)(x_2-x_1)         (x_2-x_3)
#
#     (x - x_0)(x - x_1)(x - x_2)
# y_3 ------------------------------------
#     (x_3-x_0)(x_3-x_1)(x_3-x_2)
#
# Program output:
#
#   0.00000 -0.34510  1.00000 -1.3451e+00
#   0.00391 -0.32857  1.00391 -1.3325e+00
#   0.00781 -0.31211  1.00784 -1.3200e+00
#   ...
#   0.98828  2.68467  2.68661 -1.9441e-03
#   0.99219  2.69585  2.69713 -1.2788e-03
#   0.99609  2.70705  2.70768 -6.3083e-04
#   1.00000  2.71828  2.71828  0.0000e+00
#   1.00391  2.72953  2.72892  6.1388e-04
#   1.00781  2.74081  2.73960  1.2110e-03
#   1.01172  2.75212  2.75032  1.7915e-03
#   ...
#   1.23828  3.46298  3.44968  1.3299e-02
#   1.24219  3.47640  3.46318  1.3215e-02
#   1.24609  3.48986  3.47674  1.3125e-02
#   1.25000  3.50337  3.49034  1.3029e-02
#   1.25391  3.51693  3.50400  1.2927e-02
#   1.25781  3.53054  3.51772  1.2819e-02
#   1.26172  3.54419  3.53149  1.2705e-02
#   ...
#   1.48828  4.43016  4.42948  6.8910e-04
#   1.49219  4.44727  4.44681  4.5811e-04
#   1.49609  4.46445  4.46422  2.2838e-04
#   1.50000  4.48169  4.48169  0.0000e+00
#   1.50391  4.49900  4.49923 -2.2694e-04
#   1.50781  4.51639  4.51684 -4.5236e-04
#   1.51172  4.53384  4.53452 -6.7618e-04
#   ...
#   1.98828  7.30223  7.30297 -7.4339e-04
#   1.99219  7.33106  7.33155 -4.9811e-04
#   1.99609  7.36000  7.36025 -2.5028e-04
#   2.00000  7.38906  7.38906  0.0000e+00
#   2.00391  7.41823  7.41798  2.5265e-04
#   2.00781  7.44752  7.44701  5.0758e-04
#   2.01562  7.50644  7.50542  1.0239e-03
#   ...
#   2.48828 12.04297 12.04056  2.4043e-03
#   2.49219 12.08932 12.08769  1.6277e-03
#   2.49609 12.13583 12.13500  8.2641e-04
#   2.50000 12.18249 12.18249  0.0000e+00
#   2.50391 12.22932 12.23017 -8.5189e-04
#   2.50781 12.27631 12.27804 -1.7296e-03
#   2.51172 12.32346 12.32610 -2.6336e-03
#   ...
#   2.98828 19.39450 19.85153 -4.5704e-01
#   2.99219 19.46416 19.92923 -4.6507e-01
#   2.99609 19.53404 20.00723 -4.7319e-01
#
# Note that the fit is good at the interpolation points, but worse in between
# them and even worse at the edges.
#
# When the data is plotted, examine the error curve and note that it has
# appoximately the classic W shape of a quartic.  This makes sense, since we've
# used a cubic fit; the next-order term is quartic.
# ================================================================
# This software is released under the terms of the GNU GPL.
# Please see LICENSE.txt in the same directory as this file.
# ================================================================

from __future__ import division # 1/2 = 0.5, not 0.
from math     import *
from kerlutil import *
import lagrinterp_m

# ----------------------------------------------------------------
# Input mesh:
xs = [0.20, 0.40, 0.60, 0.80]

# Function to be approximated:
f = atan

# Interpolation mesh:
xlo = 0.0
xhi = 1.0
dx  = 1.0 / 100

# Coefficients for interpolating polynomial:
cs = lagrinterp_m.li_make_coeffs_from_xs_and_f(xs, f)

# Evaluate the function and the interpolant over the mesh:
for x in mfrange(xlo, dx, xhi):
	y = f(x)
	z = lagrinterp_m.li_eval(x, xs, cs)
	print "%10.5f %10.5f %10.5f %10.3e" % (x, y, z, y-z)