'''
Version: 12/03/25
--> Version definitiva del analisis por integracion numerica
'''

# %% LIBRERIAS

from tqdm import tqdm
import numpy as np
import scipy.integrate as spi
import scipy.interpolate as interp

# %% METODOS DE INGRACION NUMERICAS (TRAPZ,CUADRATURA CADA N PUNTOS)

def gauss_quad_discrete_complex(x, y, n=2):

    integral = 0.0 + 0.0j  # Iniciar como número complejo

    for i in range(len(x) - 1):
        x_local = x[i:i+n+1]  
        y_local = y[i:i+n+1]

        # Ajustar un polinomio de Lagrange
        poly = interp.BarycentricInterpolator(x_local, y_local)

        # Integrar usando quad_vec, que maneja números complejos
        result = spi.quad_vec(poly, x[i], x[i+1])[0]  

        integral += result

    return integral

def trapecios_complejo(x, y):

    h = np.diff(x)  # Diferencias entre los puntos
    integral = np.sum(h * (y[:-1] + y[1:]) / 2)  # Aplicamos la fórmula del trapecio
    return integral

# %% CALCULO (NO SE CUENTA EL MODO M = 0)

def f2(MatrizSph,limitemodos,positive="Yes"):
    
    # DEFINICION DE PARAMETROS
    
    time = MatrizSph.shape[0]
    R = MatrizSph[:,:,0]  # MATRIZ TIEMPO X PARTICULAS
    PHI = MatrizSph[0,:,2]  # Vector angular
    
    # CAMBIO DE COORDENADAS
    
    VMR = np.mean(R)
    U = R - VMR  # matriz
    Z = PHI*VMR  # vector
    dz = abs(Z[1]-Z[0]) # util para trapz (dx)
    
    # DATOS
    
    N = int(limitemodos) # digamos 51
    n = np.linspace(1,N,N) # nunca se incluye el modo 0
    modes = n
    cn = np.zeros((time, N), dtype=complex) # coeficientes
    
    for i in tqdm(range(time), desc="Calculando coeficientes"):
        for j in range(N):
            cn[i,j] = np.trapz(U[i, :] * np.exp(-1j * (n[j] / VMR) * Z),dx = dz)
     
            #cn[i, j] = trapecios_complejo(U[i, :] * np.exp(-1j * (n[j] / VMR) * Z), Z)  # método de trapecios complejos
   
    F1 = (abs(cn))**2  # módulo de cn elevado al cuadrado
    VMF1 = np.mean(F1, axis=0)  # valor medio del módulo de cn al cuadrado
    
    F2 = abs(cn)  # módulo de cn
    VMF2 = np.mean(F2, axis=0)  # valor medio de los coeficientes
    VMF3 = VMF2**2  # valor medio del módulo de coeficientes --> todo elevado al cuadrado
    
    diference = VMF1 - VMF3
    
    VM_MODUCUA = (1/(2*np.pi*VMR)) * diference  # Con esto compararemos qx tal como dice la fórmula
    qx = modes / VMR
    
    return qx,VM_MODUCUA