1#!/usr/bin/env python3
2"""
3Comprehensive Particle Orbit Simulator
4
5This script simulates charged particle orbits in various electromagnetic field
6configurations using the Boris algorithm.
7
8Field configurations:
91. Uniform B: gyration
102. Uniform E + B: E×B drift
113. Gradient B: grad-B drift
124. Curved B: curvature drift
135. Mirror B: bounce motion
146. Tokamak-like: banana orbit
15
16Author: Plasma Physics Examples
17License: MIT
18"""
19
20import numpy as np
21import matplotlib.pyplot as plt
22from matplotlib.gridspec import GridSpec
23from mpl_toolkits.mplot3d import Axes3D
24
25# Physical constants
26QE = 1.602176634e-19 # C
27ME = 9.10938356e-31 # kg
28MP = 1.672621898e-27 # kg
29
30class ParticleOrbitSimulator:
31 """Particle orbit simulator using Boris algorithm."""
32
33 def __init__(self, q, m, dt=1e-10):
34 """
35 Initialize simulator.
36
37 Parameters:
38 -----------
39 q : float
40 Particle charge [C]
41 m : float
42 Particle mass [kg]
43 dt : float
44 Time step [s]
45 """
46 self.q = q
47 self.m = m
48 self.dt = dt
49
50 def boris_push(self, r, v, E, B):
51 """
52 Boris algorithm for particle pushing.
53
54 Parameters:
55 -----------
56 r : array (3,)
57 Position [m]
58 v : array (3,)
59 Velocity [m/s]
60 E : array (3,)
61 Electric field [V/m]
62 B : array (3,)
63 Magnetic field [T]
64
65 Returns:
66 --------
67 r_new, v_new : updated position and velocity
68 """
69 # Half acceleration from E field
70 v_minus = v + 0.5 * (self.q / self.m) * E * self.dt
71
72 # Rotation from B field
73 t = 0.5 * (self.q / self.m) * B * self.dt
74 s = 2 * t / (1 + np.dot(t, t))
75
76 v_prime = v_minus + np.cross(v_minus, t)
77 v_plus = v_minus + np.cross(v_prime, s)
78
79 # Half acceleration from E field
80 v_new = v_plus + 0.5 * (self.q / self.m) * E * self.dt
81
82 # Update position
83 r_new = r + v_new * self.dt
84
85 return r_new, v_new
86
87 def simulate(self, r0, v0, E_func, B_func, t_max, save_every=1):
88 """
89 Run simulation.
90
91 Parameters:
92 -----------
93 r0 : array (3,)
94 Initial position [m]
95 v0 : array (3,)
96 Initial velocity [m/s]
97 E_func : callable
98 E_func(r) returns E field at position r
99 B_func : callable
100 B_func(r) returns B field at position r
101 t_max : float
102 Maximum simulation time [s]
103 save_every : int
104 Save every N steps
105
106 Returns:
107 --------
108 t, r, v : time and trajectory arrays
109 """
110 n_steps = int(t_max / self.dt)
111
112 # Storage
113 r_traj = [r0]
114 v_traj = [v0]
115 t_traj = [0]
116
117 r = r0.copy()
118 v = v0.copy()
119
120 for step in range(1, n_steps + 1):
121 E = E_func(r)
122 B = B_func(r)
123
124 r, v = self.boris_push(r, v, E, B)
125
126 if step % save_every == 0:
127 r_traj.append(r.copy())
128 v_traj.append(v.copy())
129 t_traj.append(step * self.dt)
130
131 return np.array(t_traj), np.array(r_traj), np.array(v_traj)
132
133# Field configuration functions
134
135def uniform_B_field(B0):
136 """Uniform magnetic field in z direction."""
137 def B_func(r):
138 return np.array([0, 0, B0])
139 return B_func
140
141def uniform_E_and_B(E0, B0):
142 """Uniform E (x-dir) and B (z-dir) for E×B drift."""
143 def E_func(r):
144 return np.array([E0, 0, 0])
145 def B_func(r):
146 return np.array([0, 0, B0])
147 return E_func, B_func
148
149def gradient_B_field(B0, L):
150 """Gradient B field: B = B0(1 + x/L) ẑ."""
151 def B_func(r):
152 return np.array([0, 0, B0 * (1 + r[0] / L)])
153 return B_func
154
155def curved_B_field(B0, R_c):
156 """Curved B field (toroidal-like)."""
157 def B_func(r):
158 x, y, z = r
159 R = np.sqrt(x**2 + y**2)
160 if R < 1e-10:
161 return np.array([0, 0, B0])
162 # Toroidal field: B_phi = B0 * R0/R
163 B_mag = B0 * R_c / (R + R_c)
164 B_x = -B_mag * y / R
165 B_y = B_mag * x / R
166 return np.array([B_x, B_y, 0])
167 return B_func
168
169def mirror_B_field(B0, L_mirror):
170 """Mirror field: B = B0(1 + (z/L)²) ẑ."""
171 def B_func(r):
172 return np.array([0, 0, B0 * (1 + (r[2] / L_mirror)**2)])
173 return B_func
174
175def tokamak_field(B_tor, B_pol, R0, r_minor):
176 """
177 Simplified tokamak field (toroidal + poloidal).
178
179 Parameters:
180 -----------
181 B_tor : float
182 Toroidal field at major radius [T]
183 B_pol : float
184 Poloidal field strength [T]
185 R0 : float
186 Major radius [m]
187 r_minor : float
188 Minor radius [m]
189 """
190 def B_func(r):
191 x, y, z = r
192 R = np.sqrt(x**2 + y**2)
193
194 if R < 1e-10:
195 return np.array([0, 0, B_tor])
196
197 # Toroidal field (1/R dependence)
198 B_t = B_tor * R0 / R
199 B_phi_x = -B_t * y / R
200 B_phi_y = B_t * x / R
201
202 # Poloidal field (simplified: radial gradient)
203 theta = np.arctan2(z, R - R0)
204 B_p_r = B_pol * np.sin(theta)
205 B_p_z = B_pol * np.cos(theta)
206
207 B_x = B_phi_x + B_p_r * (R - R0) / R if R > 0 else 0
208 B_y = B_phi_y
209 B_z = B_p_z
210
211 return np.array([B_x, B_y, B_z])
212
213 return B_func
214
215def plot_orbit_comparison():
216 """
217 Compare particle orbits in different field configurations.
218 """
219 print("=" * 70)
220 print("Particle Orbit Simulator: Field Configuration Comparison")
221 print("=" * 70)
222
223 # Particle parameters
224 q = QE # Proton
225 m = MP
226 E_kin = 100 # eV
227 v0_mag = np.sqrt(2 * E_kin * QE / m)
228
229 # Field parameters
230 B0 = 1.0 # T
231 E0 = 1000 # V/m
232
233 # Create simulator
234 dt = 1e-9 # 1 ns
235 sim = ParticleOrbitSimulator(q, m, dt)
236
237 # Create figure
238 fig = plt.figure(figsize=(16, 12))
239 gs = GridSpec(3, 3, figure=fig, hspace=0.4, wspace=0.4)
240
241 configurations = []
242
243 # Configuration 1: Uniform B (gyration)
244 print("\n1. Uniform B field (gyration)...")
245 r0 = np.array([0, 0, 0])
246 v0 = np.array([v0_mag, 0, 0])
247 B_func = uniform_B_field(B0)
248 E_func = lambda r: np.array([0, 0, 0])
249
250 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=1e-7, save_every=10)
251
252 # Compute theoretical Larmor radius
253 rho_L = m * v0_mag / (q * B0)
254 omega_c = q * B0 / m
255
256 configurations.append({
257 'name': '1. Uniform B: Gyration',
258 'r': r,
259 't': t,
260 'theory': f'ρL = {rho_L*1e3:.2f} mm, fc = {omega_c/(2*np.pi)/1e6:.1f} MHz',
261 'drift': np.array([0, 0, 0])
262 })
263
264 # Configuration 2: E×B drift
265 print("2. E×B drift...")
266 r0 = np.array([0, 0, 0])
267 v0 = np.array([v0_mag, 0, 0])
268 E_func, B_func = uniform_E_and_B(E0, B0)
269
270 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=1e-7, save_every=10)
271
272 v_ExB = E0 / B0
273
274 configurations.append({
275 'name': '2. E×B Drift',
276 'r': r,
277 't': t,
278 'theory': f'vE×B = {v_ExB/1e3:.1f} km/s',
279 'drift': np.array([0, v_ExB, 0])
280 })
281
282 # Configuration 3: Grad-B drift
283 print("3. Grad-B drift...")
284 L_grad = 1.0 # m
285 r0 = np.array([0, 0, 0])
286 v0 = np.array([v0_mag, 0, 0])
287 B_func = gradient_B_field(B0, L_grad)
288 E_func = lambda r: np.array([0, 0, 0])
289
290 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=5e-7, save_every=10)
291
292 # Theoretical grad-B drift (perpendicular energy)
293 v_perp = v0_mag
294 v_gradB = -m * v_perp**2 / (2 * q * B0**2 * L_grad)
295
296 configurations.append({
297 'name': '3. Grad-B Drift',
298 'r': r,
299 't': t,
300 'theory': f'v∇B ≈ {v_gradB:.1f} m/s',
301 'drift': np.array([0, v_gradB, 0])
302 })
303
304 # Configuration 4: Curved B (curvature drift)
305 print("4. Curvature drift...")
306 R_c = 1.0 # m
307 r0 = np.array([R_c, 0, 0])
308 v0 = np.array([0, 0, v0_mag]) # Parallel velocity
309 B_func = curved_B_field(B0, R_c)
310 E_func = lambda r: np.array([0, 0, 0])
311
312 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=1e-6, save_every=10)
313
314 # Theoretical curvature drift
315 v_parallel = v0_mag
316 v_curv = m * v_parallel**2 / (q * B0 * R_c**2)
317
318 configurations.append({
319 'name': '4. Curvature Drift',
320 'r': r,
321 't': t,
322 'theory': f'vcurv ≈ {v_curv:.1f} m/s',
323 'drift': None
324 })
325
326 # Configuration 5: Mirror field (bounce motion)
327 print("5. Mirror field (bounce)...")
328 L_mirror = 0.5 # m
329 r0 = np.array([0.01, 0, 0]) # Small offset
330 v_parallel = 0.5 * v0_mag
331 v_perp = np.sqrt(v0_mag**2 - v_parallel**2)
332 v0 = np.array([v_perp, 0, v_parallel])
333 B_func = mirror_B_field(B0, L_mirror)
334 E_func = lambda r: np.array([0, 0, 0])
335
336 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=5e-7, save_every=10)
337
338 configurations.append({
339 'name': '5. Mirror Field: Bounce',
340 'r': r,
341 't': t,
342 'theory': f'Mirror ratio = 2.0',
343 'drift': None
344 })
345
346 # Configuration 6: Tokamak (banana orbit)
347 print("6. Tokamak field (banana orbit)...")
348 R0 = 1.0 # m
349 r_minor = 0.3 # m
350 B_tor = 2.0 # T
351 B_pol = 0.2 # T
352
353 r0 = np.array([R0 + 0.1, 0, 0])
354 v0 = np.array([v0_mag * 0.3, v0_mag * 0.3, v0_mag * 0.9])
355 B_func = tokamak_field(B_tor, B_pol, R0, r_minor)
356 E_func = lambda r: np.array([0, 0, 0])
357
358 t, r, v = sim.simulate(r0, v0, E_func, B_func, t_max=2e-6, save_every=20)
359
360 configurations.append({
361 'name': '6. Tokamak: Banana Orbit',
362 'r': r,
363 't': t,
364 'theory': f'R0={R0:.1f}m, a={r_minor:.1f}m',
365 'drift': None
366 })
367
368 # Plot all configurations
369 for idx, config in enumerate(configurations):
370 row = idx // 3
371 col = idx % 3
372
373 if idx < 5:
374 ax = fig.add_subplot(gs[row, col], projection='3d')
375 else:
376 ax = fig.add_subplot(gs[row, col])
377
378 r = config['r']
379
380 if idx < 5:
381 # 3D plot
382 ax.plot(r[:, 0] * 1e3, r[:, 1] * 1e3, r[:, 2] * 1e3,
383 'b-', linewidth=1, alpha=0.7)
384 ax.scatter(r[0, 0] * 1e3, r[0, 1] * 1e3, r[0, 2] * 1e3,
385 c='green', s=50, marker='o', label='Start')
386 ax.scatter(r[-1, 0] * 1e3, r[-1, 1] * 1e3, r[-1, 2] * 1e3,
387 c='red', s=50, marker='x', label='End')
388
389 ax.set_xlabel('x (mm)', fontsize=9)
390 ax.set_ylabel('y (mm)', fontsize=9)
391 ax.set_zlabel('z (mm)', fontsize=9)
392 ax.set_title(config['name'] + '\n' + config['theory'],
393 fontsize=10, fontweight='bold')
394 ax.legend(fontsize=8)
395 else:
396 # 2D projection for tokamak
397 R = np.sqrt(r[:, 0]**2 + r[:, 1]**2)
398 Z = r[:, 2]
399
400 ax.plot(R, Z, 'b-', linewidth=1, alpha=0.7)
401 ax.scatter(R[0], Z[0], c='green', s=50, marker='o', label='Start')
402 ax.scatter(R[-1], Z[-1], c='red', s=50, marker='x', label='End')
403
404 # Draw tokamak boundary
405 theta = np.linspace(0, 2 * np.pi, 100)
406 R_boundary = R0 + r_minor * np.cos(theta)
407 Z_boundary = r_minor * np.sin(theta)
408 ax.plot(R_boundary, Z_boundary, 'k--', linewidth=1, label='Boundary')
409
410 ax.set_xlabel('R (m)', fontsize=9)
411 ax.set_ylabel('Z (m)', fontsize=9)
412 ax.set_title(config['name'] + '\n' + config['theory'],
413 fontsize=10, fontweight='bold')
414 ax.legend(fontsize=8)
415 ax.set_aspect('equal')
416 ax.grid(True, alpha=0.3)
417
418 plt.suptitle('Particle Orbit Simulator: Drift and Bounce Motions',
419 fontsize=16, fontweight='bold', y=0.995)
420
421 plt.savefig('orbit_simulator.png', dpi=150, bbox_inches='tight')
422 print("\nPlot saved as 'orbit_simulator.png'")
423 print("=" * 70)
424
425 plt.show()
426
427if __name__ == "__main__":
428 plot_orbit_comparison()