Trajectory equations:



In [81]:

    
%matplotlib inline
import matplotlib.pyplot as plt



In [82]:

    
from sympy import *
init_printing()



In [83]:

    
Bx, By, Bz, B = symbols("B_x, B_y, B_z, B")
x, y, z = symbols("x, y, z" )
x_0, y_0, z_0 = symbols("x_0, y_0, z_0")
vx, vy, vz, v = symbols("v_x, v_y, v_z, v")
vx_0, vy_0, vz_0 = symbols("v_x0, v_y0, v_z0")
t = symbols("t")
q, m = symbols("q, m")
c, eps0 = symbols("c, epsilon_0")

The equation of motion: $$ \begin{gather*} m \frac{d^2 \vec{r} }{dt^2} = \frac{q}{c} [ \vec{v} \vec{B} ] \end{gather*} $$

For the case of a uniform magnetic field along the $z$-axis: $$ \vec{B} = B_z = B, \quad B_x = 0, \quad B_y = 0 $$

In Cortesian coordinates:



In [84]:

    
eq_x = Eq( Derivative(x(t), t, 2), q / c / m * Bz * Derivative(y(t),t) )
eq_y = Eq( Derivative(y(t), t, 2), - q / c / m * Bz * Derivative(x(t),t) )
eq_z = Eq( Derivative(z(t), t, 2), 0 )
display( eq_x, eq_y, eq_z )









    





$$\frac{d^{2}}{d t^{2}}  x{\left (t \right )} = \frac{B_{z} q}{c m} \frac{d}{d t} y{\left (t \right )}$$






    





$$\frac{d^{2}}{d t^{2}}  y{\left (t \right )} = - \frac{B_{z} q}{c m} \frac{d}{d t} x{\left (t \right )}$$






    





$$\frac{d^{2}}{d t^{2}}  z{\left (t \right )} = 0$$

Motion is uniform along the $z$-axis:



In [85]:

    
z_eq = dsolve( eq_z, z(t) )
vz_eq = Eq( z_eq.lhs.diff(t), z_eq.rhs.diff(t) )
display( z_eq, vz_eq )









    





$$z{\left (t \right )} = C_{1} + C_{2} t$$






    





$$\frac{d}{d t} z{\left (t \right )} = C_{2}$$

The constants of integration can be found from the initial conditions $z(0) = z_0$ and $v_z(0) = v_{z0}$:



In [86]:

    
c1_c2_system = []
initial_cond_subs = [(t, 0), (z(0), z_0), (diff(z(t),t).subs(t,0), vz_0) ]
c1_c2_system.append( z_eq.subs( initial_cond_subs ) )
c1_c2_system.append( vz_eq.subs( initial_cond_subs ) )

c1, c2 = symbols("C1, C2")
c1_c2 = solve( c1_c2_system, [c1, c2] )
c1_c2









    Out[86]:





$$\left \{ C_{1} : z_{0}, \quad C_{2} : v_{z0}\right \}$$

So that



In [87]:

    
z_sol = z_eq.subs( c1_c2 )
vz_sol = vz_eq.subs( c1_c2 ).subs( [( diff(z(t),t), vz(t) ) ] )
display( z_sol, vz_sol )









    





$$z{\left (t \right )} = t v_{z0} + z_{0}$$






    





$$\operatorname{v_{z}}{\left (t \right )} = v_{z0}$$

For some reason I have not been able to solve the system of differential equations for $x$ and $y$ directly with Sympy's dsolve function:



In [88]:

    
#dsolve( [eq_x, eq_y], [x(t),y(t)] )

It is necessary to resort to the manual solution. The method is to differentiate one of them over time and substitute the other. This will result in oscillator-type second-order equations for $v_y$ and $v_x$. Their solution is known. Integrating one more time, it is possible to obtain laws of motion $x(t)$ and $y(t)$.



In [89]:

    
v_subs = [ (Derivative(x(t),t), vx(t)), (Derivative(y(t),t), vy(t)) ]
eq_vx = eq_x.subs( v_subs )
eq_vy = eq_y.subs( v_subs )
display( eq_vx, eq_vy )

eq_d2t_vx = Eq( diff(eq_vx.lhs,t), diff(eq_vx.rhs,t))
eq_d2t_vx = eq_d2t_vx.subs( [(eq_vy.lhs, eq_vy.rhs)] )
display( eq_d2t_vx )









    





$$\frac{d}{d t} \operatorname{v_{x}}{\left (t \right )} = \frac{B_{z} q}{c m} \operatorname{v_{y}}{\left (t \right )}$$






    





$$\frac{d}{d t} \operatorname{v_{y}}{\left (t \right )} = - \frac{B_{z} q}{c m} \operatorname{v_{x}}{\left (t \right )}$$






    





$$\frac{d^{2}}{d t^{2}}  \operatorname{v_{x}}{\left (t \right )} = - \frac{B_{z}^{2} q^{2}}{c^{2} m^{2}} \operatorname{v_{x}}{\left (t \right )}$$

The solution of the last equation is



In [90]:

    
C1, C2, Omega = symbols( "C1, C2, Omega" )
vx_eq = Eq( vx(t), C1 * cos( Omega * t ) + C2 * sin( Omega * t ))
display( vx_eq )
omega_eq = Eq( Omega, Bz * q / c / m )
display( omega_eq )









    





$$\operatorname{v_{x}}{\left (t \right )} = C_{1} \cos{\left (\Omega t \right )} + C_{2} \sin{\left (\Omega t \right )}$$






    





$$\Omega = \frac{B_{z} q}{c m}$$

where $\Omega$ is a cyclotron frequency.



In [91]:

    
display( vx_eq )

vy_eq = Eq( vy(t), solve( Eq( diff(vx_eq.rhs,t), eq_vx.rhs ), ( vy(t) ) )[0] )
vy_eq = vy_eq.subs( [(Omega*c*m / Bz / q, omega_eq.rhs * c * m / Bz / q)]).simplify()
display( vy_eq )









    





$$\operatorname{v_{x}}{\left (t \right )} = C_{1} \cos{\left (\Omega t \right )} + C_{2} \sin{\left (\Omega t \right )}$$






    





$$\operatorname{v_{y}}{\left (t \right )} = - C_{1} \sin{\left (\Omega t \right )} + C_{2} \cos{\left (\Omega t \right )}$$

For initial conditions $v_x(0) = v_{x0}, v_y(0) = v_{y0}$:



In [92]:

    
initial_cond_subs = [(t,0), (vx(0), vx_0), (vy(0), vy_0) ]
vx0_eq = vx_eq.subs( initial_cond_subs )
vy0_eq = vy_eq.subs( initial_cond_subs )
display( vx0_eq, vy0_eq )

c1_c2 = solve( [vx0_eq, vy0_eq] )
c1_c2_subs = [ ("C1", c1_c2[c1]), ("C2", c1_c2[c2]) ]
vx_eq = vx_eq.subs( c1_c2_subs )
vy_eq = vy_eq.subs( c1_c2_subs )
display( vx_eq, vy_eq )









    





$$v_{x0} = C_{1}$$






    





$$v_{y0} = C_{2}$$






    





$$\operatorname{v_{x}}{\left (t \right )} = v_{x0} \cos{\left (\Omega t \right )} + v_{y0} \sin{\left (\Omega t \right )}$$






    





$$\operatorname{v_{y}}{\left (t \right )} = - v_{x0} \sin{\left (\Omega t \right )} + v_{y0} \cos{\left (\Omega t \right )}$$

These equations can be integrated to obtain the laws of motion:



In [93]:

    
x_eq = vx_eq.subs( vx(t), diff(x(t),t))
x_eq = dsolve( x_eq )
y_eq = vy_eq.subs( vy(t), diff(y(t),t))
y_eq = dsolve( y_eq ).subs( C1, C2 )
display( x_eq, y_eq )









    





$$x{\left (t \right )} = C_{1} + v_{x0} \left(\begin{cases} t & \text{for}\: \Omega = 0 \\\frac{1}{\Omega} \sin{\left (\Omega t \right )} & \text{otherwise} \end{cases}\right) + v_{y0} \left(\begin{cases} 0 & \text{for}\: \Omega = 0 \\- \frac{1}{\Omega} \cos{\left (\Omega t \right )} & \text{otherwise} \end{cases}\right)$$






    





$$y{\left (t \right )} = C_{2} - v_{x0} \left(\begin{cases} 0 & \text{for}\: \Omega = 0 \\- \frac{1}{\Omega} \cos{\left (\Omega t \right )} & \text{otherwise} \end{cases}\right) + v_{y0} \left(\begin{cases} t & \text{for}\: \Omega = 0 \\\frac{1}{\Omega} \sin{\left (\Omega t \right )} & \text{otherwise} \end{cases}\right)$$

For nonzero $\Omega$:



In [94]:

    
x_eq = x_eq.subs( [(Omega, 123)] ).subs( [(123, Omega)] ).subs( [(Rational(1,123), 1/Omega)] )
y_eq = y_eq.subs( [(Omega, 123)] ).subs( [(123, Omega)] ).subs( [(Rational(1,123), 1/Omega)] )
display( x_eq, y_eq )









    





$$x{\left (t \right )} = C_{1} + \frac{v_{x0}}{\Omega} \sin{\left (\Omega t \right )} - \frac{v_{y0}}{\Omega} \cos{\left (\Omega t \right )}$$






    





$$y{\left (t \right )} = C_{2} + \frac{v_{x0}}{\Omega} \cos{\left (\Omega t \right )} + \frac{v_{y0}}{\Omega} \sin{\left (\Omega t \right )}$$

For initial conditions $x(0) = x_0, y(0) = y_0$:



In [95]:

    
initial_cond_subs = [(t,0), (x(0), x_0), (y(0), y_0) ]
x0_eq = x_eq.subs( initial_cond_subs )
y0_eq = y_eq.subs( initial_cond_subs )
display( x0_eq, y0_eq )

c1_c2 = solve( [x0_eq, y0_eq] )
c1_c2_subs = [ ("C1", c1_c2[0][c1]), ("C2", c1_c2[0][c2]) ]
x_eq = x_eq.subs( c1_c2_subs )
y_eq = y_eq.subs( c1_c2_subs )
display( x_eq, y_eq )









    





$$x_{0} = C_{1} - \frac{v_{y0}}{\Omega}$$






    





$$y_{0} = C_{2} + \frac{v_{x0}}{\Omega}$$






    





$$x{\left (t \right )} = \frac{v_{x0}}{\Omega} \sin{\left (\Omega t \right )} - \frac{v_{y0}}{\Omega} \cos{\left (\Omega t \right )} + \frac{1}{\Omega} \left(\Omega x_{0} + v_{y0}\right)$$






    





$$y{\left (t \right )} = \frac{v_{x0}}{\Omega} \cos{\left (\Omega t \right )} + \frac{v_{y0}}{\Omega} \sin{\left (\Omega t \right )} + \frac{1}{\Omega} \left(\Omega y_{0} - v_{x0}\right)$$



In [106]:

    
x_eq = x_eq.simplify()
y_eq = y_eq.simplify()
x_eq = x_eq.expand().collect(Omega)
y_eq = y_eq.expand().collect(Omega)
display( x_eq, y_eq )









    





$$x{\left (t \right )} = x_{0} + \frac{v_{x0}}{\Omega} \sin{\left (\Omega t \right )} - \frac{v_{y0}}{\Omega} \cos{\left (\Omega t \right )} + \frac{v_{y0}}{\Omega}$$






    





$$y{\left (t \right )} = y_{0} + \frac{v_{x0}}{\Omega} \cos{\left (\Omega t \right )} - \frac{v_{x0}}{\Omega} + \frac{v_{y0}}{\Omega} \sin{\left (\Omega t \right )}$$

Finally



In [108]:

    
display( x_eq, y_eq, z_sol )
display( vx_eq, vy_eq, vz_sol )
display( omega_eq )









    





$$x{\left (t \right )} = x_{0} + \frac{v_{x0}}{\Omega} \sin{\left (\Omega t \right )} - \frac{v_{y0}}{\Omega} \cos{\left (\Omega t \right )} + \frac{v_{y0}}{\Omega}$$






    





$$y{\left (t \right )} = y_{0} + \frac{v_{x0}}{\Omega} \cos{\left (\Omega t \right )} - \frac{v_{x0}}{\Omega} + \frac{v_{y0}}{\Omega} \sin{\left (\Omega t \right )}$$






    





$$z{\left (t \right )} = t v_{z0} + z_{0}$$






    





$$\operatorname{v_{x}}{\left (t \right )} = v_{x0} \cos{\left (\Omega t \right )} + v_{y0} \sin{\left (\Omega t \right )}$$






    





$$\operatorname{v_{y}}{\left (t \right )} = - v_{x0} \sin{\left (\Omega t \right )} + v_{y0} \cos{\left (\Omega t \right )}$$






    





$$\operatorname{v_{z}}{\left (t \right )} = v_{z0}$$






    





$$\Omega = \frac{B_{z} q}{c m}$$



In [ ]: