"""Este script simula el movimiento de una particula cargada en un campo electromagnetico constante, debido a la fuerza de Lorentz. Fisica 2 (Q), Cat. Laura Estrada, 1er cuatrimestre 2018, FCEyN, UBA, Argentina. Autor: Rodrigo Lugones """ ## Llamo a las bibliotecas/librerias que voy a utilizar import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import matplotlib.animation as animation import scipy.integrate as integrate ## Parametros modificables. t_max = 60 # tiempo final (s) de la simulacion particula_v = np.array([1, 0, 1]) particula_x = np.array([0, 1, 0]) #particula_masa = 1 #particula_carga = 1 campo_magnetico = np.array([0, 0, 1]) campo_electrico = np.array([0.25, 0, 0]) # Clases particula y campo. Se podrian complejizar. class Particula(object): """ Clase particula. """ def __init__(self, carga=1, masa=1, inicial_pos=np.array([0, 0, 0]), inicial_vel=np.array([1, 0, 1])): self.q = carga self.m = masa self.x = inicial_pos self.v = inicial_vel def Campo(Fx=0, Fy=0, Fz=1): """ Clase campo. """ return np.array([Fx, Fy, Fz]) def derivs(state, t): """ Esta funcion 'derivs' encuentra numericamente la aceleracion y la velocidad: |-- | dx/dt = v | q _ _ _ | dv/dt = --- . (E + v x B) | m_q |-- """ dim = 3 dxdv = np.zeros_like(state) dxdv[:dim] = state[0+dim:dim+dim] dxdv[dim:] = particula.q / particula.m * (E + np.cross(state[dim:], B)) return dxdv # creo la particula y los campos. particula = Particula(inicial_pos=particula_x, inicial_vel=particula_v) B = Campo(campo_magnetico[0], campo_magnetico[1], campo_magnetico[2]) E = Campo(campo_electrico[0], campo_electrico[1], campo_electrico[2]) # creo un array de tiempos desde 0 hasta t_max, con paso dt dt = 0.05 t = np.arange(0.0, t_max, dt) # array con el estado del sistema a tiempo inicial state = np.concatenate([particula.x, particula.v]) # Encuentra numericamente la posicion. Noten que utiliza la funcion # derivs, el estado inicial y el array de tiempos. El array x es la # posicion no da la posicion x y la velocidad a cada tiempo. A # nosostros, solo nos va a interesar la posicion (por eso, nos # quedamos solo con la columna 0). aux = integrate.odeint(derivs, state, t) x = aux[:, :3] v = aux[:, 3:] # Creo el grafico que se va a ir actualizando fig = plt.figure(figsize=(9,9)) # Este es el subgrafico de arriba #ax = fig.add_subplot(111, projection='3d', autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2), zlim=(0, 10)) ax = fig.add_subplot(3, 1, (1,2), projection='3d', autoscale_on=False, xlim=(min(x[:, 0]), max(x[:, 0])), ylim=(min(x[:, 1]), max(x[:, 1])), zlim=(min(x[:, 2]), max(x[:, 2]))) ax.grid() posicion, = ax.plot([], [], [], 'o-r', lw=2) trayectoria, = ax.plot([], [], [], '-b', alpha=0.85) proyecciones = [ax.plot([], [], ':r', alpha=0.5)[0], ax.plot([], [], ':r', alpha=0.5)[0], ax.plot([], [], ':r', alpha=0.5)[0]] ax.set_xlabel('x') ax.set_ylabel('y') ax.set_zlabel('z') time_template = 'time = %.1fs' time_text = ax.text(0.05, 0.9, 50, '', transform=ax.transAxes) # Estos son los subgraficos de abajo ax2 = fig.add_subplot(4,4,13, autoscale_on=False, xlim=(0, max(t)), ylim=(0, max(0.5*np.sum(v**2, axis=1))*1.1)) ax2.grid() energia, = ax2.plot(t, 0.5*np.sum(v**2, axis=1), '-r', lw=2) ax2.set_xlabel('Tiempo') ax2.set_ylabel('Energía cinética') x_size = (max(x[:, 0]) - min(x[:, 0])) y_size = (max(x[:, 1]) - min(x[:, 1])) z_size = (max(x[:, 2]) - min(x[:, 2])) ax3 = fig.add_subplot(4,4,16, autoscale_on=False, xlim=(min(x[:, 0])-x_size/10, max(x[:, 0])+x_size/10), ylim=(min(x[:, 1])-y_size/10, max(x[:, 1])+y_size/10)) ax3.grid() proyeccion_pos, = ax3.plot([], [], 'o-r', lw=2) proyeccion, = ax3.plot([], [], '-b', lw=2) ax3.set_xlabel('x'); ax3.set_ylabel('y') ax4 = fig.add_subplot(4,3,11, autoscale_on=False, xlim=(min(x[:, 0])-x_size/10, max(x[:, 0])+x_size/10), ylim=(min(x[:, 2])-z_size/10, max(x[:, 2])+z_size/10)) ax4.grid() proyeccion2_pos, = ax4.plot([], [], 'o-r', lw=2) proyeccion2, = ax4.plot([], [], '-b', lw=2) ax4.set_xlabel('x'); ax4.set_ylabel('z') # Creo la animacion def animate(i): posicion.set_data(x[i, 0], x[i, 1]) posicion.set_3d_properties(x[i, 2]) trayectoria.set_data(x[:i, 0], x[:i, 1]) trayectoria.set_3d_properties(x[:i, 2]) #ax.arrow(x[i,0], x[i,1], x[i,2], 1, 1) energia.set_data(t[:i], 0.5*np.sum(v[:i]**2, axis=1)) proyeccion_pos.set_data(x[i, 0], x[i, 1]) proyeccion.set_data(x[:i, 0], x[:i, 1]) proyeccion2_pos.set_data(x[i, 0], x[i, 2]) proyeccion2.set_data(x[:i, 0], x[:i, 2]) #proj_x, = ax.plot(x[:i, 0], x[:i, 1], ':r', zdir='z', zs=0)[0] proyecciones[0] = ax.plot([min(x[:, 0])]*i, x[:i, 1], ':r', alpha=0.5)[0] proyecciones[0].set_3d_properties(x[:i, 2]) proyecciones[1] = ax.plot(x[:i, 0], [max(x[:, 1])]*i, ':r', alpha=0.5)[0] proyecciones[1].set_3d_properties(x[:i, 2]) proyecciones[2] = ax.plot(x[:i, 0], x[:i, 1], ':r', alpha=0.5)[0] proyecciones[2].set_3d_properties([min(x[:, 2])]*i) time_text.set_text(time_template % (i*dt)) return posicion, trayectoria, energia, time_text,\ proyeccion_pos, proyeccion, proyeccion2_pos, proyeccion2, \ proyecciones[0], proyecciones[1], proyecciones[2] # Corro la animacion ani = animation.FuncAnimation(fig, animate, np.arange(1, len(x[:,0])), interval=25, blit=True, repeat=False) #ani.save('double_pendulum.mp4', fps=15) # Con esta linea, puedo guardar un video mp4 plt.show()