1#!/usr/bin/env python3
2"""
3Magnetic Mirror Simulation
4
5This script simulates particle motion in a magnetic mirror/bottle configuration.
6Demonstrates trapping, bounce motion, loss cone, and conservation of the
7magnetic moment (adiabatic invariant).
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, m_e
16
17
18def boris_push(x, v, q, m, E, B, dt):
19 """Boris algorithm for particle pushing."""
20 v_minus = v + (q * E / m) * (dt / 2)
21
22 t = (q * dt / (2 * m)) * B
23 t_mag2 = np.dot(t, t)
24 s = 2 * t / (1 + t_mag2)
25
26 v_prime = v_minus + np.cross(v_minus, t)
27 v_plus = v_minus + np.cross(v_prime, s)
28
29 v_new = v_plus + (q * E / m) * (dt / 2)
30 x_new = x + v_new * dt
31
32 return x_new, v_new
33
34
35def magnetic_mirror_field(z, B0, L):
36 """
37 Magnetic mirror field: B_z(z) = B0 * (1 + (z/L)²)
38
39 Parameters:
40 -----------
41 z : float or array
42 Position along z [m]
43 B0 : float
44 Minimum field at center [T]
45 L : float
46 Mirror scale length [m]
47
48 Returns:
49 --------
50 B : array (3,) or (N, 3)
51 Magnetic field [T]
52 """
53 B_z = B0 * (1 + (z / L)**2)
54
55 if np.isscalar(z):
56 return np.array([0, 0, B_z])
57 else:
58 return np.column_stack([np.zeros_like(z), np.zeros_like(z), B_z])
59
60
61def simulate_mirror(q, m, B0, L, v0, x0, dt, n_steps, max_z=None):
62 """
63 Simulate particle in magnetic mirror.
64
65 Parameters:
66 -----------
67 q, m : float
68 Charge and mass
69 B0, L : float
70 Mirror parameters
71 v0, x0 : array
72 Initial velocity and position
73 dt : float
74 Timestep
75 n_steps : int
76 Maximum steps
77 max_z : float
78 Maximum |z| before particle is considered lost
79
80 Returns:
81 --------
82 trajectory : dict
83 """
84 x = np.zeros((n_steps, 3))
85 v = np.zeros((n_steps, 3))
86 t = np.zeros(n_steps)
87 B_field = np.zeros((n_steps, 3))
88 mu = np.zeros(n_steps) # Magnetic moment
89
90 x[0] = x0
91 v[0] = v0
92
93 E = np.array([0, 0, 0])
94 lost = False
95 n_actual = n_steps
96
97 for i in range(n_steps - 1):
98 B = magnetic_mirror_field(x[i, 2], B0, L)
99 B_field[i] = B
100 B_mag = np.linalg.norm(B)
101
102 # Calculate magnetic moment μ = m*v_perp²/(2*B)
103 v_para = np.dot(v[i], B) / B_mag
104 v_perp = np.sqrt(np.dot(v[i], v[i]) - v_para**2)
105 mu[i] = m * v_perp**2 / (2 * B_mag)
106
107 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, E, B, dt)
108 t[i+1] = t[i] + dt
109
110 # Check if particle is lost
111 if max_z is not None and abs(x[i+1, 2]) > max_z:
112 lost = True
113 n_actual = i + 1
114 break
115
116 if not lost:
117 B = magnetic_mirror_field(x[-1, 2], B0, L)
118 B_field[-1] = B
119 B_mag = np.linalg.norm(B)
120 v_para = np.dot(v[-1], B) / B_mag
121 v_perp = np.sqrt(np.dot(v[-1], v[-1]) - v_para**2)
122 mu[-1] = m * v_perp**2 / (2 * B_mag)
123
124 return {
125 'x': x[:n_actual, 0], 'y': x[:n_actual, 1], 'z': x[:n_actual, 2],
126 'vx': v[:n_actual, 0], 'vy': v[:n_actual, 1], 'vz': v[:n_actual, 2],
127 't': t[:n_actual], 'B': B_field[:n_actual],
128 'mu': mu[:n_actual], 'lost': lost
129 }
130
131
132def plot_trapped_vs_lost():
133 """Demonstrate trapped particles vs lost particles."""
134
135 # Parameters
136 B0 = 0.1 # Tesla (minimum field at center)
137 L = 0.5 # meters
138 B_mirror = B0 * (1 + 1)**2 # Field at z = ±L
139
140 mirror_ratio = B_mirror / B0
141
142 q = -e
143 m = m_e
144
145 # Initial position
146 x0 = np.array([0, 0, 0])
147
148 # Test different pitch angles
149 v_total = 1e7 # m/s (total speed)
150
151 # Loss cone angle: sin²(θ_lc) = B0 / B_mirror
152 sin_theta_lc = np.sqrt(B0 / B_mirror)
153 theta_lc = np.arcsin(sin_theta_lc) * 180 / np.pi
154
155 print(f"\nMagnetic Mirror Configuration:")
156 print(f" B0 (center) = {B0} T")
157 print(f" B_mirror (at z=±L) = {B_mirror} T")
158 print(f" Mirror ratio R = {mirror_ratio:.2f}")
159 print(f" Loss cone angle θ_lc = {theta_lc:.2f}°")
160
161 # Pitch angles to test
162 pitch_angles = [20, 40, 50, 60, 70, 80] # degrees
163
164 omega_c = abs(q) * B0 / m
165 T_c = 2 * np.pi / omega_c
166 dt = T_c / 100
167 n_steps = 5000
168 max_z = 2 * L
169
170 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
171
172 colors = plt.cm.rainbow(np.linspace(0, 1, len(pitch_angles)))
173
174 # Storage for loss cone plot
175 v_perp_list = []
176 v_para_list = []
177 lost_list = []
178
179 for pitch_angle, color in zip(pitch_angles, colors):
180 theta_rad = pitch_angle * np.pi / 180
181 v_perp = v_total * np.sin(theta_rad)
182 v_para = v_total * np.cos(theta_rad)
183
184 v0 = np.array([v_perp, 0, v_para])
185
186 traj = simulate_mirror(q, m, B0, L, v0, x0, dt, n_steps, max_z)
187
188 v_perp_list.append(v_perp)
189 v_para_list.append(v_para)
190 lost_list.append(traj['lost'])
191
192 # Plot z vs t
193 label = f"{pitch_angle}° {'(lost)' if traj['lost'] else '(trapped)'}"
194 linestyle = '--' if traj['lost'] else '-'
195 ax = axes[0, 0]
196 ax.plot(traj['t'] * 1e6, traj['z'] * 1e2, linestyle=linestyle,
197 linewidth=2, color=color, label=label)
198
199 # Finish first plot
200 ax = axes[0, 0]
201 ax.axhline(L * 1e2, color='red', linestyle=':', linewidth=2, label='Mirror location')
202 ax.axhline(-L * 1e2, color='red', linestyle=':', linewidth=2)
203 ax.set_xlabel('Time [μs]', fontsize=12)
204 ax.set_ylabel('z Position [cm]', fontsize=12)
205 ax.set_title('Trapped vs Lost Particles', fontsize=13, fontweight='bold')
206 ax.legend(loc='best', fontsize=9)
207 ax.grid(True, alpha=0.3)
208
209 # Plot 2: Loss cone in velocity space
210 ax = axes[0, 1]
211
212 # Plot particles
213 for v_perp, v_para, lost in zip(v_perp_list, v_para_list, lost_list):
214 marker = 'x' if lost else 'o'
215 color = 'red' if lost else 'blue'
216 size = 100
217 ax.scatter(v_para / 1e6, v_perp / 1e6, marker=marker, s=size,
218 color=color, edgecolors='black', linewidth=1.5)
219
220 # Draw loss cone boundary
221 v_para_cone = np.linspace(0, v_total / 1e6, 100)
222 v_perp_cone = v_para_cone * np.tan(theta_lc * np.pi / 180)
223 ax.plot(v_para_cone, v_perp_cone, 'g--', linewidth=3,
224 label=f'Loss cone θ_lc = {theta_lc:.1f}°')
225 ax.plot(v_para_cone, -v_perp_cone, 'g--', linewidth=3)
226
227 # Draw total velocity circle
228 theta = np.linspace(0, 2*np.pi, 100)
229 ax.plot((v_total / 1e6) * np.cos(theta), (v_total / 1e6) * np.sin(theta),
230 'k:', linewidth=2, alpha=0.5, label='|v| = const')
231
232 ax.set_xlabel('v_∥ [10⁶ m/s]', fontsize=12)
233 ax.set_ylabel('v_⊥ [10⁶ m/s]', fontsize=12)
234 ax.set_title('Loss Cone in Velocity Space', fontsize=13, fontweight='bold')
235 ax.legend(loc='best', fontsize=10)
236 ax.grid(True, alpha=0.3)
237 ax.set_aspect('equal')
238
239 # Plot 3: Magnetic field profile
240 ax = axes[1, 0]
241 z_profile = np.linspace(-1.5*L, 1.5*L, 200)
242 B_profile = B0 * (1 + (z_profile / L)**2)
243
244 ax.plot(z_profile * 1e2, B_profile, 'b-', linewidth=3)
245 ax.axvline(L * 1e2, color='red', linestyle='--', linewidth=2, label='Mirror points')
246 ax.axvline(-L * 1e2, color='red', linestyle='--', linewidth=2)
247 ax.axhline(B0, color='green', linestyle=':', linewidth=2, label=f'B0 = {B0} T')
248 ax.axhline(B_mirror, color='orange', linestyle=':', linewidth=2,
249 label=f'B_mirror = {B_mirror} T')
250
251 ax.set_xlabel('z Position [cm]', fontsize=12)
252 ax.set_ylabel('|B| [T]', fontsize=12)
253 ax.set_title('Magnetic Field Profile', fontsize=13, fontweight='bold')
254 ax.legend(loc='best', fontsize=10)
255 ax.grid(True, alpha=0.3)
256
257 # Plot 4: Example trapped particle bounce
258 theta_trapped = 60 # degrees (well above loss cone)
259 v_perp = v_total * np.sin(theta_trapped * np.pi / 180)
260 v_para = v_total * np.cos(theta_trapped * np.pi / 180)
261 v0 = np.array([v_perp, 0, v_para])
262
263 traj = simulate_mirror(q, m, B0, L, v0, x0, dt, n_steps, max_z)
264
265 ax = axes[1, 1]
266 # 2D projection
267 ax.plot(traj['z'] * 1e2, traj['x'] * 1e2, 'b-', linewidth=1, alpha=0.7)
268 ax.plot(traj['z'][0] * 1e2, traj['x'][0] * 1e2, 'go', markersize=10, label='Start')
269 ax.axvline(L * 1e2, color='red', linestyle='--', linewidth=2, label='Mirror')
270 ax.axvline(-L * 1e2, color='red', linestyle='--', linewidth=2)
271
272 ax.set_xlabel('z [cm]', fontsize=12)
273 ax.set_ylabel('x [cm]', fontsize=12)
274 ax.set_title(f'Trapped Particle (θ = {theta_trapped}°)', fontsize=13, fontweight='bold')
275 ax.legend(loc='best', fontsize=10)
276 ax.grid(True, alpha=0.3)
277 ax.set_aspect('equal')
278
279 plt.tight_layout()
280 plt.savefig('magnetic_mirror_trapped_lost.png', dpi=300, bbox_inches='tight')
281 plt.show()
282
283
284def plot_adiabatic_invariant():
285 """Verify conservation of magnetic moment (adiabatic invariant)."""
286
287 # Parameters
288 B0 = 0.1 # Tesla
289 L = 0.5 # meters
290
291 q = -e
292 m = m_e
293
294 # Trapped particle
295 v_total = 1e7 # m/s
296 theta = 60 * np.pi / 180 # Pitch angle
297 v_perp = v_total * np.sin(theta)
298 v_para = v_total * np.cos(theta)
299
300 v0 = np.array([v_perp, 0, v_para])
301 x0 = np.array([0, 0, 0])
302
303 omega_c = abs(q) * B0 / m
304 T_c = 2 * np.pi / omega_c
305 dt = T_c / 100
306 n_steps = 5000
307
308 traj = simulate_mirror(q, m, B0, L, v0, x0, dt, n_steps)
309
310 # Calculate B magnitude along trajectory
311 B_mag = np.linalg.norm(traj['B'], axis=1)
312
313 # Calculate v_perp and v_para
314 v_perp_traj = np.zeros(len(traj['t']))
315 v_para_traj = np.zeros(len(traj['t']))
316
317 for i in range(len(traj['t'])):
318 v_vec = np.array([traj['vx'][i], traj['vy'][i], traj['vz'][i]])
319 B_vec = traj['B'][i]
320 B_norm = np.linalg.norm(B_vec)
321
322 v_para_traj[i] = np.dot(v_vec, B_vec) / B_norm
323 v_perp_traj[i] = np.sqrt(np.dot(v_vec, v_vec) - v_para_traj[i]**2)
324
325 # Create figure
326 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
327
328 # Plot 1: Bounce trajectory
329 ax = axes[0, 0]
330 ax.plot(traj['t'] * 1e6, traj['z'] * 1e2, 'b-', linewidth=2)
331 ax.axhline(L * 1e2, color='red', linestyle='--', linewidth=2, label='Mirror points')
332 ax.axhline(-L * 1e2, color='red', linestyle='--', linewidth=2)
333
334 ax.set_xlabel('Time [μs]', fontsize=12)
335 ax.set_ylabel('z Position [cm]', fontsize=12)
336 ax.set_title('Bounce Motion', fontsize=13, fontweight='bold')
337 ax.legend(loc='best', fontsize=10)
338 ax.grid(True, alpha=0.3)
339
340 # Plot 2: Magnetic moment conservation
341 ax = axes[0, 1]
342 mu_normalized = traj['mu'] / traj['mu'][0]
343 ax.plot(traj['t'] * 1e6, mu_normalized, 'b-', linewidth=2)
344 ax.axhline(1.0, color='red', linestyle='--', linewidth=2, label='Perfect conservation')
345
346 ax.set_xlabel('Time [μs]', fontsize=12)
347 ax.set_ylabel('μ / μ₀', fontsize=12)
348 ax.set_title('Magnetic Moment Conservation', fontsize=13, fontweight='bold')
349 ax.legend(loc='best', fontsize=10)
350 ax.grid(True, alpha=0.3)
351 ax.set_ylim([0.995, 1.005])
352
353 # Print conservation statistics
354 mu_error = np.abs(traj['mu'] - traj['mu'][0]) / traj['mu'][0]
355 print(f"\nMagnetic Moment Conservation:")
356 print(f" Maximum error: {np.max(mu_error) * 100:.4f}%")
357 print(f" Mean error: {np.mean(mu_error) * 100:.4f}%")
358
359 # Plot 3: v_perp vs B (should follow v_perp² ∝ B)
360 ax = axes[1, 0]
361 ax.plot(B_mag, v_perp_traj / 1e6, 'b.', markersize=2, alpha=0.5, label='Simulation')
362
363 # Theoretical: v_perp² = 2μB/m, so v_perp = sqrt(2μB/m)
364 B_theory = np.linspace(B_mag.min(), B_mag.max(), 100)
365 mu_avg = np.mean(traj['mu'])
366 v_perp_theory = np.sqrt(2 * mu_avg * B_theory / m)
367 ax.plot(B_theory, v_perp_theory / 1e6, 'r-', linewidth=3, label='Theory: v_⊥ ∝ √B')
368
369 ax.set_xlabel('|B| [T]', fontsize=12)
370 ax.set_ylabel('v_⊥ [10⁶ m/s]', fontsize=12)
371 ax.set_title('v_⊥ vs B (Adiabatic Relation)', fontsize=13, fontweight='bold')
372 ax.legend(loc='best', fontsize=10)
373 ax.grid(True, alpha=0.3)
374
375 # Plot 4: v_para vs z (reverses at mirror points)
376 ax = axes[1, 1]
377 ax.plot(traj['z'] * 1e2, v_para_traj / 1e6, 'b-', linewidth=2)
378 ax.axvline(L * 1e2, color='red', linestyle='--', linewidth=2, label='Mirror points')
379 ax.axvline(-L * 1e2, color='red', linestyle='--', linewidth=2)
380 ax.axhline(0, color='gray', linestyle='-', linewidth=1)
381
382 ax.set_xlabel('z Position [cm]', fontsize=12)
383 ax.set_ylabel('v_∥ [10⁶ m/s]', fontsize=12)
384 ax.set_title('Parallel Velocity (Reverses at Mirror)', fontsize=13, fontweight='bold')
385 ax.legend(loc='best', fontsize=10)
386 ax.grid(True, alpha=0.3)
387
388 plt.tight_layout()
389 plt.savefig('magnetic_mirror_adiabatic_invariant.png', dpi=300, bbox_inches='tight')
390 plt.show()
391
392
393if __name__ == '__main__':
394 print("\n" + "="*80)
395 print("MAGNETIC MIRROR SIMULATION")
396 print("="*80)
397
398 print("\nPart 1: Trapped vs Lost Particles")
399 print("-" * 80)
400 plot_trapped_vs_lost()
401
402 print("\nPart 2: Adiabatic Invariant (Magnetic Moment Conservation)")
403 print("-" * 80)
404 plot_adiabatic_invariant()
405
406 print("\nDone! Generated 2 figures:")
407 print(" - magnetic_mirror_trapped_lost.png")
408 print(" - magnetic_mirror_adiabatic_invariant.png")