#!/usr/bin/env -S python3
import numpy
from numpy.polynomial import Polynomial, Chebyshev
from matplotlib.pyplot import *
from scipy.signal import windows
from scipy import interpolate

measured = numpy.asarray([
    [1.0000, 0.9241, 0.9978, 0.9802],
    [0.9666, 0.9209, 0.9966, 0.9774],
    [0.9330, 0.9174, 0.9952, 0.9746],
    [0.9000, 0.9135, 0.9936, 0.9715],
    [0.8666, 0.9098, 0.9916, 0.9680],
    [0.8333, 0.9057, 0.9893, 0.9646],
    [0.8000, 0.9011, 0.9867, 0.9609],
    [0.7666, 0.8968, 0.9834, 0.9569],
    [0.7333, 0.8920, 0.9801, 0.9528],
    [0.7000, 0.8869, 0.9764, 0.9483],
    [0.6666, 0.8813, 0.9724, 0.9433],
    [0.6333, 0.8753, 0.9681, 0.9384],
    [0.6000, 0.8691, 0.9632, 0.9325],
    [0.5666, 0.8622, 0.9580, 0.9264],
    [0.5333, 0.8551, 0.9522, 0.9199],
    [0.5000, 0.8471, 0.9464, 0.9129],
    [0.4666, 0.8370, 0.9393, 0.9051],
    [0.4333, 0.8234, 0.9321, 0.8964],
    [0.4000, 0.8053, 0.9239, 0.8871],
    [0.3666, 0.7829, 0.9143, 0.8763],
    [0.3333, 0.7557, 0.9038, 0.8641],
    [0.3000, 0.7223, 0.8917, 0.8500],
    [0.2666, 0.6789, 0.8777, 0.8295],
    [0.2333, 0.6224, 0.8610, 0.7973],
    [0.2000, 0.5434, 0.8399, 0.7504],
    [0.1666, 0.4629, 0.7992, 0.6823],
    [0.1333, 0.3816, 0.7318, 0.5687],
    [0.1000, 0.3009, 0.6064, 0.4406],
    [0.0666, 0.2207, 0.4530, 0.3140],
    [0.0333, 0.1400, 0.2403, 0.1872],
    [0.0000, 0.0625, 0.0649, 0.0629],
])

xs  = measured[:,0]
xs_ = numpy.linspace(-0.10, 1.15, 100)

tukey = windows.tukey(measured.shape[0]);

int_r = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,1]), extrapolate=True)
int_g = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,2]), extrapolate=True)
int_b = interpolate.PchipInterpolator(numpy.flip(xs), numpy.flip(measured[:,3]), extrapolate=True)

off_r = int_r(0)
off_g = int_g(0)
off_b = int_b(0)

k = 1./max( int_r(1)-off_r, int_g(1)-off_g, int_b(1)-off_b )

fit_r = Polynomial.fit(xs_, k*(int_r(xs_)-off_r), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))
fit_g = Polynomial.fit(xs_, k*(int_g(xs_)-off_g), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))
fit_b = Polynomial.fit(xs_, k*(int_b(xs_)-off_b), deg=16, domain=[0,1], window=[0,1], w=windows.tukey(xs_.shape[0]))

def c_arr_print(label, arr):
    print( label+"[] = {" + ', '.join([f'\n\t{v:e}' if 0 == i%4 else f'{v:e}' for i,v in enumerate(arr)]) + "};\n" )

c_arr_print('GLfloat const coef_r', fit_r.convert().coef[1:])
c_arr_print('GLfloat const coef_g', fit_g.convert().coef[1:])
c_arr_print('GLfloat const coef_b', fit_b.convert().coef[1:])

def polyeval(c, xs):
    C_ = numpy.flip(c)
    ys = numpy.ones(xs.size)*C_[0]
    for c in C_[1:]:
        ys = ys*xs + c
    return ys*xs

if True:
    plot(xs, measured[:,1], '+')
    plot(xs_, fit_r(xs_))
    plot(xs_, polyeval(fit_r.convert().coef[1:], xs_), '*')

    plot(xs, measured[:,2], '+')
    plot(xs_, fit_g(xs_))
    plot(xs_, polyeval(fit_g.convert().coef[1:], xs_), '*')

    plot(xs, measured[:,3], '+')
    plot(xs_, fit_b(xs_))
    plot(xs_, polyeval(fit_b.convert().coef[1:], xs_), '*')

    xlim(0,1)
    ylim(0,1)
    show()