#!/usr/bin/env python # coding: utf-8 # # Movimiento de una partícula cargada en el campo de un dipolo magnético # ## Computa numéricamente las órbitas de un electrón en un campo magnético externo de forma dipolar. # # El Lagrangiano de la partícula libre en presencia de un campo con simetría asimutal es: # # \begin{equation} # \mathcal{L} = -mc^2\sqrt{1-v^2/c^2} + \frac{e}{c}A(r,\theta)r\sin{\theta} v^{\phi}, # \end{equation} # # donde $A$ es la componente del co-vector campo de gauge. # Recordemos que en estas coordenadas tenemos $v^2 = \dot{r}^2 + r^2 \dot{\theta}^2 + r^2\sin^2(\theta) \dot{\phi}^2$. # Es conveniente redefinir $A \to A/c$ $m=1$, $e \to e/m$, ya que las cantidades solo dependen de estos parámetros. # # La independencia del lagrangiano con respecto al tiempo implica que la energía es conservada, en este caso (sin potencial electrostático), viene dada por: # # \begin{equation} # H = \frac{c^2}{ \sqrt{1-v^2/c^2}} = \gamma c^2 , # \end{equation} # lo que implica que también el módulo de la velocidad es conservado. # La invariancia asimutal nos dice que el momento angular asociado es también conservado, # # \begin{equation} # P := \frac{\partial \mathcal{L}}{\partial \dot{\phi}} = \frac{\dot{\phi} r^2 \sin^2{\theta} }{\sqrt{1-v^2/c^2} }+ e A(r,\theta) r \sin{\theta} = \gamma \dot{\phi} r^2 \sin^2{\theta} + e A(r,\theta) r \sin{\theta} # \end{equation} # # En coordenadas esféricas el sistema completo es: # # \begin{eqnarray*} # \frac{\partial \phi}{dt} &=& \frac{h}{\gamma r\sin{\theta}} \;\;\;\;\;\;\; h = [\frac{P}{r\sin{\theta}} - e A] \\ # \frac{\partial \dot{r}}{dt} &=& e \frac{1}{ \gamma}\frac{\partial(rA)}{\partial r} \sin{\theta} \dot{\phi} + \frac{v^2 - \dot{r}^2}{r} \\ # \frac{\partial \dot{\theta}}{dt} &=& [\frac{e}{\gamma} \frac{1}{r}\frac{\partial(\sin{\theta}A)}{\partial\theta} + \sin{\theta}\cos{\theta} \dot{\phi}]\dot{\phi} - 2\frac{\dot{r}}{r}\dot{\theta} # \end{eqnarray*} # # La conservación de $H$ también implica la conservación de $v^2$, usando la expresión para el momento angular, obtenemos, # # \begin{equation} # v^2 = \dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2{\theta} \dot{\phi}^2 \;\;\;\;\; \dot{r}^2 + r^2 \dot{\theta}^2 = v^2 - \frac{h^2}{\gamma^2} . # \end{equation} # # Por lo tanto el movimiento es solo posible cuando $h^2 = [\frac{P}{r\sin{\theta}} - e A] < v^2\gamma^2$ # # Para un multipolo genérico con simetría axial tenemos, # # \begin{equation} # A(r,\theta) = \frac{a_n}{r^{n+1}}P^1_n(\cos{\theta}) # \end{equation} # # # ## Caso dipolar # # En particular tenemos, $P^1_1(\cos{\theta}) = -\sin(\theta)$. # # Definiendo $e a_1 = a_0$, # # \begin{eqnarray*} # \frac{\partial \phi}{dt} &=& \frac{h}{\gamma r\sin{\theta}} \;\;\;\;\;\;\; h = [\frac{P}{r\sin{\theta}} + \frac{a_0}{r^2}\sin{\theta} ] \\ # \frac{\partial \dot{r}}{dt} &=& \frac{1}{\gamma}\frac{a_0}{r^2} \sin^2{\theta} \dot{\phi} + \frac{v^2 - \dot{r}^2}{r} \\ # \frac{\partial \dot{\theta}}{dt} &=& [- \frac{2a_0}{\gamma r^3} + \dot{\phi}]\dot{\phi}\sin{\theta}\cos{\theta} - 2\frac{\dot{r}}{r}\dot{\theta} # \end{eqnarray*} # # Donde # \begin{equation} # v^2 = \frac{\dot{r}^2 + r^2 \dot{\theta}^2 + h^2}{1+h^2}, # \end{equation} # es constante a lo largo del movimiento. # # # In[1]: import numpy as np from scipy.integrate import ode #%matplotlib get_ipython().run_line_magic('matplotlib', 'inline') import matplotlib as mpl import matplotlib.pyplot as plt #import matplotlib.animation as animation from matplotlib import animation, rc from IPython.display import HTML import matplotlib.animation as animation from mpl_toolkits.mplot3d import Axes3D # In[2]: p=-0.1 # momento angular (reescaleado con m) a = 0.1 # amplitud del campo (reescaleado con e, m y c) # In[3]: N=5 # número de variables # In[4]: # Primero las condiciones iniciales U=(rho,rho_t,theta, theta_t,phi) x0 = 0. # tiempo inicial rho0 = 1. #rho_t0 = -0.02 rho_t0 = -0.0 #pert = 0.000000000000000 pert = 0.01 theta0 = np.pi/2.0 theta_t0 = pert phi0 = 0 h0 = (p/rho0/np.sin(theta0) + a*np.sin(theta0)/rho0/rho0) h20 = h0*h0 v2 = (rho_t0*rho_t0 + rho0*rho0*(theta_t0*theta_t0) + h20) / (1.+h20) if (v2 > 0.99): print('velocidad demasiado alta! v2=%f', v2) exit() gamma = np.sqrt(1.-v2) u0 = [rho0, rho_t0, theta0, theta_t0, phi0] # Definimos la funcion del integrador, o sea el sistema de ecuaciones: # In[5]: def f(x, u): sint = np.sin(u[2]) cost = np.cos(u[2]) phi_t = (p/u[0]/sint + a/u[0]/u[0]*sint)/u[0]/sint/gamma return [u[1], a/u[0]/u[0]*sint*sint*phi_t/gamma + (v2 - u[1]*u[1])/u[0],u[3], (-2.*a/u[0]/u[0]/u[0]/gamma + phi_t)*phi_t*sint*cost - 2.*u[1]*u[3]/u[0],phi_t] # Definimos cual integrador usar. Damos distintos tipos, para este problema que en algunos casos es muy sensible, # es recomendable usar la dop853, pero pruebe también la dopri5 # # In[6]: #r = ode(f).set_integrator('zvode', method='bdf') #r = ode(f).set_integrator('dopri5') # Runge Kutta de orden 4/5 r = ode(f).set_integrator('dop853') # Runge Kutta de orden 8 r.rtol = 1.e-6 # Damos los valores iniciales # In[7]: r.set_initial_value(u0, x0) # El valor final de integracion y dividimos el intervalo en trozos iguales a los fines de graficar puntos intermedios. # In[8]: #x1 = 10000 # x1 = 5000. # Esto es solo a los fines de graficar, si da un error: out of bounds poner K mas grande dx = 0.1 K=50001 K_anim = 5000 # puede ser igual a K o KK, pero en el notebook tarda mucho en graficar. #K=100000 #dx=1. #K=10000 KK=0 # para tener solo los puntos justos en las gráficas. U = np.zeros((K,N)) X = np.zeros(K) # In[ ]: # Integramos K veces separadas (para graficar) i = 0 while r.successful() and r.t < x1 and r.y[0] < rho0*12.: r.integrate(r.t+dx) U[i,:] = r.y X[i] = r.t i=i+1 KK=KK+1 # In[10]: # Finalmente graficamos la solucion plt.ioff() ax = plt.subplot(111, projection='3d') ax.grid(True) ax.set_title("general orbits", va='bottom') # cambiar para ajustar mejor el gráfico ax.set_xlim((-1.5, 1.5)) ax.set_ylim((-1.5, 1.5)) ax.set_zlim((-0.5, 0.5)) ax.plot(U[0:KK,0]*np.sin(U[0:KK,4])*np.sin(U[0:KK,2]), U[0:KK,0]*np.cos(U[0:KK,4])*np.sin(U[0:KK,2]), U[0:KK,0]*np.cos(U[0:KK,2]), color='r', linewidth=1) plt.show() # Ahora hacemos una animación del gráfico anterior. Suele demorar, paciencia! # In[11]: fig = plt.figure() #ax1 = fig.add_subplot(111, projection='3d') ax1 = fig.gca(projection='3d') ax1.grid(True) line, = ax1.plot([], [], [], color='r' ,lw=1) time_template = 'time = %.1fs' time_text1 = ax1.text(0.01, 0.01, 0.95, '', transform=ax1.transAxes) ax1.set_xlim((-1.5, 1.5)) ax1.set_ylim((-1.5, 1.5)) ax1.set_zlim((-1.5, 1.5)) def init(): line.set_data([], []) line.set_3d_properties([]) time_text1.set_text('') return line, time_text1 def animate(i): line.set_data(U[0:i,0]*np.sin(U[0:i,4])*np.sin(U[0:i,2]), U[0:i,0]*np.cos(U[0:i,4])*np.sin(U[0:i,2])) line.set_3d_properties(U[0:i,0]*np.cos(U[0:i,2])) time_text1.set_text(time_template % (i*dx)) return line, time_text1 ani = animation.FuncAnimation(fig, animate, K_anim, interval=15, blit=True, init_func=init, repeat=False) #rc('animation', html='html5') HTML(ani.to_html5_video()) #ani.save('3d-dipole.mp4', fps=15) # In[ ]: