1#!/usr/bin/env python3
2"""
3Gradient and Curvature Drift Simulation
4
5This script demonstrates grad-B drift and curvature drift in non-uniform
6magnetic fields. Shows that these drifts are charge-dependent (opposite
7directions for electrons and ions).
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, m_p
16
17
18def boris_push(x, v, q, m, E, B, dt):
19 """
20 Boris algorithm for particle pushing.
21
22 Parameters:
23 -----------
24 x, v : array (3,)
25 Position [m] and velocity [m/s]
26 q, m : float
27 Charge [C] and mass [kg]
28 E, B : array (3,)
29 Electric [V/m] and magnetic [T] fields
30 dt : float
31 Timestep [s]
32
33 Returns:
34 --------
35 x_new, v_new : Updated position and velocity
36 """
37 # Half acceleration
38 v_minus = v + (q * E / m) * (dt / 2)
39
40 # Rotation
41 t = (q * dt / (2 * m)) * B
42 t_mag2 = np.dot(t, t)
43 s = 2 * t / (1 + t_mag2)
44
45 v_prime = v_minus + np.cross(v_minus, t)
46 v_plus = v_minus + np.cross(v_prime, s)
47
48 # Half acceleration
49 v_new = v_plus + (q * E / m) * (dt / 2)
50
51 # Update position
52 x_new = x + v_new * dt
53
54 return x_new, v_new
55
56
57def magnetic_field_gradient(x, B0, L):
58 """
59 Linearly increasing magnetic field in z direction.
60 B = B0 * (1 + x/L) * ẑ
61
62 Parameters:
63 -----------
64 x : array (3,)
65 Position [m]
66 B0 : float
67 Reference field [T]
68 L : float
69 Gradient scale length [m]
70
71 Returns:
72 --------
73 B : array (3,)
74 Magnetic field [T]
75 """
76 B_magnitude = B0 * (1 + x[0] / L)
77 return np.array([0, 0, B_magnitude])
78
79
80def magnetic_field_dipole(x, B0, R0):
81 """
82 Simplified dipole-like field (2D approximation).
83
84 In cylindrical coordinates centered on z-axis:
85 B_r ≈ -B0 * (R0/r)³ * sin(θ)
86 B_z ≈ B0 * (R0/r)³ * 2*cos(θ)
87
88 Simplified for |z| << R:
89 B_x ≈ -3 * B0 * R0² * x * z / r⁵
90 B_y ≈ -3 * B0 * R0² * y * z / r⁵
91 B_z ≈ B0 * R0² * (2*z² - x² - y²) / r⁵
92
93 Parameters:
94 -----------
95 x : array (3,)
96 Position [m]
97 B0 : float
98 Reference field at equator [T]
99 R0 : float
100 Reference radius [m]
101
102 Returns:
103 --------
104 B : array (3,)
105 Magnetic field [T]
106 """
107 r2 = x[0]**2 + x[1]**2 + x[2]**2
108 r2 = max(r2, (0.01 * R0)**2) # Avoid singularity
109 r5 = r2**2.5
110
111 R02 = R0**2
112
113 Bx = -3 * B0 * R02 * x[0] * x[2] / r5
114 By = -3 * B0 * R02 * x[1] * x[2] / r5
115 Bz = B0 * R02 * (2*x[2]**2 - x[0]**2 - x[1]**2) / r5
116
117 return np.array([Bx, By, Bz])
118
119
120def simulate_grad_b_drift(q, m, B0, L, v0, x0, dt, n_steps):
121 """Simulate particle in gradient-B field."""
122
123 x = np.zeros((n_steps, 3))
124 v = np.zeros((n_steps, 3))
125 t = np.zeros(n_steps)
126 B_field = np.zeros((n_steps, 3))
127
128 x[0] = x0
129 v[0] = v0
130
131 E = np.array([0, 0, 0]) # No electric field
132
133 for i in range(n_steps - 1):
134 B = magnetic_field_gradient(x[i], B0, L)
135 B_field[i] = B
136 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, E, B, dt)
137 t[i+1] = t[i] + dt
138
139 B_field[-1] = magnetic_field_gradient(x[-1], B0, L)
140
141 return {
142 'x': x[:, 0], 'y': x[:, 1], 'z': x[:, 2],
143 'vx': v[:, 0], 'vy': v[:, 1], 'vz': v[:, 2],
144 't': t, 'B': B_field
145 }
146
147
148def simulate_dipole_drift(q, m, B0, R0, v0, x0, dt, n_steps):
149 """Simulate particle in dipole-like field."""
150
151 x = np.zeros((n_steps, 3))
152 v = np.zeros((n_steps, 3))
153 t = np.zeros(n_steps)
154
155 x[0] = x0
156 v[0] = v0
157
158 E = np.array([0, 0, 0])
159
160 for i in range(n_steps - 1):
161 B = magnetic_field_dipole(x[i], B0, R0)
162 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, E, B, dt)
163 t[i+1] = t[i] + dt
164
165 return {
166 'x': x[:, 0], 'y': x[:, 1], 'z': x[:, 2],
167 'vx': v[:, 0], 'vy': v[:, 1], 'vz': v[:, 2],
168 't': t
169 }
170
171
172def plot_grad_b_drift():
173 """Plot grad-B drift showing opposite drift for electron and ion."""
174
175 # Parameters
176 B0 = 1.0 # Tesla
177 L = 0.1 # meters (gradient scale length)
178
179 # Particles
180 q_e = -e
181 q_p = e
182 m_e_kg = m_e
183 m_p_kg = m_p
184
185 # Initial conditions (perpendicular velocity)
186 v_perp = 1e6 # m/s
187 v0_e = np.array([0, v_perp, 0])
188 v0_p = np.array([0, v_perp, 0])
189 x0 = np.array([0, 0, 0])
190
191 # Time parameters
192 omega_ce = abs(q_e) * B0 / m_e_kg
193 T_ce = 2 * np.pi / omega_ce
194 dt = T_ce / 100
195 n_steps = 2000
196
197 # Simulate
198 traj_e = simulate_grad_b_drift(q_e, m_e_kg, B0, L, v0_e, x0, dt, n_steps)
199 traj_p = simulate_grad_b_drift(q_p, m_p_kg, B0, L, v0_p, x0, dt, n_steps)
200
201 # Theoretical drift velocity (grad-B drift)
202 # v_∇B = (m * v_perp²) / (2 * q * B²) * (B × ∇B) / B
203 # For B = B0(1 + x/L)ẑ, ∇B = (B0/L)x̂
204 # B × ∇B = B0²(1+x/L)/L * (ẑ × x̂) = B0²(1+x/L)/L * ŷ
205 # At x=0: v_∇B ≈ (m * v_perp²) / (2 * q * B0 * L) * ŷ
206
207 v_drift_e_theory = (m_e_kg * v_perp**2) / (2 * abs(q_e) * B0**2 * L)
208 v_drift_p_theory = (m_p_kg * v_perp**2) / (2 * abs(q_p) * B0**2 * L)
209
210 # Calculate drift from simulation
211 t_start = int(0.3 * n_steps)
212 drift_e_y = np.polyfit(traj_e['t'][t_start:], traj_e['y'][t_start:], 1)[0]
213 drift_p_y = np.polyfit(traj_p['t'][t_start:], traj_p['y'][t_start:], 1)[0]
214
215 print(f"\nGrad-B Drift:")
216 print(f" B = B0(1 + x/L)ẑ, B0 = {B0} T, L = {L} m")
217 print(f" v_perp = {v_perp/1e6:.1f} Mm/s")
218 print(f"\nTheoretical drift velocities:")
219 print(f" Electron: {v_drift_e_theory:.2f} m/s (negative y)")
220 print(f" Proton: {v_drift_p_theory:.2f} m/s (positive y)")
221 print(f"\nSimulated drift velocities:")
222 print(f" Electron: {drift_e_y:.2f} m/s")
223 print(f" Proton: {drift_p_y:.2f} m/s")
224
225 # Create figure
226 fig, axes = plt.subplots(2, 2, figsize=(14, 12))
227
228 # Plot 1: Trajectories in x-y plane
229 ax1 = axes[0, 0]
230 ax1.plot(traj_e['x'] * 1e3, traj_e['y'] * 1e3, 'b-',
231 linewidth=1.5, alpha=0.7, label='Electron')
232 ax1.plot(traj_p['x'] * 1e3, traj_p['y'] * 1e3, 'r-',
233 linewidth=1.5, alpha=0.7, label='Proton')
234 ax1.plot(0, 0, 'go', markersize=10, label='Start')
235
236 # Show gradient direction
237 ax1.arrow(20, -20, 10, 0, head_width=3, head_length=2,
238 fc='purple', ec='purple', linewidth=2)
239 ax1.text(35, -20, '∇B', fontsize=12, fontweight='bold', color='purple')
240
241 ax1.set_xlabel('x [mm]', fontsize=12)
242 ax1.set_ylabel('y [mm]', fontsize=12)
243 ax1.set_title('Grad-B Drift Trajectories', fontsize=13, fontweight='bold')
244 ax1.legend(loc='best', fontsize=10)
245 ax1.grid(True, alpha=0.3)
246 ax1.set_aspect('equal')
247
248 # Plot 2: y position vs time
249 ax2 = axes[0, 1]
250 ax2.plot(traj_e['t'] * 1e6, traj_e['y'] * 1e3, 'b-', linewidth=2, label='Electron (drifts down)')
251 ax2.plot(traj_p['t'] * 1e6, traj_p['y'] * 1e3, 'r-', linewidth=2, label='Proton (drifts up)')
252
253 ax2.set_xlabel('Time [μs]', fontsize=12)
254 ax2.set_ylabel('y Position [mm]', fontsize=12)
255 ax2.set_title('Opposite Drift Directions', fontsize=13, fontweight='bold')
256 ax2.legend(loc='best', fontsize=10)
257 ax2.grid(True, alpha=0.3)
258
259 # Plot 3: Magnetic field strength along trajectory
260 ax3 = axes[1, 0]
261 B_mag_e = np.linalg.norm(traj_e['B'], axis=1)
262 ax3.plot(traj_e['x'] * 1e3, B_mag_e, 'b-', linewidth=2, label='|B| along trajectory')
263 ax3.axhline(B0, color='gray', linestyle='--', linewidth=2, label=f'B0 = {B0} T')
264
265 ax3.set_xlabel('x Position [mm]', fontsize=12)
266 ax3.set_ylabel('|B| [T]', fontsize=12)
267 ax3.set_title('Magnetic Field Strength', fontsize=13, fontweight='bold')
268 ax3.legend(loc='best', fontsize=10)
269 ax3.grid(True, alpha=0.3)
270
271 # Plot 4: Drift velocity comparison
272 ax4 = axes[1, 1]
273
274 # Calculate instantaneous drift by averaging over gyro-orbits
275 window = 50
276 from scipy.ndimage import uniform_filter1d
277 smooth_y_e = uniform_filter1d(traj_e['y'], window)
278 smooth_y_p = uniform_filter1d(traj_p['y'], window)
279 v_drift_e = np.gradient(smooth_y_e, traj_e['t'])
280 v_drift_p = np.gradient(smooth_y_p, traj_p['t'])
281
282 ax4.plot(traj_e['t'] * 1e6, v_drift_e, 'b-', linewidth=2, alpha=0.7, label='Electron')
283 ax4.plot(traj_p['t'] * 1e6, v_drift_p, 'r-', linewidth=2, alpha=0.7, label='Proton')
284 ax4.axhline(-v_drift_e_theory, color='blue', linestyle='--',
285 linewidth=2, label=f'Theory (e⁻): {-v_drift_e_theory:.2f} m/s')
286 ax4.axhline(v_drift_p_theory, color='red', linestyle='--',
287 linewidth=2, label=f'Theory (p⁺): {v_drift_p_theory:.2f} m/s')
288
289 ax4.set_xlabel('Time [μs]', fontsize=12)
290 ax4.set_ylabel('Drift Velocity v_y [m/s]', fontsize=12)
291 ax4.set_title('Drift Velocity vs Theory', fontsize=13, fontweight='bold')
292 ax4.legend(loc='best', fontsize=9)
293 ax4.grid(True, alpha=0.3)
294
295 plt.tight_layout()
296 plt.savefig('grad_b_drift.png', dpi=300, bbox_inches='tight')
297 plt.show()
298
299
300def plot_curvature_drift():
301 """Plot curvature drift in dipole-like field."""
302
303 # Parameters
304 B0 = 1.0 # Tesla at equator
305 R0 = 0.5 # Reference radius [m]
306
307 # Particles
308 q_e = -e
309 q_p = e
310 m_e_kg = m_e
311 m_p_kg = m_p
312
313 # Initial conditions (in x-z plane, with parallel velocity)
314 v_para = 5e5 # m/s
315 v_perp = 1e6 # m/s
316 x0_e = np.array([R0, 0, 0])
317 x0_p = np.array([R0, 0, 0])
318 v0_e = np.array([0, v_perp, v_para])
319 v0_p = np.array([0, v_perp, v_para])
320
321 # Time parameters
322 omega_ce = abs(q_e) * B0 / m_e_kg
323 T_ce = 2 * np.pi / omega_ce
324 dt = T_ce / 50
325 n_steps = 1500
326
327 # Simulate
328 traj_e = simulate_dipole_drift(q_e, m_e_kg, B0, R0, v0_e, x0_e, dt, n_steps)
329 traj_p = simulate_dipole_drift(q_p, m_p_kg, B0, R0, v0_p, x0_p, dt, n_steps)
330
331 # Create figure
332 fig = plt.figure(figsize=(14, 10))
333
334 # 3D trajectory
335 ax1 = fig.add_subplot(2, 2, 1, projection='3d')
336 ax1.plot(traj_e['x'] * 1e2, traj_e['y'] * 1e2, traj_e['z'] * 1e2,
337 'b-', linewidth=1.5, label='Electron')
338 ax1.plot([traj_e['x'][0] * 1e2], [traj_e['y'][0] * 1e2], [traj_e['z'][0] * 1e2],
339 'go', markersize=10, label='Start')
340
341 ax1.set_xlabel('x [cm]', fontsize=11)
342 ax1.set_ylabel('y [cm]', fontsize=11)
343 ax1.set_zlabel('z [cm]', fontsize=11)
344 ax1.set_title('3D Trajectory in Dipole Field', fontsize=13, fontweight='bold')
345 ax1.legend(loc='best', fontsize=9)
346
347 # x-y projection
348 ax2 = fig.add_subplot(2, 2, 2)
349 ax2.plot(traj_e['x'] * 1e2, traj_e['y'] * 1e2, 'b-',
350 linewidth=1.5, alpha=0.7, label='Electron')
351 ax2.plot(traj_p['x'] * 1e2, traj_p['y'] * 1e2, 'r-',
352 linewidth=1.5, alpha=0.7, label='Proton')
353 ax2.plot(0, 0, 'k*', markersize=15, label='Center')
354 ax2.plot(traj_e['x'][0] * 1e2, traj_e['y'][0] * 1e2,
355 'go', markersize=10, label='Start')
356
357 ax2.set_xlabel('x [cm]', fontsize=12)
358 ax2.set_ylabel('y [cm]', fontsize=12)
359 ax2.set_title('x-y Projection (Curvature Drift)', fontsize=13, fontweight='bold')
360 ax2.legend(loc='best', fontsize=10)
361 ax2.grid(True, alpha=0.3)
362 ax2.set_aspect('equal')
363
364 # x-z projection
365 ax3 = fig.add_subplot(2, 2, 3)
366 ax3.plot(traj_e['x'] * 1e2, traj_e['z'] * 1e2, 'b-',
367 linewidth=1.5, alpha=0.7, label='Electron')
368 ax3.plot(traj_p['x'] * 1e2, traj_p['z'] * 1e2, 'r-',
369 linewidth=1.5, alpha=0.7, label='Proton')
370 ax3.plot(0, 0, 'k*', markersize=15, label='Center')
371
372 ax3.set_xlabel('x [cm]', fontsize=12)
373 ax3.set_ylabel('z [cm]', fontsize=12)
374 ax3.set_title('x-z Projection', fontsize=13, fontweight='bold')
375 ax3.legend(loc='best', fontsize=10)
376 ax3.grid(True, alpha=0.3)
377 ax3.set_aspect('equal')
378
379 # Radial distance vs time
380 ax4 = fig.add_subplot(2, 2, 4)
381 r_e = np.sqrt(traj_e['x']**2 + traj_e['y']**2 + traj_e['z']**2)
382 r_p = np.sqrt(traj_p['x']**2 + traj_p['y']**2 + traj_p['z']**2)
383
384 ax4.plot(traj_e['t'] * 1e6, r_e * 1e2, 'b-', linewidth=2, label='Electron')
385 ax4.plot(traj_p['t'] * 1e6, r_p * 1e2, 'r-', linewidth=2, label='Proton')
386
387 ax4.set_xlabel('Time [μs]', fontsize=12)
388 ax4.set_ylabel('Radial Distance [cm]', fontsize=12)
389 ax4.set_title('Distance from Center', fontsize=13, fontweight='bold')
390 ax4.legend(loc='best', fontsize=10)
391 ax4.grid(True, alpha=0.3)
392
393 plt.tight_layout()
394 plt.savefig('curvature_drift.png', dpi=300, bbox_inches='tight')
395 plt.show()
396
397
398if __name__ == '__main__':
399 print("\n" + "="*80)
400 print("GRADIENT AND CURVATURE DRIFT SIMULATION")
401 print("="*80)
402
403 print("\nPart 1: Grad-B Drift")
404 print("-" * 80)
405 plot_grad_b_drift()
406
407 print("\nPart 2: Curvature Drift in Dipole Field")
408 print("-" * 80)
409 plot_curvature_drift()
410
411 print("\nKey observations:")
412 print(" - Grad-B drift: Opposite directions for e⁻ and p⁺")
413 print(" - E×B drift: Same direction for all species")
414 print(" - Curvature drift: Related to parallel motion along curved field lines")
415
416 print("\nDone! Generated 2 figures:")
417 print(" - grad_b_drift.png")
418 print(" - curvature_drift.png")