1#!/usr/bin/env python3
2"""
3Generalized Ohm's Law: Term-by-Term Comparison
4
5This script compares the relative magnitude of different terms in
6generalized Ohm's law for various plasma regimes.
7
8Generalized Ohm's law:
9E + v×B = ηJ + (J×B)/(ne) - ∇pe/(ne) + (me/ne²)(dJ/dt)
10
11Key Physics:
12- Ideal MHD: E + v×B = 0 (only this term)
13- Resistive MHD: + ηJ
14- Hall MHD: + J×B/(ne)
15- Electron pressure: - ∇pe/(ne)
16- Electron inertia: + (me/ne²)dJ/dt
17
18Author: Plasma Physics Examples
19License: MIT
20"""
21
22import numpy as np
23import matplotlib.pyplot as plt
24from matplotlib.gridspec import GridSpec
25
26# Physical constants
27QE = 1.602176634e-19 # C
28ME = 9.10938356e-31 # kg
29MP = 1.672621898e-27 # kg
30EPS0 = 8.854187817e-12 # F/m
31MU0 = 4 * np.pi * 1e-7 # H/m
32
33def compute_plasma_parameters(ne, B0, Te, L, v):
34 """
35 Compute characteristic parameters for a plasma.
36
37 Parameters:
38 -----------
39 ne : float
40 Electron density [m^-3]
41 B0 : float
42 Magnetic field [T]
43 Te : float
44 Electron temperature [eV]
45 L : float
46 Characteristic length scale [m]
47 v : float
48 Characteristic velocity [m/s]
49
50 Returns:
51 --------
52 dict : Plasma parameters
53 """
54 mi = MP
55
56 # Characteristic current density (from force balance)
57 J = B0 / (MU0 * L)
58
59 # Ion skin depth
60 di = np.sqrt(ME * MU0 / (QE * ne)) / np.sqrt(MP / ME) # Corrected
61 di = np.sqrt(MP / (MU0 * ne * QE**2))
62
63 # Resistivity (Spitzer)
64 Te_joule = Te * QE
65 ln_Lambda = 15 # Coulomb logarithm
66 eta = 5.2e-5 * ln_Lambda / (Te**1.5) # Ω·m
67
68 # Electron thermal velocity
69 vth_e = np.sqrt(2 * Te_joule / ME)
70
71 # Electron plasma frequency
72 omega_pe = np.sqrt(ne * QE**2 / (ME * EPS0))
73
74 # Characteristic time
75 tau = L / v
76
77 return {
78 'ne': ne,
79 'B0': B0,
80 'Te': Te,
81 'L': L,
82 'v': v,
83 'J': J,
84 'di': di,
85 'eta': eta,
86 'vth_e': vth_e,
87 'omega_pe': omega_pe,
88 'tau': tau
89 }
90
91def ideal_mhd_term(params):
92 """E + v×B term magnitude."""
93 return params['v'] * params['B0']
94
95def resistive_term(params):
96 """η·J term magnitude."""
97 return params['eta'] * params['J']
98
99def hall_term(params):
100 """(J×B)/(ne) term magnitude."""
101 return params['J'] * params['B0'] / (params['ne'] * QE)
102
103def pressure_term(params):
104 """∇pe/(ne) term magnitude."""
105 Te_joule = params['Te'] * QE
106 # ∇p ~ p/L = nkT/L
107 return Te_joule / params['L']
108
109def inertia_term(params):
110 """(me/ne²)(dJ/dt) term magnitude."""
111 # dJ/dt ~ J/τ
112 dJ_dt = params['J'] / params['tau']
113 return (ME / (params['ne'] * QE)) * dJ_dt
114
115def plot_ohms_law_comparison():
116 """
117 Compare Ohm's law terms for different plasma regimes.
118 """
119 print("=" * 70)
120 print("Generalized Ohm's Law: Term Comparison")
121 print("=" * 70)
122
123 # Create figure
124 fig = plt.figure(figsize=(14, 12))
125 gs = GridSpec(3, 2, figure=fig, hspace=0.4, wspace=0.35)
126
127 # Define plasma regimes
128 regimes = {
129 'Tokamak Core': {
130 'ne': 1e20, # m^-3
131 'B0': 5.0, # T
132 'Te': 10000, # eV (10 keV)
133 'L': 1.0, # m
134 'v': 1e5, # m/s
135 'color': 'red'
136 },
137 'Solar Wind': {
138 'ne': 1e7, # m^-3
139 'B0': 5e-9, # T (5 nT)
140 'Te': 10, # eV
141 'L': 1e6, # m (1000 km)
142 'v': 4e5, # m/s (400 km/s)
143 'color': 'orange'
144 },
145 'Magnetopause': {
146 'ne': 1e7, # m^-3
147 'B0': 50e-9, # T (50 nT)
148 'Te': 100, # eV
149 'L': 1e5, # m (100 km)
150 'v': 1e5, # m/s
151 'color': 'blue'
152 },
153 'Ionosphere': {
154 'ne': 1e12, # m^-3
155 'B0': 50e-6, # T (50 μT)
156 'Te': 0.1, # eV
157 'L': 1e4, # m (10 km)
158 'v': 1e3, # m/s
159 'color': 'green'
160 },
161 'Lab Plasma': {
162 'ne': 1e18, # m^-3
163 'B0': 0.1, # T
164 'Te': 10, # eV
165 'L': 0.1, # m
166 'v': 1e4, # m/s
167 'color': 'purple'
168 }
169 }
170
171 # Plot 1: Term comparison for each regime (bar chart)
172 ax1 = fig.add_subplot(gs[0, :])
173
174 regime_names = list(regimes.keys())
175 x_pos = np.arange(len(regime_names))
176 width = 0.15
177
178 terms_data = {name: [] for name in regime_names}
179 term_names = ['Ideal', 'Resistive', 'Hall', 'Pressure', 'Inertia']
180
181 for regime_name, regime_params in regimes.items():
182 params = compute_plasma_parameters(**regime_params)
183
184 # Compute terms
185 E_ideal = ideal_mhd_term(params)
186 E_resist = resistive_term(params)
187 E_hall = hall_term(params)
188 E_press = pressure_term(params)
189 E_inert = inertia_term(params)
190
191 # Normalize by ideal term
192 terms_data[regime_name] = [
193 1.0, # Ideal
194 E_resist / E_ideal,
195 E_hall / E_ideal,
196 E_press / E_ideal,
197 E_inert / E_ideal
198 ]
199
200 # Plot bars
201 for i, term_name in enumerate(term_names):
202 values = [terms_data[name][i] for name in regime_names]
203 ax1.bar(x_pos + i * width, values, width, label=term_name)
204
205 ax1.set_ylabel('Magnitude (normalized to Ideal MHD)', fontsize=11)
206 ax1.set_title('Relative Magnitude of Ohm\'s Law Terms by Regime',
207 fontsize=13, fontweight='bold')
208 ax1.set_xticks(x_pos + width * 2)
209 ax1.set_xticklabels(regime_names, rotation=15, ha='right')
210 ax1.set_yscale('log')
211 ax1.legend(fontsize=9, ncol=5, loc='upper right')
212 ax1.grid(True, alpha=0.3, axis='y')
213 ax1.axhline(y=1, color='k', linestyle='--', linewidth=1)
214
215 # Plot 2: Terms vs scale length (tokamak)
216 ax2 = fig.add_subplot(gs[1, 0])
217
218 L_array = np.logspace(-2, 2, 100) # 1 cm to 100 m
219 tokamak_base = regimes['Tokamak Core'].copy()
220
221 terms_vs_L = {name: [] for name in term_names}
222
223 for L in L_array:
224 tokamak_base['L'] = L
225 params = compute_plasma_parameters(**tokamak_base)
226
227 E_ideal = ideal_mhd_term(params)
228 terms_vs_L['Ideal'].append(E_ideal)
229 terms_vs_L['Resistive'].append(resistive_term(params))
230 terms_vs_L['Hall'].append(hall_term(params))
231 terms_vs_L['Pressure'].append(pressure_term(params))
232 terms_vs_L['Inertia'].append(inertia_term(params))
233
234 # Normalize by ideal
235 for term_name in term_names[1:]: # Skip ideal
236 normalized = np.array(terms_vs_L[term_name]) / np.array(terms_vs_L['Ideal'])
237 ax2.loglog(L_array / params['di'], normalized, linewidth=2, label=term_name)
238
239 ax2.axhline(y=1, color='k', linestyle='--', linewidth=1,
240 label='Equal to Ideal')
241 ax2.set_xlabel(r'$L / d_i$ (scale length / ion skin depth)', fontsize=11)
242 ax2.set_ylabel('Term / Ideal MHD Term', fontsize=11)
243 ax2.set_title('Tokamak: Terms vs Scale Length', fontsize=12, fontweight='bold')
244 ax2.legend(fontsize=9)
245 ax2.grid(True, alpha=0.3, which='both')
246
247 # Plot 3: Terms vs scale length (magnetopause)
248 ax3 = fig.add_subplot(gs[1, 1])
249
250 L_array_mp = np.logspace(3, 7, 100) # 1 km to 10,000 km
251 mp_base = regimes['Magnetopause'].copy()
252
253 terms_vs_L_mp = {name: [] for name in term_names}
254
255 for L in L_array_mp:
256 mp_base['L'] = L
257 params = compute_plasma_parameters(**mp_base)
258
259 E_ideal = ideal_mhd_term(params)
260 terms_vs_L_mp['Ideal'].append(E_ideal)
261 terms_vs_L_mp['Resistive'].append(resistive_term(params))
262 terms_vs_L_mp['Hall'].append(hall_term(params))
263 terms_vs_L_mp['Pressure'].append(pressure_term(params))
264 terms_vs_L_mp['Inertia'].append(inertia_term(params))
265
266 # Normalize by ideal
267 for term_name in term_names[1:]:
268 normalized = np.array(terms_vs_L_mp[term_name]) / np.array(terms_vs_L_mp['Ideal'])
269 ax3.loglog(L_array_mp / params['di'], normalized, linewidth=2, label=term_name)
270
271 ax3.axhline(y=1, color='k', linestyle='--', linewidth=1,
272 label='Equal to Ideal')
273 ax3.set_xlabel(r'$L / d_i$ (scale length / ion skin depth)', fontsize=11)
274 ax3.set_ylabel('Term / Ideal MHD Term', fontsize=11)
275 ax3.set_title('Magnetopause: Terms vs Scale Length',
276 fontsize=12, fontweight='bold')
277 ax3.legend(fontsize=9)
278 ax3.grid(True, alpha=0.3, which='both')
279
280 # Plot 4: Summary table
281 ax4 = fig.add_subplot(gs[2, :])
282 ax4.axis('off')
283
284 # Create table data
285 table_data = [
286 ['Term', 'Formula', 'When Important', 'Example'],
287 ['Ideal MHD', 'E + v×B = 0', 'L >> di, collisional', 'MHD turbulence'],
288 ['Resistive', '+ η·J', 'Collisional plasma', 'Resistive instabilities'],
289 ['Hall', '+ J×B/(ne)', 'L ~ di', 'Magnetic reconnection'],
290 ['Pressure', '- ∇pe/(ne)', 'Strong gradients', 'Current sheets'],
291 ['Inertia', '+ (me/ne²)·dJ/dt', 'Fast dynamics (ω ~ ωpe)', 'Beam instabilities'],
292 ]
293
294 table = ax4.table(cellText=table_data, cellLoc='left',
295 bbox=[0, 0, 1, 1])
296 table.auto_set_font_size(False)
297 table.set_fontsize(9)
298 table.scale(1, 2.5)
299
300 # Style header
301 for i in range(4):
302 table[(0, i)].set_facecolor('#2196F3')
303 table[(0, i)].set_text_props(weight='bold', color='white')
304
305 # Color rows
306 for i in range(1, len(table_data)):
307 for j in range(4):
308 if i % 2 == 0:
309 table[(i, j)].set_facecolor('#f5f5f5')
310
311 ax4.set_title('Generalized Ohm\'s Law: When Each Term Matters',
312 fontsize=12, fontweight='bold', pad=10)
313
314 plt.suptitle('Generalized Ohm\'s Law Term-by-Term Analysis',
315 fontsize=16, fontweight='bold', y=0.995)
316
317 plt.savefig('ohms_law_comparison.png', dpi=150, bbox_inches='tight')
318 print("\nPlot saved as 'ohms_law_comparison.png'\n")
319
320 # Print detailed results
321 print("\nDetailed Results by Regime:")
322 print("-" * 70)
323
324 for regime_name, regime_params in regimes.items():
325 params = compute_plasma_parameters(**regime_params)
326
327 print(f"\n{regime_name}:")
328 print(f" ne = {params['ne']:.2e} m^-3")
329 print(f" B0 = {params['B0']:.2e} T")
330 print(f" Te = {params['Te']:.2e} eV")
331 print(f" L = {params['L']:.2e} m")
332 print(f" L/di = {params['L']/params['di']:.2f}")
333 print(f" η = {params['eta']:.2e} Ω·m")
334
335 E_ideal = ideal_mhd_term(params)
336 E_resist = resistive_term(params)
337 E_hall = hall_term(params)
338 E_press = pressure_term(params)
339 E_inert = inertia_term(params)
340
341 print(f"\n Term magnitudes (V/m):")
342 print(f" Ideal MHD: {E_ideal:.2e}")
343 print(f" Resistive: {E_resist:.2e} ({E_resist/E_ideal:.2e} × Ideal)")
344 print(f" Hall: {E_hall:.2e} ({E_hall/E_ideal:.2e} × Ideal)")
345 print(f" Pressure: {E_press:.2e} ({E_press/E_ideal:.2e} × Ideal)")
346 print(f" Inertia: {E_inert:.2e} ({E_inert/E_ideal:.2e} × Ideal)")
347
348 # Determine dominant terms
349 terms = {
350 'Resistive': E_resist / E_ideal,
351 'Hall': E_hall / E_ideal,
352 'Pressure': E_press / E_ideal,
353 'Inertia': E_inert / E_ideal
354 }
355 important_terms = [name for name, ratio in terms.items() if ratio > 0.1]
356
357 if important_terms:
358 print(f"\n Important corrections: {', '.join(important_terms)}")
359 else:
360 print(f"\n Regime: Pure Ideal MHD")
361
362 print("\n" + "=" * 70)
363
364 plt.show()
365
366if __name__ == "__main__":
367 plot_ohms_law_comparison()