import numpy as np
from scipy.interpolate import barycentric_interpolate

def runge(x):
    """Función de Runge."""
    return 1 / (1 + x ** 2)

N = 11  # Nodos de interpolación

xp = np.arange(11) - 5  # -5, -4, -3, ..., 3, 4, 5
fp = runge(xp)

x = np.linspace(-5, 5, num=1000)

y = barycentric_interpolate(xp, fp, x)

%matplotlib inline

import matplotlib.pyplot as plt

l, = plt.plot(x, y)
plt.plot(x, runge(x), '--', c=l.get_color())
plt.plot(xp, fp, 'o', c=l.get_color())
leg = plt.legend(['Interpolación', 'Real'])
leg.get_frame().set_facecolor('#fafafa')

from numpy.polynomial import chebyshev

coeffs_cheb = [0] * 11 + [1]
T11 = chebyshev.Chebyshev(coeffs_cheb, [-5, 5])

xp_ch = T11.roots()
fp_ch = runge(xp_ch)

y_ch = barycentric_interpolate(xp_ch, fp_ch, x)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

l1, = ax1.plot(x, y)
ax1.plot(x, runge(x), '--', c=l1.get_color())
ax1.plot(xp, fp, 'o', c=l1.get_color())
ax1.set_title("Nodos equiespaciados")

leg1 = ax1.legend(['Interpolación', 'Real'])
leg1.get_frame().set_facecolor('#fafafa')

l2, = ax2.plot(x, y_ch)
ax2.plot(x, runge(x), '--', c=l2.get_color())
ax2.plot(xp_ch, fp_ch, 'o', c=l2.get_color())
ax2.set_ylim(ax1.get_ylim())
ax2.set_title("Nodos de Chebyshev")

leg2 = ax2.legend(['Interpolación', 'Real'])
leg2.get_frame().set_facecolor('#fafafa')

from scipy.interpolate import InterpolatedUnivariateSpline

# Pato
P = [(0.9, 1.3), (1.3, 1.5), (1.9, 1.8), (2.1,2.1), (2.6, 2.6), (3.0, 2.7),
     (3.9, 2.3), (4.4, 2.1), (4.8, 2.0), (5.0, 2.1), (6, 2.2), (7, 2.3),
     (8, 2.2), (9.1, 1.9), (10.5, 1.4), (11.2, 0.9), (11.6, 0.8), (12, 0.6),
     (12.6, 0.5), (13, 0.4), (13.2, 0.2)]
xi, yi = zip(*P)

print("Número de puntos: {}".format(len(xi)))

x = np.linspace(min(xi), max(xi), num=1001)

y1d = np.interp(x, xi, yi)
ysp = InterpolatedUnivariateSpline(xi, yi)(x)

plt.plot(xi, yi, 'x', mew=2)
plt.plot(x, y1d)
plt.plot(x, ysp)
leg = plt.legend(['Puntos', 'Lineal', 'Cúbico'])
leg.get_frame().set_facecolor('#fafafa')

import numpy.polynomial as P

# Cargamos los datos
data = np.loadtxt("polar.dat")
C_L, C_D = data

# Descarto los datos que no me sirven
stall_idx = np.argmax(C_L)

y = C_D[:stall_idx + 1]
x = C_L[:stall_idx + 1] ** 2

# Ajuste lineal, devuelve los coeficientes en orden creciente
C_D0, k = P.polynomial.polyfit(x, y, deg=1)

print(C_D0, k)

C_L_dom = np.linspace(C_L[0], C_L[stall_idx], num=1001)
C_D_int = P.polynomial.polyval(C_L_dom ** 2, (C_D0, k))

l, = plt.plot(C_D_int, C_L_dom)
plt.plot(C_D, C_L, 's', mew=2, mec=l.get_color(), mfc="none")
plt.title("Curva polar")
plt.xlabel("$C_D$")
plt.ylabel("$C_L$")
leg = plt.legend(["Ajuste parabólico", "Datos reales"], 4)
leg.get_frame().set_facecolor('#fafafa')

def func(x, A, B, C):
    """Modelo para nuestros datos."""
    return A * np.exp(-B * x ** 2) + C

from scipy.optimize import curve_fit

Adat, Bdat, Cdat = 2.5, 1.3, 0.5

xdat = np.linspace(-2, 4, 12)
ydat = func(xdat, Adat, Bdat, Cdat) + 0.2 * np.random.normal(size=len(xdat))

(A, B, C), _ = curve_fit(func, xdat, ydat)
print(A, B, C)

x = np.linspace(-2, 4, 1001)

l, = plt.plot(x, func(x, Adat, Bdat, Cdat), '--', label="Sin ruido")
plt.plot(xdat, ydat, 'o', mew=2, mfc="none", mec=l.get_color(), label="Datos")
plt.plot(x, func(x, A, B, C), label="Ajuste")
plt.ylim(0, 3.5)
leg = plt.legend()
leg.get_frame().set_facecolor('#fafafa')