commit 43af600d31f2886ffb6c616f6d42ca5c38903055 Author: Gabriel Quispe Date: Wed Apr 9 01:43:45 2025 +0200 Subida inicial del proyecto diff --git a/Codigo2D_AjusteHelfrich.py b/Codigo2D_AjusteHelfrich.py new file mode 100644 index 0000000..ed10291 --- /dev/null +++ b/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 + + diff --git a/Codigo2D_Base.py b/Codigo2D_Base.py new file mode 100644 index 0000000..0ddc7dd --- /dev/null +++ b/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 + + + diff --git a/Codigo2D_CodigoBaseFiltro.py b/Codigo2D_CodigoBaseFiltro.py new file mode 100644 index 0000000..ef036b1 --- /dev/null +++ b/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 diff --git a/Codigo2D_Coefs.py b/Codigo2D_Coefs.py new file mode 100644 index 0000000..ec90070 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Codigo2D_Ejemplo.py b/Codigo2D_Ejemplo.py new file mode 100644 index 0000000..796e4c5 --- /dev/null +++ b/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() \ No newline at end of file diff --git a/Codigo2D_FiltroZ.py b/Codigo2D_FiltroZ.py new file mode 100644 index 0000000..d02ca73 --- /dev/null +++ b/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 + diff --git a/Codigo2D_Interpolacion.py b/Codigo2D_Interpolacion.py new file mode 100644 index 0000000..ff4837f --- /dev/null +++ b/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 diff --git a/Codigo3D_AjusteHelfrich.py b/Codigo3D_AjusteHelfrich.py new file mode 100644 index 0000000..3441abe --- /dev/null +++ b/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 + + + diff --git a/Codigo3D_Coefs.py b/Codigo3D_Coefs.py new file mode 100644 index 0000000..e71504d --- /dev/null +++ b/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) : Elevamos la matriz al cuadrado, despues promediamos temporalmente + 2) : 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 (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) + + diff --git a/Codigo3D_Ejemplo.py b/Codigo3D_Ejemplo.py new file mode 100644 index 0000000..0f5250e --- /dev/null +++ b/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() + diff --git a/Codigo3D_FuncionGeneral.py b/Codigo3D_FuncionGeneral.py new file mode 100644 index 0000000..d408f92 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/SPH.log b/SPH.log new file mode 100644 index 0000000..3cba4fc --- /dev/null +++ b/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... diff --git a/__pycache__/Codigo3D_AjusteHelfrich.cpython-311.pyc b/__pycache__/Codigo3D_AjusteHelfrich.cpython-311.pyc new file mode 100644 index 0000000..8110fb9 Binary files /dev/null and b/__pycache__/Codigo3D_AjusteHelfrich.cpython-311.pyc differ diff --git a/__pycache__/Codigo3D_Coefs.cpython-311.pyc b/__pycache__/Codigo3D_Coefs.cpython-311.pyc new file mode 100644 index 0000000..054247b Binary files /dev/null and b/__pycache__/Codigo3D_Coefs.cpython-311.pyc differ diff --git a/__pycache__/Codigo3D_FuncionGeneral.cpython-311.pyc b/__pycache__/Codigo3D_FuncionGeneral.cpython-311.pyc new file mode 100644 index 0000000..e066b59 Binary files /dev/null and b/__pycache__/Codigo3D_FuncionGeneral.cpython-311.pyc differ diff --git a/matrizcartesianas.npy b/matrizcartesianas.npy new file mode 100644 index 0000000..f9d03b5 Binary files /dev/null and b/matrizcartesianas.npy differ