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}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.
import numpy as np
from scipy.integrate import ode
#%matplotlib
%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
p=-0.1 # momento angular (reescaleado con m)
a = 0.1 # amplitud del campo (reescaleado con e, m y c)
N=5 # número de variables
# 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:
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
#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
r.set_initial_value(u0, x0)
<scipy.integrate._ode.ode at 0xe1ff1bd0>
El valor final de integracion y dividimos el intervalo en trozos iguales a los fines de graficar puntos intermedios.
#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)
# 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
# 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!
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)