# %% LIBRERIAS

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from Codigo3D_FuncionGeneral import calculate_al

import Codigo3D_AjusteHelfrich as c3

# %% FUNCION DE AJUSTE (Para generar continuo de interpolacion)

kb = 1.380649e-23  # Constante de Boltzmann
T = 300  # Temperatura en Kelvin
kbT = kb * T  # Factor kb * T

def fx(x,kappa,b):
    y = kbT/(8*b + kappa*x) # NOTA x incluye r0 y el termino de l 
    return y

# %% ANALISIS

Datos_simul_cart = np.load('matrizcartesianas.npy') 

Grid,ProyecGrid,terminol,a2l,r0 = calculate_al(Datos_simul_cart,20)

# %% DATOS PARA GRAFICAS

minimo = 2 # corresponde a l = 2
maximo = 6 # corresponde a l = 5, pero luego sumamos 1 y es l = 6

coefs = a2l[minimo:maximo+1] # coeficientes entre l = 2, l = 6
valorx = terminol[minimo:maximo+1]/(r0*r0) # x de la funcion

kappa,b = c3.aproxHelfrich(valorx,coefs,T) 

bending = kappa/kbT
stretching = b

print("-------------------")
print(f"Para modos entre l = {minimo} y l = {maximo} se obtiene:")
print("Bending (en uds kbT): ",bending)
print(" (Aprox) Stretching [N/m]: ",stretching)

# %% AJUSTE PAPER 

lvector = np.linspace(0,20,21) # vector de l 
lfiltro = np.linspace(2,6,100) # vector de l entre l = 2, l = 6

terminolfiltro = lfiltro*(lfiltro+1)*(lfiltro-1)*(lfiltro+2) # termino en l filtrado
valorxfiltro = terminolfiltro/(r0*r0) # x filtrado

Y = fx(valorxfiltro,kappa,b) # se genera el continuo de lo que nos da el ajuste (generado con lfiltro!)

# %% GRAFICAS

#### GRAFICA 1: Representacion de los coeficientes en funcion del termino en l

plt.figure()
plt.plot(valorxfiltro*r0**2,Y,linestyle="-.",color = "red") # misma dim
plt.plot(valorx*r0**2,coefs,linestyle = "none", color = "blue",marker = "o")  # misma dim
plt.xlabel("l*(l+1)*(l-1)*(l+2)",fontweight = "bold")
plt.ylabel(r'$\langle |a_\ell|^2 \rangle$',fontweight = "bold")
plt.title("Coeficientes en funcion de l*(l+1)*(l-1)*(l+2)",fontsize = 14,fontweight = "bold")

#### GRAFICA 2: Representacion de los coeficientes en funcion de l (todos)

plt.figure()
plt.plot(lfiltro,Y,linestyle="-.",color = "red")
plt.plot(lvector,a2l,linestyle = "none", color = "blue",marker = "o")
plt.yscale('log')
plt.xlabel("l",fontweight = "bold")
plt.ylabel(r'$\langle |a_\ell|^2 \rangle$',fontweight = "bold")
plt.xticks(np.arange(0, 22, 2))
plt.title("Coeficientes en funcion de l",fontsize = 14,fontweight = "bold")
plt.show()