Subida inicial del proyecto

Codigo2D_AjusteHelfrich.py

@@ -0,0 +1,40 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+
+kb = 1.380649e-23  # Constante de Boltzmann
+T = 300  # Temperatura en Kelvin
+kbT = kb * T  # Factor kb * T
+
+def helfrich_model(x, bend, stretch):
+    
+    bend = abs(bend)
+    stretch = abs(stretch)
+    
+    term = stretch / bend + x**2  
+
+    result = kbT * (0.5 / stretch) * ((1 / x) - term ** -0.5)
+    return result
+
+def aproxHelfrich(x, y):
+    x = np.array(x)
+    y = np.array(y)
+
+    try:
+
+        popt, _ = curve_fit(helfrich_model, x, y, p0=[1.4e-19, 1e-7]) # valores referencia
+        bend, stretch = popt
+        
+        # Aseguramos que bend y stretch sean siempre positivos
+        bend = abs(bend)
+        stretch = abs(stretch)
+
+        # # Establecemos un límite inferior para el valor de 'bend'
+        # bend = max(bend, 1e-20)
+        return bend, stretch
+    except RuntimeError:
+        print("Error: La optimización no convergió")
+        return None, None
+
+

Codigo2D_Base.py

@@ -0,0 +1,84 @@
+'''
+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
+   
+  
+

Codigo2D_CodigoBaseFiltro.py

@@ -0,0 +1,293 @@
+'''
+- Codigos para filtrar datos simulados en un intervalo de -z a z (valor)
+- Graficas para matrices con ceros y cambio de coordenadas en estas matrices
+- Interpolacion de puntos 
+'''
+
+import numpy as np
+import orthopoly
+import numpy as np
+import time as timeglobal
+from scipy.spatial import KDTree
+from tqdm import tqdm
+from colorama import Fore, Style
+from scipy.interpolate import splprep, splev
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+# from Codigo1 import mcart_msph
+
+def public_c2s(posicionescart):
+    '''
+    De cartesianas a esfericas (c2s) de forma publica
+    '''
+    x,y,z       = posicionescart.T
+    r,theta,phi = orthopoly.spherical_harmonic.cart2sph(-x,-y,z) #Bug in library
+    return np.asarray([r,theta,phi]).T
+
+def public_s2c(posicionesesf):
+    '''
+    De esfericas a cartesianas (s2c) de forma publica
+    '''
+    r,theta,phi = posicionesesf.T
+    x,y,z       = orthopoly.spherical_harmonic.sph2cart(r,theta,phi)
+    return np.asarray([x,y,z]).T
+
+def mcart_msph(matrizcart):
+    tiempo = matrizcart.shape[0]
+    matrizsph = (public_c2s(matrizcart.reshape(-1, 3))).reshape(tiempo, -1, 3)
+    return matrizsph
+
+def msph_mcart_filtrado(matrizsph):
+    """
+    Transforma las filas de una matriz de coordenadas cartesianas a esféricas,
+    ignorando las filas con valores [0, 0, 0].
+
+    Parameters:
+    matrizcart (numpy.ndarray): Matriz de coordenadas cartesianas [tiempo, n, 3].
+
+    Returns:
+    numpy.ndarray: Matriz transformada a coordenadas esféricas, 
+                   manteniendo las filas [0, 0, 0] sin cambiar.
+    """
+    tiempo = matrizsph.shape[0]
+    matrizcart = np.zeros_like(matrizsph)  
+
+    for t in range(tiempo):
+        frame = matrizsph[t]
+        
+        mask = ~np.all(frame == 0, axis=1) # solo seleccionamos filas que sean distintas de 0 (sino sale error en arcos)
+
+        matrizcart[t, mask] = public_s2c(frame[mask])
+
+    return matrizcart
+
+def mcart_msph_filtrado(matrizcart):
+    """
+    Transforma las filas de una matriz de coordenadas cartesianas a esféricas,
+    ignorando las filas con valores [0, 0, 0].
+
+    Parameters:
+    matrizcart (numpy.ndarray): Matriz de coordenadas cartesianas [tiempo, n, 3].
+
+    Returns:
+    numpy.ndarray: Matriz transformada a coordenadas esféricas, 
+                   manteniendo las filas [0, 0, 0] sin cambiar.
+    """
+    tiempo = matrizcart.shape[0]
+    matrizsph = np.zeros_like(matrizcart)  
+
+    for t in range(tiempo):
+        frame = matrizcart[t]
+        
+        mask = ~np.all(frame == 0, axis=1) # solo seleccionamos filas que sean distintas de 0 (sino sale error en arcos)
+
+        matrizsph[t, mask] = public_c2s(frame[mask])
+
+    return matrizsph
+
+def contarpuntos(valor,matrix):
+
+    
+    tiempo, particulas, _ = matrix.shape
+
+    contador = np.zeros(tiempo, dtype=int)
+
+    for t in range(tiempo):
+
+        z_coords = matrix[t, :, 2]
+        contador[t] = np.sum((-valor <= z_coords) & (z_coords <= valor))
+      
+    return contador
+
+def datosfiltrados(valor,matrix):
+    
+    contador = contarpuntos(valor,matrix)
+    
+    fecha_actual = timeglobal.strftime("%d-%b-%Y %H:%M:%S")
+    print(f"{fecha_actual} - Análisis3 - Contando datos con Z ∈ [-{valor}, {valor}] ...")
+    print("El numero de puntos esta entre ",np.min(contador)," y ",np.max(contador))
+    print(f"{fecha_actual} - Análisis3 - Filtrando datos de las simulaciones...")
+
+    tiempo = matrix.shape[0]
+    max_contador = np.max(contador)
+    matriz_homogenea = np.zeros((tiempo, max_contador, 3), dtype=float)
+
+    for t in range(tiempo):
+        puntos = matrix[t, :, :]
+        puntos_en_rango = puntos[(-valor <= puntos[:, 2]) & (puntos[:, 2] <= valor)]
+        matriz_homogenea[t, :len(puntos_en_rango), :] = puntos_en_rango
+
+    return matriz_homogenea
+
+def distribuirtraj_sobregrid(traj_sph, grid_sph):
+    
+    '''
+    Funcion analoga a la de distribucion de Codigo1
+    '''
+    fecha_actual = timeglobal.strftime("%d-%b-%Y %H:%M:%S")
+    n_traj = traj_sph.shape[1]
+    n_grid_points = grid_sph.shape[0] 
+    
+    print("------------------------------")
+    print(f"{fecha_actual} - Análisis3 - Distribuyendo {Style.BRIGHT + Fore.GREEN}{n_traj}{Style.RESET_ALL} puntos de trayectoria en {Style.BRIGHT + Fore.BLUE}{n_grid_points}{Style.RESET_ALL} puntos del grid.")
+
+    print(f"Se está usando {Style.BRIGHT + Fore.GREEN}{traj_sph.shape[0]}{Style.RESET_ALL} valores de tiempo")
+
+    
+    theta_grid = grid_sph[:, 1]
+    phi_grid = grid_sph[:, 2]
+
+    time = traj_sph.shape[0]
+    # n_grid_points = grid_sph.shape[0]
+
+    grid_kdtree = KDTree(grid_sph)
+
+    traj_grid_sph = np.empty((time, n_grid_points, 3))
+
+    for t, frame_pos in enumerate(tqdm(traj_sph, total=time, desc="Procesando")): # Para cada instante de tiempo
+        
+        frame_pos = frame_pos[~np.all(frame_pos == 0, axis=1)] # Filtrar valores distintos de cero
+
+        if frame_pos.size == 0:
+            raise Exception(f"Ningun punto encontrado para el tiempo {t}. Revisar los datos de trayectoria.")
+
+        r, theta, phi = frame_pos.T
+
+        r = np.ones_like(r) # radio 1
+        traj_over_sphere = np.vstack([r, theta, phi]).T
+
+        frame_grid_positions = grid_kdtree.query(traj_over_sphere)[1]
+
+        for i in range(n_grid_points):
+            mask = (frame_grid_positions == i)
+
+            if not np.any(mask):
+                raise Exception(f"Punto del grid {i} sin datos. Pruebe a reducir el numero de puntos del grid.")
+
+            r_vals, theta_vals, phi_vals = frame_pos[mask].T
+
+            traj_grid_sph[t, i] = np.array([np.mean(r_vals), theta_grid[i], phi_grid[i]])
+
+    return traj_grid_sph # devuelve en esfericas
+
+# %% GRAFICAS
+
+def graficapuntual(matriz_homogenea, t, rangoz,title):
+    """
+    Grafica los puntos almacenados en la matriz homogénea en 3D para un instante de tiempo específico,
+    ignorando los valores cero.
+
+    Parameters:
+    matriz_homogenea (numpy.ndarray): Matriz de dimensiones [tiempo, max_contador, 3].
+    t (int): Instante de tiempo para el cual se desea graficar los puntos.
+    """
+    if t < 0 or t >= matriz_homogenea.shape[0]:
+        raise ValueError("El tiempo especificado está fuera del rango válido.")
+    
+    puntos = matriz_homogenea[t]
+    
+    puntos_filtrados = puntos[~np.all(puntos == 0, axis=1)]
+    
+    if len(puntos_filtrados) > 0:
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        
+        x = puntos_filtrados[:, 0]
+        y = puntos_filtrados[:, 1]
+        z = puntos_filtrados[:, 2]
+        
+        ax.scatter(x, y, z, color = "blue", linestyle='None')
+        ax.set_xlabel('X')
+        ax.set_ylabel('Y')
+        ax.set_zlabel('Z')
+        ax.set_zlim(-rangoz, rangoz)
+        ax.set_title(f"{title} en t = {t}", fontweight="bold")
+        
+        plt.show()
+    else:
+        print(f"No hay puntos válidos para el tiempo {t}.")
+        
+def gengraficav2(matrix, t, title, valor):
+    '''
+    Funcion identica a gengrafica solo que ahora tiene un input para los limites en z
+    --> Util para comprobar que la funcion de filtrado actue adecuadamente
+    '''
+
+    a3 = matrix[t][:, 0]
+    b3 = matrix[t][:, 1]
+    c3 = matrix[t][:, 2]
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+    ax.scatter(a3, b3, c3, marker='o', color="blue")
+
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
+    ax.set_title(f"{title} en t = {t}", fontweight="bold")
+
+    ax.set_zlim(-valor, valor)
+
+    plt.show()
+    
+# %% CONTINUO
+
+# from Codigo1 import coefs,gencircunferencia
+
+# def expandir(MSph,npoints,Lmax):
+#     '''
+#     # Es como interpolar mediante armonicos esfericos (no funciona!)
+#     '''
+#     _,coe = coefs(MSph,Lmax)
+#     Matrizcart,Matrizsph = gencircunferencia(coe,npoints) # selecciono las esfericas en concreto
+#     return Matrizcart,Matrizsph # cartesianas
+
+# def interpolarpuntos(matrizcart, n):
+#     '''
+#     Interpola los puntos en coordenadas cartesianas para cada instante de tiempo.
+#     ver: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html
+#     '''
+    
+#     print("------------------")
+#     fecha_actual = timeglobal.strftime("%d-%b-%Y %H:%M:%S")
+#     print(f"{fecha_actual} - Análisis3 - Interpolando {matrizcart.shape[1]} .")
+#     print(f"Se han obtenido {n} puntos tras la interpolación.")
+
+    
+#     tiempo = matrizcart.shape[0]
+#     M1 = np.empty((tiempo, n+1, 3)) # cerrada con termino repetido
+#     minterpolada_cart = np.empty((tiempo,n,3)) # lo que buscamos
+#     u_new = np.linspace(0, 1, n+1) # siempre es cte, creamos un punto mas de lo normal (para forzar cierre) 
+#     znew = np.zeros((n+1)) # vector nulo (ya que se supone que el input es este mismo)
+    
+#     for t in range(tiempo):
+
+#         frame = matrizcart[t]
+
+#         mask = ~np.all(frame == 0, axis=1)
+#         frame = frame[mask]
+
+#         if frame.shape[0] < 2:
+#             raise ValueError(f"No hay suficientes puntos para interpolar en el tiempo {t}.")
+        
+#         # Datos sin cerrar
+#         x = frame[:,0]
+#         y = frame[:,1]
+        
+#         # Forzamos cierre
+#         x_closed = np.append(x, x[0])
+#         y_closed = np.append(y, y[0])
+        
+#         tck,_ = splprep([x_closed, y_closed], s=0) # tck generado a partir de los valores cerrados
+
+#         newpoints = splev(u_new, tck) # generamos puntos de dimension n+1, pero el ultimo termino es el primero
+
+#         M1[t] = np.vstack([newpoints[0], newpoints[1], znew]).T
+        
+#         minterpolada_cart[t][:,:] = M1[t][0:-1,:] # no contamos el ultimo termino (repetido)
+    
+#         if np.mean(abs(minterpolada_cart[t, :, 2])) <= 0.00001: # para que la grafica salga exactamente en 0
+#             minterpolada_cart[t,:,2] = 0
+        
+#     minterpolada_sph = mcart_msph(minterpolada_cart)
+#     return minterpolada_cart,minterpolada_sph

Codigo2D_Coefs.py

@@ -0,0 +1,51 @@
+'''
+Calculo de bending y stretching sobre datos en una circunferencia (se asume equiespaciamiento angular)
+Ultima Fecha Modificacion: 30/03/25
+'''
+
+# %% LIBRERIAS
+
+import numpy as np
+import matplotlib.pyplot as plt
+import Codigo2D_Base as c7
+import Codigo2D_AjusteHelfrich as c8
+
+# %% FUNCION TEORICA
+
+kb = 1.380649*10**-23
+T = 300
+kbT = kb*T
+
+valor = 10*kbT # este es el valor que deberia tener el bending
+
+def fx(x,stretch,bend):
+    y = kbT*(0.5/stretch)*((1/x) - (stretch/bend + x**2)**(-0.5))
+    return y
+
+# %% CODIGO 
+
+def calculatecoefs(MSph, nmin, nmax):
+    
+    '''
+    Se utiliza el metodo de trapecios
+    '''
+    
+    qx,u = c7.f2(MSph,21,positive = "Yes") # TODOS LOS COEFICIENTES (U) del 1 al 20
+    
+    '''
+    El bending solo se calcula para ciertos modos
+    '''
+    
+    qxnew = qx[nmin:nmax]
+    unew = u[nmin:nmax]
+    
+    bend,stretch = c8.aproxHelfrich(qxnew,unew)
+    qxcont = np.linspace(qxnew[0],qxnew[-1],100) # limites coinciden con qx, pero Y no tiene porque es un ajuste no hay restriccion
+    Ycont = fx(qxcont,stretch,bend)  # Al representar sera en el eje x: qxnew y en el y: Y
+    
+    print("----------------")
+    print(f"Para modos entre {nmin+1} y {nmax+1} se obtiene:")
+    print("Bending (en uds kbT): ",bend/kbT)
+    print("Stretching [N/m]: ",stretch)
+    
+    return qx,u,qxcont,Ycont,bend,stretch

Codigo2D_Ejemplo.py

@@ -0,0 +1,70 @@
+'''
+Ejemplo de implementacion del Ánalisis en el plano ecuatorial
+'''
+
+import numpy as np
+import matplotlib.pyplot as plt
+# from Codigo1 import mcart_msph,msph_mcart,gengraficasph
+import Codigo2D_FiltroZ as c1
+import Codigo2D_Interpolacion as c2
+import Codigo2D_Coefs as c3 # Este incluye a ajuste de Helfrich y a Base
+
+time = 50 # tiempo para graficar
+
+# %% PROGRAMA ARIN 
+
+normalizacion = 7*10**5 # No es necesaria son uds de los datos
+MATRIZ = np.load('matrizcartesianas.npy') # DATOS PROGRAMA ARIN
+
+MATRIZESF = c1.mcart_msph(MATRIZ)
+MATRIZESF[:,:,0] = MATRIZESF[:,:,0] / normalizacion # cambiamos el radio / normalizamos (Ver luego que tambien intorud)
+MATRIZCART = c1.msph_mcart(MATRIZESF)
+
+# %% FILTRO EN Z
+
+RMatrix = MATRIZESF[:,:,0] # radios
+VMR = np.mean(RMatrix) # Valor medio del radio
+
+Z1 = (VMR/10)*0.1 # Una centesima del radio medio
+FiltroCart,FiltroEsf = c1.filtroz(MATRIZ,Z1,normalizacion) # Filtramos
+
+# %% INTERPOLACION 
+
+MSph = c2.newinterpolate2(FiltroEsf,2001) # Interpolamos
+MCart = c1.msph_mcart(MSph)
+MCart[:,:,2] = 0
+
+# %% CALCULO DE BENDING Y STRETCH
+nmax = 19
+nmin = 4 # NOTA nmin = 0, es el modo 1
+qx,u,qxnew,unew,bend,stretch = c3.calculatecoefs(MSph,nmin,nmax)
+
+# %% GRAFICA
+
+# plt.figure()
+# plt.loglog(qx,u,marker = "o", color = "blue",linestyle = "none")
+# plt.loglog(qxnew,unew,color="red", linestyle="-.")
+
+plt.figure()
+plt.plot(qx[nmin:nmax],u[nmin:nmax],marker = "o",color = "blue")
+plt.plot(qxnew,unew,color="red", linestyle="-.")
+
+# GRAFICA FILTRO (Simplemente para ver forma)
+
+mask = FiltroEsf[time,:,0] != 0
+FiltroEsfmask = FiltroEsf[time,mask] # Eliminar 0's para ver mejor la grafica
+
+plt.figure()
+
+plt.plot(MSph[time, :, 2], MSph[time, :, 0], marker=".", linestyle="--", color="blue", label=r"Interpolación datos filtro ($\Delta z_1$)")
+
+z1_label = f"$\Delta z_1$ = {np.format_float_scientific(Z1, precision=2)}"
+
+plt.plot(FiltroEsfmask[:,2], FiltroEsfmask[:,0], marker="h", linestyle="none", color="red", label=z1_label)
+
+plt.title("Datos simulaciones filtro", fontweight="bold")
+plt.xlabel("qx", fontweight="bold")
+plt.ylabel(r"$\langle \mathbf{|u|^2} \rangle$", fontsize=14) 
+
+plt.legend()
+plt.show()

Codigo2D_FiltroZ.py

@@ -0,0 +1,70 @@
+'''
+Codigo para seleccionar en un rango de +-z
+Los datos de filtro estan en el plano z = 0
+'''
+import orthopoly
+import numpy as np
+import Codigo2D_CodigoBaseFiltro as c4 
+# from Codigo1 import msph_mcart,mcart_msph
+
+def public_c2s(posicionescart):
+    '''
+    De cartesianas a esfericas (c2s) de forma publica
+    '''
+    x,y,z       = posicionescart.T
+    r,theta,phi = orthopoly.spherical_harmonic.cart2sph(-x,-y,z) #Bug in library
+    return np.asarray([r,theta,phi]).T
+
+def public_s2c(posicionesesf):
+    '''
+    De esfericas a cartesianas (s2c) de forma publica
+    '''
+    r,theta,phi = posicionesesf.T
+    x,y,z       = orthopoly.spherical_harmonic.sph2cart(r,theta,phi)
+    return np.asarray([x,y,z]).T
+
+def msph_mcart(matrizsph):
+    tiempo = matrizsph.shape[0]
+    matrizcart = (public_s2c(matrizsph.reshape(-1, 3))).reshape(tiempo, -1, 3)
+    return matrizcart
+
+def mcart_msph(matrizcart):
+    tiempo = matrizcart.shape[0]
+    matrizsph = (public_c2s(matrizcart.reshape(-1, 3))).reshape(tiempo, -1, 3)
+    return matrizsph
+
+def funciondefiltro(Datoscart,selectz,cte):
+    
+    '''
+    # cte es para ajustar las unidades, el codigo da del orden de 10**1
+    '''
+    
+    Datosesf = mcart_msph(Datoscart)
+
+    Datosesf[:,:,0] = Datosesf[:,:,0]/(cte) # Normalizamos
+    Datoscart = msph_mcart(Datosesf)
+    
+    contador = c4.contarpuntos(selectz,Datoscart) # Contamos el numero de puntos que cumplen la condicion en cada tiempo
+    filtro1 = c4.datosfiltrados(selectz,Datoscart) # Filtramos los puntos (nota: la matriz contiene filas de ceros, el numero de filas de ceros varian en el tiempo)
+    filtroesf1 = c4.mcart_msph_filtrado(filtro1)
+
+    return filtro1,filtroesf1
+
+def filtroz(Datoscart,selectz,cte):
+    '''
+    Codigo para hacer la proyeccion (sigue manteniendo 0's  pero interpolar ignora los 0's)
+    '''
+    f1,esf1 = funciondefiltro(Datoscart,selectz,cte)
+    
+    FiltroCart = np.zeros_like(f1)
+    FiltroEsf = np.zeros_like(f1)
+    
+    FiltroEsf[:,:,0] = np.sqrt(f1[:,:,0]**2 + f1[:,:,1]**2)
+    FiltroEsf[:,:,1][esf1[:,:,1] != 0] = np.pi/2 # da igual poner el esf1[1]
+    FiltroEsf[:,:,2] = esf1[:,:,2]
+    
+    FiltroCart = c4.msph_mcart_filtrado(FiltroEsf)
+    FiltroCart[:,:,2] = 0
+    
+    return FiltroCart,FiltroEsf
+

Codigo2D_Interpolacion.py

@@ -0,0 +1,52 @@
+'''
+Fecha: 15/03/25
+Codigo definitivo para la interpolacion de datos en coordenadas esfericas (polares)
+
+--> El metodo de ajuste devuelve datos equiespaciados (por tanto no hay que proyectar constantemente)
+
+--> Se añade un punto al final y principio, esto siempre mejora la precision del metodo porque es una condicion necesaria (alpha = 2pi + alpha)
+
+NOTA: Este codigo ignora los 0's por tanto es compatible con filtroz, porque ademas se encarga de ordenar ascendentemente phi
+'''
+
+import numpy as np
+from scipy.interpolate import PchipInterpolator
+
+def newinterpolate2(matrizsph,n):
+    
+    tiempo = matrizsph.shape[0]
+    MSph = np.empty((tiempo,n,3)) # OBJETIVO
+    
+    theta = np.zeros(n) + np.pi/2
+    MSph[:,:,1] = theta
+    
+    for t in range(tiempo):
+        
+        frame = matrizsph[t]
+        mask = ~np.all(frame == 0, axis=1) # seleccionamos valores distintos de 0
+        frame = frame[mask]
+        
+        if frame.shape[0] < 2:
+            raise ValueError(f"No hay suficientes puntos para interpolar en el tiempo {t}.")
+        
+        x = frame[:,2] # esto es phi
+        y = frame[:,0] # esto es r
+        
+        indices_ordenados = np.argsort(x)
+        
+        x = x[indices_ordenados] # ordenados ascendentemente
+        y = y[indices_ordenados]
+        
+        x = np.concatenate(([x[-1] - 2*np.pi], x, [x[0] + 2*np.pi]))
+        y = np.concatenate(([y[-1]], y, [y[0]]))
+        
+        interp = PchipInterpolator(x, y)
+        
+        #x_new = np.linspace(0, 2*np.pi, n)  # 500 puntos en el intervalo de x EQUIESPACIADO
+        x_new = np.linspace(0, 2*np.pi, n,endpoint = False)  # 500 puntos en el intervalo de x EQUIESPACIADO
+        y_new = interp(x_new)
+        
+        MSph[t,:,2] = x_new # phi
+        MSph[t,:,0] = y_new # r(phi)
+          
+    return MSph

Codigo3D_AjusteHelfrich.py

@@ -0,0 +1,42 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+
+# %% AJUSTE
+
+kb = 1.380649e-23  # Constante de Boltzmann
+T = 300  # Temperatura en Kelvin
+kbT = kb * T  # Factor kb * T
+
+def helfrich_model(x, kappa, b):
+    
+    # b = abs(b)
+    kappa = abs(kappa)
+    b = abs(b)
+
+    result = kbT / (8*b + kappa*x)
+    return result
+
+def aproxHelfrich(x, y):
+    x = np.array(x)
+    y = np.array(y)
+
+    try:
+
+        popt, _ = curve_fit(helfrich_model, x, y, p0=[1.4e-19, 1e-7]) # valores referencia
+        kappa, b = popt
+        
+        # Aseguramos que bend y stretch sean siempre positivos
+        b = abs(b)
+        kappa = abs(kappa)
+
+        # # Establecemos un límite inferior para el valor de 'bend'
+        # bend = max(bend, 1e-20)
+        return kappa, b
+    except RuntimeError:
+        print("Error: La optimización no convergió")
+        return None, None
+
+
+

Codigo3D_Coefs.py

@@ -0,0 +1,323 @@
+'''
+Version del codigo completa para obtener los coeficientes
+PASOS:
+    1. Creamos un grid icosahedrico, y proyectamos sobre el 
+    2. Sobre esta proyeccion que "conserva" el radio de los datos de simulaciones se realiza expansion por arm.es
+    3. De alm se convierte a al
+    FALTA. Se eligen los modos del l = 2 a l = 6, y se realiza ajuste por la funcion (ecuacion 7)
+    
+    NOTA: No es necesario especifiar Lmin, parte del 0 el programa
+    
+'''
+
+# %% LIBRERIAS
+
+import sys,os
+
+from itertools import combinations
+from tqdm import tqdm
+
+import logging
+
+import numpy as np
+
+import orthopoly
+
+import MDAnalysis as mda
+from MDAnalysis.analysis import align
+
+from scipy import spatial
+from scipy.optimize import least_squares
+from scipy.optimize import lsq_linear
+
+
+# %% FUNCIONES PREVIAS
+
+def public_c2s(posicionescart):
+    '''
+    De cartesianas a esfericas (c2s) de forma publica
+    '''
+    x,y,z       = posicionescart.T
+    r,theta,phi = orthopoly.spherical_harmonic.cart2sph(-x,-y,z) #Bug in library
+    return np.asarray([r,theta,phi]).T
+
+def public_s2c(posicionesesf):
+    '''
+    De esfericas a cartesianas (s2c) de forma publica
+    '''
+    r,theta,phi = posicionesesf.T
+    x,y,z       = orthopoly.spherical_harmonic.sph2cart(r,theta,phi)
+    return np.asarray([x,y,z]).T
+
+def msph_mcart(matrizsph):
+    tiempo = matrizsph.shape[0]
+    matrizcart = (public_s2c(matrizsph.reshape(-1, 3))).reshape(tiempo, -1, 3)
+    return matrizcart
+
+def mcart_msph(matrizcart):
+    tiempo = matrizcart.shape[0]
+    matrizsph = (public_c2s(matrizcart.reshape(-1, 3))).reshape(tiempo, -1, 3)
+    return matrizsph
+
+# %% DEFINICION DE LA CLASE
+
+'''
+NOTA: phi [0,2pi]
+    theta [0,pi]
+'''
+
+class SHD:
+
+    def __cart2sph(self,cartPositions):
+        x,y,z       = cartPositions.T
+        r,theta,phi = orthopoly.spherical_harmonic.cart2sph(-x,-y,z) #Bug in library
+        return np.asarray([r,theta,phi]).T
+
+    def __sph2cart(self,sphPositions):
+
+        r,theta,phi = sphPositions.T
+        x,y,z       = orthopoly.spherical_harmonic.sph2cart(r,theta,phi)
+        return np.asarray([x,y,z]).T
+            
+    def __degreeOfTruncationToLM(self,n):
+        A = orthopoly.spherical_harmonic.Tnm(self.__Lmax)
+        l,m = A[0][n],A[1][n]
+        return l,m
+
+    def __init__(self,name="SHD",Lmin=2,Lmax=10,debug=False,
+                 radiusMode    = "geometric",
+                 expansionMode = "abs"):
+
+        self.__availableRadiusModes    = ["expansion","geometric"]
+        self.__availableExpansionModes = ["fluct","abs"]
+
+        self.__DEBUG = debug
+
+        self.__Lmin = Lmin
+        self.__Lmax = Lmax
+        self.__name = name
+
+        self.logger = logging.getLogger("sphAnalysis")
+        #Clean up the logger
+        for handler in self.logger.handlers[:]:
+            self.logger.removeHandler(handler)
+        self.logger.setLevel(logging.DEBUG)
+        #Formatter, ascii time to day, month, year, hour, minute, second
+        formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s',"%d-%b-%Y %H:%M:%S")
+
+        #Set logger to console, with level INFO
+        ch = logging.StreamHandler()
+        ch.setLevel(logging.INFO)
+        ch.setFormatter(formatter)
+
+        #Create other logger. Log to file, with level DEBUG
+        fh = logging.FileHandler(f'{self.__name}.log')
+        fh.setLevel(logging.DEBUG)
+        fh.setFormatter(formatter)
+
+        self.logger.addHandler(ch)
+        self.logger.addHandler(fh)
+
+        ######################################################
+
+        if radiusMode not in self.__availableRadiusModes:
+            self.logger.error("Radius mode {} not available. Available modes are: {}".format(radiusMode,self.__availableRadiusModes))
+            raise ValueError
+        else:
+            self.__radiusMode = radiusMode
+
+        if expansionMode not in self.__availableExpansionModes:
+            self.logger.error("Traj mode {} not available. Available modes are: {}".format(expansionMode,self.__availableExpansionModes))
+            raise ValueError
+        else:
+            self.__expansionMode = expansionMode
+
+        if self.__radiusMode == "expansion" and self.__expansionMode == "fluct":
+            #Not compatible
+            self.logger.error("Radius mode {} and traj mode {} not compatible".format(self.__radiusMode,self.__expansionMode))
+            raise ValueError
+    
+   
+    def __setGrid(self):
+        self.gridSize = self.thetaGrid.shape[0]
+        self.grid = np.asarray([self.thetaGrid,self.phiGrid]).T
+
+        r=[1.0]*self.gridSize
+        self.gridVectors = self.__sph2cart(np.asarray([r,self.thetaGrid,self.phiGrid]).T)
+        self.gridVectorssph = self.__cart2sph(self.gridVectors)
+
+        self.gridKDtree = spatial.KDTree(self.gridVectors)
+     
+    def generateIcosahedralGrid(self,n):
+        
+        print("----------")
+        self.logger.info("Generating icosahedral grid with {} subdivisions ...".format(n))
+
+        self.thetaGrid,self.phiGrid=orthopoly.spherical_harmonic.grid_icosahedral(n)
+        self.__setGrid()
+        
+
+    def distributeTrajPointsAlongGrid(self):
+        """
+        Distribute the trajectory points along the grid.
+        """
+        self.traj = self.traj2 # --> self.traj se deberia obtener desde loadsimulation
+        
+        self.time = int(self.traj.shape[0]) # Numero de "tiempos"
+        
+        self.trajSph = self.__cart2sph(self.traj.reshape(-1,3)).reshape(self.time,-1,3)
+        
+        #Check if the trajectory has been loaded
+        if self.traj is None:
+            self.logger.error("Trying to distribute trajectory points along grid, but no trajectory has been loaded.")
+            raise RuntimeError("No trajectory loaded.")
+
+        #Check if the grid has been generated
+        if self.grid is None:
+            self.logger.error("Trying to distribute trajectory points along grid, but no grid has been generated.")
+            raise RuntimeError("No grid generated.")
+
+        self.logger.info("Distributing trajectory points along grid ...")
+
+        r,theta,phi = self.__cart2sph(self.traj.reshape(-1,3)).T
+
+        #Projects points over sphere
+        r=np.asarray([1.0]*r.shape[0])
+        trajOverSphere = self.__sph2cart(np.asarray([r,theta,phi]).T)
+
+        self.trajBeadsGridPosition = self.gridKDtree.query(trajOverSphere)[1].reshape(self.time,-1)
+
+        self.trajGridSph = np.empty((self.time,self.gridSize,3))
+        for f,[framePos,frameGrid] in enumerate(tqdm(zip(self.trajSph,self.trajBeadsGridPosition),total=self.time)):
+            for i in range(self.gridSize):
+                mask = (frameGrid == i)
+
+                if not np.any(mask):
+                    self.logger.error("Grid point with no data. Try to decrease the number of grid points.")
+                    raise Exception("Grid point with no data")
+
+                r,theta,phi = framePos[mask].T
+                
+                # f es cada tiempo, i es cada fila del grid que se calcula con el valor medio de los particulas
+                # mas cercanas al grid
+
+                self.trajGridSph[f,i] = np.asarray([np.mean(r),self.thetaGrid[i],self.phiGrid[i]])
+        
+        self.trajGrid = self.__sph2cart(self.trajGridSph)
+
+        if self.__DEBUG:
+
+            self.logger.debug("Distributed trajectory points along grid. Writing to file ...")
+
+            with open("gridPos.sp","w") as f:
+                for framePos,frameGrid in zip(self.traj,self.trajBeadsGridPosition):
+                    f.write("#\n")
+                    for p,g in zip(framePos,frameGrid):
+                        f.write("{} {} {} 1.0 {}\n".format(p[0],p[1],p[2],g))
+
+            with open("gridTraj.sp","w") as f:
+                for frame in self.trajGridSph:
+                    f.write("#\n")
+                    for p in frame:
+                        x,y,z = self.__sph2cart(np.asarray([p[0],p[1],p[2]]).T).T
+                        f.write("{} {} {} 1.0 0\n".format(x,y,z))
+                        
+    def sphericalHarmonicExpansion(self,TGS):
+        
+        '''
+        TGS (esf): Son los datos que queremos desarrollar, entonces deberian ser los datos que nos de onlytheta_mitdpi matriz 3D
+        Objetivo --> Calcular los coeficientes
+        '''
+      
+        if self.trajGridSph is None:
+            self.logger.error("Trying to compute spherical harmonic expansion, but trajectory has not been distributed along grid.")
+            raise RuntimeError("No trajectory distributed along grid.")
+
+        self.logger.info("Performing spherical harmonic expansion ...")
+        
+        v1 = TGS.shape[1]
+        v2 = int((self.__Lmax+1)**2)
+        
+        MATRIXY = np.empty((self.time,v1,v2))
+        
+        self.aCoeffComputed      = []
+        
+        # THETA Y PHI NO CAMBIAN EN EL TIEMPO (cambia solo r)
+        
+        th1 = TGS[0,:,1]
+        ph1 = TGS[0,:,2]
+        
+        infoY = orthopoly.spherical_harmonic.sph_har_T_matrix(th1,ph1,self.__Lmax)
+        Y = infoY[0] # Matriz Ylm, como le conocemos theta y phi conocemos esta matriz
+        
+        for i in range(self.time):
+            
+            f,theta,phi = TGS[i].T
+
+            aFluct = lsq_linear(Y,f,max_iter=1000,tol=0.000001,lsq_solver='lsmr') # Sistema de ecuaciones R = Ylm*A donde conocemos R e Ylm
+            self.aCoeffComputed.append(aFluct.x)
+        
+        self.aCoeffComputed = np.asarray(self.aCoeffComputed)
+        self.MatrizY = Y;
+        
+    def valoresmedios(self,M1):
+        
+        # Mensaje de advertencia
+        
+        if self.aCoeffComputed is None:
+            self.logger.error("No se puede calcular valores medios porque no esta definido self.aCoeffComputed")
+            self.logger.error("--> Comprobar función: sphericalHarmonicExpansion")
+            raise RuntimeError("No se hizo nada")
+        self.logger.info("Realizando el calculo de valores medios...")
+        
+        '''
+        Tenemos los armonicos esfericos en formato de matriz de dimensiones (tiempo x (lmax+1)**2)
+        En esta función calculamos:
+            1) <alm**2> : Elevamos la matriz al cuadrado, despues promediamos temporalmente
+            2) <alm> : Promediamos temporalmente (por columnas) # dimensiones de Lmax+1
+        
+        Calcularemos el radio promedio tambien (espacial y temporal)
+        --> M1 es la matriz (en esfericas) de dimensiones tiempoxnumero de particulas o numero puntos grid x3 
+        matriz proyectada para calcular R0
+        '''
+        # CALCULO DE VALORES MEDIOS
+        
+        Vini = self.aCoeffComputed;
+        
+        self.VM_CoefsCuadrado = np.mean(Vini**2,axis = 0) # Valor medio a lo largo de las filas (tiempo)
+        self.VM_Coefs = np.mean(Vini,axis = 0)
+        self.VMCuadrado_Coefs = self.VM_Coefs**2
+        
+        # CALCULO DE LA LIBRERIA (nos permite conocer l y m a partir de n)
+        
+        N = orthopoly.spherical_harmonic.T2nharm(self.__Lmax)
+        
+        lm2n = {}
+        
+        for n in range(N):
+            l,m = self.__degreeOfTruncationToLM(n)
+            lm2n[(l,m)] = n
+        
+        # CALCULO DE <al**2> (ya no depende de m porque sumamos)
+        
+        L   = []
+        a2l = []
+        
+        for l in range(self.__Lmax+1):
+            
+            L.append(l)
+            
+            a2 = 0
+            for m in range(-l,l+1):
+
+                a2       += self.VM_CoefsCuadrado[lm2n[(l,m)]]/(2*l+1)
+
+            a2l.append(a2)
+            
+        va2l = np.asarray(a2l)
+        self.VM_AlCuadrado = va2l # ESTO ES LO QUE NOS INTERESA
+        
+        radial = M1[:,:,0] # dimensiones de tiempo x particulas
+        self.R0 = np.mean(radial)
+        
+

Codigo3D_Ejemplo.py

@@ -0,0 +1,78 @@
+# %% 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) 
+
+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()
+

Codigo3D_FuncionGeneral.py

@@ -0,0 +1,71 @@
+import numpy as np
+
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+from Codigo3D_Coefs import SHD,mcart_msph,msph_mcart
+from colorama import Fore, Style 
+
+import orthopoly
+import Codigo3D_AjusteHelfrich as c3
+
+def calculate_al(Datos,Lmaximo):
+    
+    cdata = SHD(name="SPH",Lmin=0,Lmax=Lmaximo,expansionMode="abs",radiusMode="expansion") # def clase
+    
+    ### CARGAMOS DATOS A LA CLASE
+    
+    factor = 7*10**5
+    d1 = mcart_msph(Datos) # datos en esfericas
+    d1[:,:,0] = d1[:,:,0]/factor # normalizamos el radio al orden de 10**-5 como los datos experimentales
+    d1cart = msph_mcart(d1) # convertimos de nuevo a cartesianas pero con el cambio
+    
+    cdata.traj2 = d1cart # asignamos la trayectoria de la clase
+    
+    print(f"Tenemos: {Style.BRIGHT + Fore.GREEN}{Datos.shape[1]}{Style.RESET_ALL} partículas")
+    print(f"En {Style.BRIGHT + Fore.GREEN}{Datos.shape[0]}{Style.RESET_ALL} instantes temporales")
+    
+    ### GENERAMOS UN GRID
+    
+    cdata.generateIcosahedralGrid(3)
+
+    G2esf = cdata.gridVectorssph  # Puntos del grid en esfericas
+    G2cart = cdata.gridVectors # Puntos del grid en cartesianas (sera exportado)
+
+    print(f"Nuestro grid tiene: {Style.BRIGHT + Fore.GREEN}{G2cart.shape[0]}{Style.RESET_ALL} puntos")
+    
+    #### DISTRIBUIMOS TRAYECTORIA SOBRE GRID
+    
+    print("------")
+    cdata.distributeTrajPointsAlongGrid()
+
+    ProyeccionGrid = cdata.trajGrid # en cartesianas (sera exportado)
+    
+    #### CALCULAMOS LA VARIACION CON RESPECTO A LA SUPERFICIE PROMEDIO
+    
+    Superficie = np.mean(cdata.trajGridSph[:,:,0],axis = 0)
+    
+    VariacionSuper = cdata.trajGridSph
+    VariacionSuper[:,:,0] = VariacionSuper[:,:,0] - Superficie # Por defecto resta en cada fila
+    
+    #### DESARROLLO EN ARMONICOS ESFERICOS
+    
+    '''
+    NOTA: los armonicos esfericos se hacen a partir de la variacion no directamente de la proyeccion del grid
+    '''
+    
+    cdata.sphericalHarmonicExpansion(VariacionSuper) 
+    
+    cdata.trajGridSph = mcart_msph(ProyeccionGrid) # Solucion a un bug, que asigna a la proyeccion variacion super
+    
+    #### CALCULO DE LOS VALORES MEDIOS Y DETERMINACION DE AL
+    
+    cdata.valoresmedios(cdata.trajGridSph) # solo tiene el input para calcular R0. Se calcula R0 del general no de la variacion
+
+    coefs_al2 = cdata.VM_AlCuadrado
+    r0 = cdata.R0
+    
+    lmain = np.linspace(0,Lmaximo,Lmaximo+1)
+    terminol = lmain*(lmain+1)*(lmain-1)*(lmain+2)
+    
+    return G2cart,ProyeccionGrid,terminol,coefs_al2,r0

SPH.log

@@ -0,0 +1,8 @@
+09-Apr-2025 01:30:39 - sphAnalysis - INFO - Generating icosahedral grid with 3 subdivisions ...
+09-Apr-2025 01:30:39 - sphAnalysis - INFO - Distributing trajectory points along grid ...
+09-Apr-2025 01:31:47 - sphAnalysis - INFO - Performing spherical harmonic expansion ...
+09-Apr-2025 01:31:51 - sphAnalysis - INFO - Realizando el calculo de valores medios...
+09-Apr-2025 01:32:35 - sphAnalysis - INFO - Generating icosahedral grid with 3 subdivisions ...
+09-Apr-2025 01:32:36 - sphAnalysis - INFO - Distributing trajectory points along grid ...
+09-Apr-2025 01:33:45 - sphAnalysis - INFO - Performing spherical harmonic expansion ...
+09-Apr-2025 01:33:50 - sphAnalysis - INFO - Realizando el calculo de valores medios...