1#!/usr/bin/env python3
2"""
3Larmor Gyration Simulation
4
5This script simulates charged particle gyration in a uniform magnetic field
6using the Boris algorithm. Demonstrates circular and helical orbits for
7electrons and protons, and verifies energy conservation.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from mpl_toolkits.mplot3d import Axes3D
16from scipy.constants import e, m_e, m_p, c
17
18
19def boris_push(x, v, q, m, B, dt):
20 """
21 Boris algorithm for particle pushing in electromagnetic fields.
22
23 This is the standard algorithm used in PIC codes for its superior
24 energy conservation properties.
25
26 Parameters:
27 -----------
28 x : array (3,)
29 Position [m]
30 v : array (3,)
31 Velocity [m/s]
32 q : float
33 Charge [C]
34 m : float
35 Mass [kg]
36 B : array (3,)
37 Magnetic field [T]
38 dt : float
39 Timestep [s]
40
41 Returns:
42 --------
43 x_new, v_new : Updated position and velocity
44 """
45 # Half acceleration (no E field in this case)
46 v_minus = v.copy()
47
48 # Rotation
49 t = (q * dt / (2 * m)) * B
50 t_mag2 = np.dot(t, t)
51 s = 2 * t / (1 + t_mag2)
52
53 v_prime = v_minus + np.cross(v_minus, t)
54 v_plus = v_minus + np.cross(v_prime, s)
55
56 # Half acceleration again
57 v_new = v_plus
58
59 # Update position
60 x_new = x + v_new * dt
61
62 return x_new, v_new
63
64
65def simulate_gyration(q, m, B, v0, x0, dt, n_steps):
66 """
67 Simulate particle gyration in uniform magnetic field.
68
69 Parameters:
70 -----------
71 q : float
72 Charge [C]
73 m : float
74 Mass [kg]
75 B : array (3,)
76 Magnetic field [T]
77 v0 : array (3,)
78 Initial velocity [m/s]
79 x0 : array (3,)
80 Initial position [m]
81 dt : float
82 Timestep [s]
83 n_steps : int
84 Number of steps
85
86 Returns:
87 --------
88 trajectory : dict with 'x', 'y', 'z', 'vx', 'vy', 'vz', 't'
89 """
90 # Initialize arrays
91 x = np.zeros((n_steps, 3))
92 v = np.zeros((n_steps, 3))
93 t = np.zeros(n_steps)
94
95 x[0] = x0
96 v[0] = v0
97
98 # Time integration
99 for i in range(n_steps - 1):
100 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, B, dt)
101 t[i+1] = t[i] + dt
102
103 return {
104 'x': x[:, 0], 'y': x[:, 1], 'z': x[:, 2],
105 'vx': v[:, 0], 'vy': v[:, 1], 'vz': v[:, 2],
106 't': t
107 }
108
109
110def theoretical_gyroradius(v_perp, q, m, B):
111 """Calculate theoretical Larmor radius."""
112 omega_c = abs(q) * B / m
113 return v_perp / omega_c
114
115
116def theoretical_gyrofrequency(q, m, B):
117 """Calculate theoretical gyrofrequency."""
118 return abs(q) * B / m
119
120
121def plot_circular_orbit():
122 """Plot circular orbit (no parallel velocity)."""
123
124 # Parameters
125 B0 = 1.0 # Tesla
126 B = np.array([0, 0, B0])
127
128 # Electron
129 q_e = -e
130 m_e_kg = m_e
131 v0_e = np.array([1e6, 0, 0]) # 1000 km/s perpendicular
132 x0 = np.array([0, 0, 0])
133
134 # Proton
135 q_p = e
136 m_p_kg = m_p
137 v0_p = np.array([1e5, 0, 0]) # 100 km/s perpendicular
138
139 # Time parameters
140 omega_ce = abs(q_e) * B0 / m_e_kg
141 T_ce = 2 * np.pi / omega_ce
142 dt = T_ce / 100
143 n_steps = 300
144
145 # Simulate
146 traj_e = simulate_gyration(q_e, m_e_kg, B, v0_e, x0, dt, n_steps)
147 traj_p = simulate_gyration(q_p, m_p_kg, B, v0_p, x0, dt, n_steps)
148
149 # Theoretical values
150 r_L_e = theoretical_gyroradius(np.linalg.norm(v0_e), q_e, m_e_kg, B0)
151 r_L_p = theoretical_gyroradius(np.linalg.norm(v0_p), q_p, m_p_kg, B0)
152
153 # Plot
154 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
155
156 # Electron orbit
157 ax1.plot(traj_e['x'] * 1e3, traj_e['y'] * 1e3, 'b-', linewidth=2, label='Electron')
158 ax1.plot(traj_e['x'][0] * 1e3, traj_e['y'][0] * 1e3, 'go', markersize=10, label='Start')
159 ax1.plot(0, 0, 'r+', markersize=15, markeredgewidth=3, label='Guiding center')
160
161 # Draw circle with theoretical radius
162 theta = np.linspace(0, 2*np.pi, 100)
163 ax1.plot(r_L_e * 1e3 * np.cos(theta), r_L_e * 1e3 * np.sin(theta),
164 'r--', alpha=0.5, linewidth=2, label=f'Theory: r_L = {r_L_e*1e3:.3f} mm')
165
166 ax1.set_xlabel('x [mm]', fontsize=12)
167 ax1.set_ylabel('y [mm]', fontsize=12)
168 ax1.set_title('Electron Larmor Gyration (B = 1 T)', fontsize=13, fontweight='bold')
169 ax1.legend(loc='best', fontsize=10)
170 ax1.grid(True, alpha=0.3)
171 ax1.set_aspect('equal')
172 ax1.arrow(0.8 * r_L_e * 1e3, 0, 0, 0.3 * r_L_e * 1e3, head_width=0.1 * r_L_e * 1e3,
173 head_length=0.05 * r_L_e * 1e3, fc='black', ec='black')
174 ax1.text(0.85 * r_L_e * 1e3, 0.2 * r_L_e * 1e3, 'B', fontsize=14, fontweight='bold')
175
176 # Proton orbit
177 ax2.plot(traj_p['x'] * 1e3, traj_p['y'] * 1e3, 'r-', linewidth=2, label='Proton')
178 ax2.plot(traj_p['x'][0] * 1e3, traj_p['y'][0] * 1e3, 'go', markersize=10, label='Start')
179 ax2.plot(0, 0, 'b+', markersize=15, markeredgewidth=3, label='Guiding center')
180
181 # Draw circle with theoretical radius
182 ax2.plot(r_L_p * 1e3 * np.cos(theta), r_L_p * 1e3 * np.sin(theta),
183 'b--', alpha=0.5, linewidth=2, label=f'Theory: r_L = {r_L_p*1e3:.1f} mm')
184
185 ax2.set_xlabel('x [mm]', fontsize=12)
186 ax2.set_ylabel('y [mm]', fontsize=12)
187 ax2.set_title('Proton Larmor Gyration (B = 1 T)', fontsize=13, fontweight='bold')
188 ax2.legend(loc='best', fontsize=10)
189 ax2.grid(True, alpha=0.3)
190 ax2.set_aspect('equal')
191 ax2.arrow(0.8 * r_L_p * 1e3, 0, 0, 0.3 * r_L_p * 1e3, head_width=2,
192 head_length=1, fc='black', ec='black')
193 ax2.text(0.85 * r_L_p * 1e3, 0.2 * r_L_p * 1e3, 'B', fontsize=14, fontweight='bold')
194
195 plt.tight_layout()
196 plt.savefig('larmor_circular_orbit.png', dpi=300, bbox_inches='tight')
197 plt.show()
198
199
200def plot_helical_orbit():
201 """Plot 3D helical orbit (with parallel velocity)."""
202
203 # Parameters
204 B0 = 1.0 # Tesla
205 B = np.array([0, 0, B0])
206
207 q_e = -e
208 m_e_kg = m_e
209
210 # Initial velocity with both perpendicular and parallel components
211 v_perp = 1e6 # m/s
212 v_para = 5e5 # m/s
213 v0 = np.array([v_perp, 0, v_para])
214 x0 = np.array([0, 0, 0])
215
216 # Time parameters
217 omega_ce = abs(q_e) * B0 / m_e_kg
218 T_ce = 2 * np.pi / omega_ce
219 dt = T_ce / 100
220 n_steps = 500
221
222 # Simulate
223 traj = simulate_gyration(q_e, m_e_kg, B, v0, x0, dt, n_steps)
224
225 # Theoretical values
226 r_L = theoretical_gyroradius(v_perp, q_e, m_e_kg, B0)
227 pitch = v_para * T_ce
228
229 # 3D plot
230 fig = plt.figure(figsize=(12, 9))
231 ax = fig.add_subplot(111, projection='3d')
232
233 ax.plot(traj['x'] * 1e3, traj['y'] * 1e3, traj['z'] * 1e3,
234 'b-', linewidth=2, label='Electron trajectory')
235 ax.plot([traj['x'][0] * 1e3], [traj['y'][0] * 1e3], [traj['z'][0] * 1e3],
236 'go', markersize=10, label='Start')
237 ax.plot([0], [0], traj['z'] * 1e3, 'r--', alpha=0.5, linewidth=2,
238 label='Guiding center line')
239
240 ax.set_xlabel('x [mm]', fontsize=12)
241 ax.set_ylabel('y [mm]', fontsize=12)
242 ax.set_zlabel('z [mm]', fontsize=12)
243 ax.set_title(f'Helical Orbit in Uniform B Field\n'
244 f'v_⊥ = {v_perp/1e3:.0f} km/s, v_∥ = {v_para/1e3:.0f} km/s',
245 fontsize=13, fontweight='bold')
246 ax.legend(loc='best', fontsize=10)
247
248 # Add annotations
249 ax.text2D(0.05, 0.95, f'Larmor radius: {r_L*1e3:.3f} mm\n'
250 f'Pitch: {pitch*1e3:.2f} mm\n'
251 f'Gyroperiod: {T_ce*1e9:.2f} ns',
252 transform=ax.transAxes, fontsize=11,
253 verticalalignment='top',
254 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
255
256 plt.tight_layout()
257 plt.savefig('larmor_helical_orbit.png', dpi=300, bbox_inches='tight')
258 plt.show()
259
260
261def plot_energy_conservation():
262 """Verify energy conservation with Boris algorithm."""
263
264 # Parameters
265 B0 = 1.0 # Tesla
266 B = np.array([0, 0, B0])
267
268 q_e = -e
269 m_e_kg = m_e
270 v0 = np.array([1e6, 5e5, 3e5]) # General velocity
271 x0 = np.array([0, 0, 0])
272
273 # Time parameters
274 omega_ce = abs(q_e) * B0 / m_e_kg
275 T_ce = 2 * np.pi / omega_ce
276
277 # Test different timesteps
278 dt_factors = [0.01, 0.05, 0.1, 0.2]
279 n_periods = 20
280
281 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
282
283 for dt_factor in dt_factors:
284 dt = T_ce * dt_factor
285 n_steps = int(n_periods * T_ce / dt)
286
287 traj = simulate_gyration(q_e, m_e_kg, B, v0, x0, dt, n_steps)
288
289 # Calculate kinetic energy
290 v_squared = traj['vx']**2 + traj['vy']**2 + traj['vz']**2
291 KE = 0.5 * m_e_kg * v_squared
292 KE_relative = (KE - KE[0]) / KE[0]
293
294 # Plot
295 ax1.plot(traj['t'] / T_ce, v_squared / v_squared[0],
296 linewidth=2, label=f'dt = {dt_factor:.2f} T_ce')
297 ax2.semilogy(traj['t'] / T_ce, np.abs(KE_relative),
298 linewidth=2, label=f'dt = {dt_factor:.2f} T_ce')
299
300 ax1.set_xlabel('Time [gyroperiods]', fontsize=12)
301 ax1.set_ylabel('|v|² / |v₀|²', fontsize=12)
302 ax1.set_title('Velocity Magnitude Conservation', fontsize=13, fontweight='bold')
303 ax1.legend(loc='best', fontsize=10)
304 ax1.grid(True, alpha=0.3)
305 ax1.set_ylim([0.999, 1.001])
306
307 ax2.set_xlabel('Time [gyroperiods]', fontsize=12)
308 ax2.set_ylabel('|ΔKE / KE₀|', fontsize=12)
309 ax2.set_title('Relative Energy Error (Boris Algorithm)', fontsize=13, fontweight='bold')
310 ax2.legend(loc='best', fontsize=10)
311 ax2.grid(True, alpha=0.3, which='both')
312
313 plt.tight_layout()
314 plt.savefig('larmor_energy_conservation.png', dpi=300, bbox_inches='tight')
315 plt.show()
316
317
318def plot_gyrofrequency_verification():
319 """Verify gyrofrequency matches theory."""
320
321 # Parameters
322 B0 = 1.0 # Tesla
323 B = np.array([0, 0, B0])
324
325 q_e = -e
326 m_e_kg = m_e
327 v0 = np.array([1e6, 0, 0])
328 x0 = np.array([0, 0, 0])
329
330 # Theoretical frequency
331 omega_ce = abs(q_e) * B0 / m_e_kg
332 f_ce = omega_ce / (2 * np.pi)
333 T_ce = 1 / f_ce
334
335 # Simulate for several periods
336 dt = T_ce / 100
337 n_steps = 500
338 traj = simulate_gyration(q_e, m_e_kg, B, v0, x0, dt, n_steps)
339
340 # Find peaks in x position (gyroperiods)
341 from scipy.signal import find_peaks
342 peaks, _ = find_peaks(traj['x'])
343
344 if len(peaks) > 1:
345 measured_periods = np.diff(traj['t'][peaks])
346 T_measured = np.mean(measured_periods)
347 f_measured = 1 / T_measured
348 else:
349 T_measured = np.nan
350 f_measured = np.nan
351
352 # Plot
353 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 9))
354
355 # Position vs time
356 ax1.plot(traj['t'] * 1e9, traj['x'] * 1e3, 'b-', linewidth=2, label='x(t)')
357 ax1.plot(traj['t'] * 1e9, traj['y'] * 1e3, 'r-', linewidth=2, label='y(t)')
358 if len(peaks) > 0:
359 ax1.plot(traj['t'][peaks] * 1e9, traj['x'][peaks] * 1e3, 'go',
360 markersize=8, label='Peaks')
361
362 ax1.set_xlabel('Time [ns]', fontsize=12)
363 ax1.set_ylabel('Position [mm]', fontsize=12)
364 ax1.set_title('Position vs Time', fontsize=13, fontweight='bold')
365 ax1.legend(loc='best', fontsize=10)
366 ax1.grid(True, alpha=0.3)
367
368 # Velocity phase space
369 ax2.plot(traj['vx'] / 1e6, traj['vy'] / 1e6, 'b-', linewidth=2)
370 ax2.plot(traj['vx'][0] / 1e6, traj['vy'][0] / 1e6, 'go', markersize=10, label='Start')
371
372 ax2.set_xlabel('v_x [10⁶ m/s]', fontsize=12)
373 ax2.set_ylabel('v_y [10⁶ m/s]', fontsize=12)
374 ax2.set_title('Velocity Phase Space', fontsize=13, fontweight='bold')
375 ax2.legend(loc='best', fontsize=10)
376 ax2.grid(True, alpha=0.3)
377 ax2.set_aspect('equal')
378
379 # Add text with frequencies
380 if not np.isnan(f_measured):
381 textstr = f'Theory: f_ce = {f_ce/1e9:.4f} GHz\n' \
382 f'Measured: f = {f_measured/1e9:.4f} GHz\n' \
383 f'Error: {abs(f_ce - f_measured)/f_ce * 100:.4f}%'
384 else:
385 textstr = f'Theory: f_ce = {f_ce/1e9:.4f} GHz'
386
387 ax2.text(0.05, 0.95, textstr, transform=ax2.transAxes, fontsize=11,
388 verticalalignment='top',
389 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
390
391 plt.tight_layout()
392 plt.savefig('larmor_frequency_verification.png', dpi=300, bbox_inches='tight')
393 plt.show()
394
395
396if __name__ == '__main__':
397 print("\n" + "="*80)
398 print("LARMOR GYRATION SIMULATION (Boris Algorithm)")
399 print("="*80 + "\n")
400
401 # Print theoretical values
402 B0 = 1.0
403 print(f"Magnetic field: B = {B0} T\n")
404
405 omega_ce = e * B0 / m_e
406 omega_cp = e * B0 / m_p
407 f_ce = omega_ce / (2 * np.pi)
408 f_cp = omega_cp / (2 * np.pi)
409
410 print("Electron:")
411 print(f" Gyrofrequency: {f_ce/1e9:.4f} GHz")
412 print(f" Gyroperiod: {1/f_ce*1e9:.4f} ns")
413
414 v_perp = 1e6 # m/s
415 r_Le = v_perp / omega_ce
416 print(f" Larmor radius (v_perp = {v_perp/1e6:.1f} Mm/s): {r_Le*1e3:.4f} mm\n")
417
418 print("Proton:")
419 print(f" Gyrofrequency: {f_cp/1e6:.4f} MHz")
420 print(f" Gyroperiod: {1/f_cp*1e6:.4f} μs")
421
422 v_perp = 1e5 # m/s
423 r_Lp = v_perp / omega_cp
424 print(f" Larmor radius (v_perp = {v_perp/1e3:.1f} km/s): {r_Lp*1e3:.2f} mm\n")
425
426 # Generate plots
427 print("Generating plots...")
428 print(" 1. Circular orbits (electron and proton)...")
429 plot_circular_orbit()
430
431 print(" 2. Helical orbit (3D)...")
432 plot_helical_orbit()
433
434 print(" 3. Energy conservation test...")
435 plot_energy_conservation()
436
437 print(" 4. Gyrofrequency verification...")
438 plot_gyrofrequency_verification()
439
440 print("\nDone! Generated 4 figures:")
441 print(" - larmor_circular_orbit.png")
442 print(" - larmor_helical_orbit.png")
443 print(" - larmor_energy_conservation.png")
444 print(" - larmor_frequency_verification.png")