1#!/usr/bin/env python3
2"""
3Diamagnetic Drift Visualization
4
5This script visualizes the diamagnetic drift, which is a fluid effect
6(not a real particle drift) arising from pressure gradients.
7
8Key Physics:
9- Diamagnetic current: J* = B × ∇p / B²
10- Ion and electron drifts in opposite directions
11- NOT a real particle drift - fluid velocity only
12- Compare with E×B drift (same for both species)
13
14Author: Plasma Physics Examples
15License: MIT
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20from matplotlib.gridspec import GridSpec
21import matplotlib.patches as mpatches
22
23# Physical constants
24QE = 1.602176634e-19 # C
25ME = 9.10938356e-31 # kg
26MP = 1.672621898e-27 # kg
27
28def density_profile(x, n0, Ln):
29 """
30 Exponential density profile.
31
32 n(x) = n0 * exp(-x/Ln)
33
34 Parameters:
35 -----------
36 x : array
37 Position [m]
38 n0 : float
39 Central density [m^-3]
40 Ln : float
41 Density scale length [m]
42
43 Returns:
44 --------
45 n : array
46 Density [m^-3]
47 """
48 return n0 * np.exp(-x / Ln)
49
50def pressure_gradient(x, n0, T, Ln):
51 """
52 Pressure gradient for exponential profile.
53
54 p = n·kT
55 ∇p = -kT·n0/Ln·exp(-x/Ln)
56
57 Parameters:
58 -----------
59 x : array
60 Position [m]
61 n0 : float
62 Central density [m^-3]
63 T : float
64 Temperature [eV]
65 Ln : float
66 Scale length [m]
67
68 Returns:
69 --------
70 grad_p : array
71 Pressure gradient [Pa/m]
72 """
73 T_joule = T * QE
74 n = density_profile(x, n0, Ln)
75 grad_p = -T_joule * n / Ln
76 return grad_p
77
78def diamagnetic_drift_velocity(grad_p, B0, charge_sign=1):
79 """
80 Compute diamagnetic drift velocity.
81
82 v* = (B × ∇p) / (qnB²)
83
84 For electrons: charge_sign = -1
85 For ions: charge_sign = +1
86
87 Parameters:
88 -----------
89 grad_p : array
90 Pressure gradient [Pa/m]
91 B0 : float
92 Magnetic field [T]
93 charge_sign : int
94 +1 for ions, -1 for electrons
95
96 Returns:
97 --------
98 v_star : array
99 Diamagnetic drift velocity [m/s]
100 """
101 # For gradient in x direction, B in z direction
102 # v* is in y direction
103 v_star = -charge_sign * grad_p / (B0**2)
104
105 return v_star
106
107def visualize_gyro_orbits(x_center, n, B0, T, m, charge_sign, num_particles=5):
108 """
109 Visualize gyro orbits at different x positions to show mechanism.
110
111 Parameters:
112 -----------
113 x_center : array
114 Center positions
115 n : array
116 Density at each position
117 B0 : float
118 Magnetic field
119 T : float
120 Temperature [eV]
121 m : float
122 Particle mass
123 charge_sign : int
124 ±1
125 num_particles : int
126 Number of particles to show
127
128 Returns:
129 --------
130 x_orbits, y_orbits : lists of orbit coordinates
131 """
132 T_joule = T * QE
133 vth = np.sqrt(2 * T_joule / m)
134 omega_c = abs(charge_sign * QE * B0 / m)
135 rho = vth / omega_c
136
137 x_orbits = []
138 y_orbits = []
139
140 # Sample particles at different x positions
141 x_samples = np.linspace(x_center.min(), x_center.max(), num_particles)
142
143 for x_c in x_samples:
144 # Gyro orbit
145 theta = np.linspace(0, 2 * np.pi, 100)
146
147 if charge_sign > 0:
148 # Ion gyrates clockwise (looking along B)
149 x_orbit = x_c + rho * np.cos(theta)
150 y_orbit = rho * np.sin(theta)
151 else:
152 # Electron gyrates counterclockwise
153 x_orbit = x_c - rho * np.cos(theta)
154 y_orbit = rho * np.sin(theta)
155
156 x_orbits.append(x_orbit)
157 y_orbits.append(y_orbit)
158
159 return x_orbits, y_orbits, rho
160
161def plot_diamagnetic_drift():
162 """
163 Create comprehensive visualization of diamagnetic drift.
164 """
165 # Plasma parameters
166 n0 = 1e19 # m^-3
167 Te = 10.0 # eV
168 Ti = 10.0 # eV
169 B0 = 1.0 # T
170 Ln = 0.05 # 5 cm scale length
171 mi = MP
172
173 print("=" * 70)
174 print("Diamagnetic Drift Parameters")
175 print("=" * 70)
176 print(f"Central density: {n0:.2e} m^-3")
177 print(f"Electron temperature: {Te:.1f} eV")
178 print(f"Ion temperature: {Ti:.1f} eV")
179 print(f"Magnetic field: {B0:.2f} T")
180 print(f"Density scale length: {Ln*100:.1f} cm")
181 print("=" * 70)
182
183 # Position array
184 x = np.linspace(0, 0.2, 1000) # 0 to 20 cm
185
186 # Density and pressure profiles
187 ne = density_profile(x, n0, Ln)
188 ni = ne # Quasineutrality
189
190 grad_pe = pressure_gradient(x, n0, Te, Ln)
191 grad_pi = pressure_gradient(x, n0, Ti, Ln)
192
193 # Diamagnetic velocities (in y-direction)
194 # For gradient in -x direction: ∇p points in -x
195 # B in +z direction
196 # B × ∇p points in ±y direction
197
198 v_star_e = diamagnetic_drift_velocity(grad_pe, B0, charge_sign=-1)
199 v_star_i = diamagnetic_drift_velocity(grad_pi, B0, charge_sign=+1)
200
201 # Diamagnetic current
202 J_star = QE * (ni * v_star_i - ne * v_star_e)
203
204 # For comparison: E×B drift (same for both species)
205 # Assume some electric field
206 E_field = 1000 # V/m (arbitrary)
207 v_ExB = E_field / B0
208
209 # Create figure
210 fig = plt.figure(figsize=(14, 12))
211 gs = GridSpec(4, 2, figure=fig, hspace=0.4, wspace=0.3)
212
213 # Plot 1: Density profile
214 ax1 = fig.add_subplot(gs[0, 0])
215 ax1.plot(x * 100, ne / 1e19, 'b-', linewidth=2)
216 ax1.set_xlabel('Position x (cm)', fontsize=11)
217 ax1.set_ylabel(r'Density ($10^{19}$ m$^{-3}$)', fontsize=11)
218 ax1.set_title('Density Profile (∇n points in -x)', fontsize=12, fontweight='bold')
219 ax1.grid(True, alpha=0.3)
220 ax1.axvline(x=Ln * 100, color='r', linestyle='--', linewidth=1,
221 label=f'Scale length = {Ln*100:.1f} cm')
222 ax1.legend(fontsize=9)
223
224 # Plot 2: Pressure gradient
225 ax2 = fig.add_subplot(gs[0, 1])
226 ax2.plot(x * 100, grad_pe / 1e6, 'b-', linewidth=2, label='Electron')
227 ax2.plot(x * 100, grad_pi / 1e6, 'r-', linewidth=2, label='Ion')
228 ax2.set_xlabel('Position x (cm)', fontsize=11)
229 ax2.set_ylabel(r'∇p (MPa/m)', fontsize=11)
230 ax2.set_title('Pressure Gradient', fontsize=12, fontweight='bold')
231 ax2.grid(True, alpha=0.3)
232 ax2.legend(fontsize=9)
233
234 # Plot 3: Diamagnetic drift velocities
235 ax3 = fig.add_subplot(gs[1, :])
236 ax3.plot(x * 100, v_star_e / 1e3, 'b-', linewidth=2.5, label='Electron v*e (fluid)')
237 ax3.plot(x * 100, v_star_i / 1e3, 'r-', linewidth=2.5, label='Ion v*i (fluid)')
238 ax3.axhline(y=v_ExB / 1e3, color='g', linestyle='--', linewidth=2,
239 label=f'E×B drift = {v_ExB/1e3:.1f} km/s (both species)')
240 ax3.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
241
242 ax3.set_xlabel('Position x (cm)', fontsize=11)
243 ax3.set_ylabel('Drift Velocity (km/s)', fontsize=11)
244 ax3.set_title('Diamagnetic Drift vs E×B Drift', fontsize=12, fontweight='bold')
245 ax3.grid(True, alpha=0.3)
246 ax3.legend(fontsize=10, loc='upper right')
247
248 # Add text annotation
249 ax3.text(0.02, 0.95, 'Diamagnetic drifts are OPPOSITE\nE×B drifts are SAME',
250 transform=ax3.transAxes, fontsize=10, verticalalignment='top',
251 bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))
252
253 # Plot 4: Diamagnetic current
254 ax4 = fig.add_subplot(gs[2, 0])
255 ax4.plot(x * 100, J_star / 1e6, 'purple', linewidth=2.5)
256 ax4.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
257 ax4.set_xlabel('Position x (cm)', fontsize=11)
258 ax4.set_ylabel(r'Current Density (MA/m$^2$)', fontsize=11)
259 ax4.set_title('Diamagnetic Current J* = (B × ∇p)/B²',
260 fontsize=12, fontweight='bold')
261 ax4.grid(True, alpha=0.3)
262
263 # Plot 5: Mechanism illustration (gyro orbits)
264 ax5 = fig.add_subplot(gs[2, 1])
265
266 # Show more particles on one side than the other
267 x_left = 0.05 # High density
268 x_right = 0.15 # Low density
269
270 n_left = density_profile(x_left, n0, Ln)
271 n_right = density_profile(x_right, n0, Ln)
272
273 num_left = int(20 * n_left / n0)
274 num_right = int(20 * n_right / n0)
275
276 # Electron Larmor radius
277 Te_joule = Te * QE
278 vth_e = np.sqrt(2 * Te_joule / ME)
279 omega_ce = QE * B0 / ME
280 rho_e = vth_e / omega_ce
281
282 # Plot electron gyro orbits
283 y_positions_left = np.linspace(-0.01, 0.01, num_left)
284 y_positions_right = np.linspace(-0.01, 0.01, num_right)
285
286 for y_pos in y_positions_left:
287 theta = np.linspace(0, 2 * np.pi, 50)
288 x_orbit = x_left - rho_e * np.cos(theta)
289 y_orbit = y_pos + rho_e * np.sin(theta)
290 ax5.plot(x_orbit * 100, y_orbit * 100, 'b-', alpha=0.3, linewidth=0.5)
291
292 for y_pos in y_positions_right:
293 theta = np.linspace(0, 2 * np.pi, 50)
294 x_orbit = x_right - rho_e * np.cos(theta)
295 y_orbit = y_pos + rho_e * np.sin(theta)
296 ax5.plot(x_orbit * 100, y_orbit * 100, 'b-', alpha=0.3, linewidth=0.5)
297
298 # Mark regions
299 ax5.axvline(x=x_left * 100, color='red', linestyle='--', linewidth=2,
300 label='High density')
301 ax5.axvline(x=x_right * 100, color='blue', linestyle='--', linewidth=2,
302 label='Low density')
303
304 # Arrow showing net drift
305 ax5.annotate('Net drift →', xy=(0.5, 0.95), xytext=(0.2, 0.95),
306 xycoords='axes fraction',
307 arrowprops=dict(arrowstyle='->', lw=2, color='green'),
308 fontsize=11, fontweight='bold')
309
310 ax5.set_xlabel('x (cm)', fontsize=11)
311 ax5.set_ylabel('y (cm)', fontsize=11)
312 ax5.set_title('Mechanism: More Orbits on Left → Net Drift',
313 fontsize=12, fontweight='bold')
314 ax5.legend(fontsize=9)
315 ax5.set_aspect('equal')
316
317 # Plot 6: Summary comparison table
318 ax6 = fig.add_subplot(gs[3, :])
319 ax6.axis('off')
320
321 # Create comparison table
322 comparison_data = [
323 ['Property', 'Diamagnetic Drift', 'E×B Drift'],
324 ['Direction', 'Opposite for e⁻ and ions', 'Same for both'],
325 ['Real drift?', 'NO (fluid only)', 'YES (particle drift)'],
326 ['Formula', 'v* = (B × ∇p)/(qnB²)', 'v = E×B/B²'],
327 ['Current', 'J* = (B × ∇p)/B²', 'J = 0 (both drift same)'],
328 ['Origin', 'Pressure gradient', 'Electric field'],
329 ]
330
331 # Draw table
332 table = ax6.table(cellText=comparison_data, cellLoc='left',
333 bbox=[0, 0, 1, 1])
334 table.auto_set_font_size(False)
335 table.set_fontsize(10)
336 table.scale(1, 2)
337
338 # Style header row
339 for i in range(3):
340 table[(0, i)].set_facecolor('#4CAF50')
341 table[(0, i)].set_text_props(weight='bold', color='white')
342
343 # Alternate row colors
344 for i in range(1, len(comparison_data)):
345 for j in range(3):
346 if i % 2 == 0:
347 table[(i, j)].set_facecolor('#f0f0f0')
348
349 ax6.set_title('Comparison: Diamagnetic vs E×B Drift',
350 fontsize=12, fontweight='bold', pad=20)
351
352 plt.suptitle('Diamagnetic Drift: A Fluid Effect',
353 fontsize=16, fontweight='bold', y=0.998)
354
355 plt.savefig('diamagnetic_drift.png', dpi=150, bbox_inches='tight')
356 print("\nPlot saved as 'diamagnetic_drift.png'")
357 print("\nKey insight:")
358 print(" Diamagnetic drift is NOT a real particle drift!")
359 print(" It's a fluid velocity arising from density gradients.")
360 print(" Individual particles don't drift diamagnetically,")
361 print(" but the fluid (average) velocity does.")
362
363 plt.show()
364
365if __name__ == "__main__":
366 plot_diamagnetic_drift()