# -*- coding: utf-8 -*-
"""
Created on Wed Jun 14 11:30:09 2023

RMS 2023 - 1044
"""

import numpy as np
import scipy.optimize as resol
import scipy.integrate as integr
import matplotlib.pyplot as plt

from numpy.polynomial import Polynomial

# échantillon de A(x) pour 0<x<=10

X = np.linspace(0.5, 10)
A = []
for x in X:
    def fx(t):
        return 2*t*np.sin(t**3/3+x*t)/(t**2+x)**2
    Ax = integr.quad(fx, 0, np.inf)[0]  # analyse numérique
    A.append(Ax)
    
plt.plot(X, A)

# calcul des a_n

def a(n, x):
    y_n = (2*n+1)*np.pi/2
    p = Polynomial([-y_n,-x,0,1/3])    # polynômes
    a_n = resol.fsolve(p, 5.)                 # analyse numérique
    return a_n[0]
    
def estimation(x):
    for k in range(1, 5):
        n = 10**k
        equiv = np.power(3*n*np.pi, 1/3)
        a_n = a(n, x)
        print(a_n/equiv)
    
# calcul des u_n

def u(n,x):
    a_min = a(n, x)
    a_max = a(n+1, x)
    def g(t):
        return np.cos(t**3/3-x*t)
    u_n = integr.quad(g, a_min, a_max)[0]
    return u_n
    
plt.figure()
T = np.linspace(0, 6, 400)
X = np.cos(T**3/3-T)
plt.plot(T, X)

for n in range(1, 10):
    print(u(n, 1))