1#!/usr/bin/env python3
2"""
3Langmuir Probe I-V Characteristic Analysis
4
5This script generates and analyzes synthetic Langmuir probe I-V curves
6to extract plasma parameters (Te, ne, Vp).
7
8Key Physics:
9- Ion saturation: I = Isat (V << Vp)
10- Electron retardation: I ∝ exp(eV/kTe) (V < Vp)
11- Electron saturation: I → Iesat (V > Vp)
12
13Probe configurations:
14- Single probe
15- Double probe
16- Triple probe
17
18Author: Plasma Physics Examples
19License: MIT
20"""
21
22import numpy as np
23import matplotlib.pyplot as plt
24from matplotlib.gridspec import GridSpec
25from scipy.optimize import curve_fit
26
27# Physical constants
28QE = 1.602176634e-19 # C
29ME = 9.10938356e-31 # kg
30MP = 1.672621898e-27 # kg
31KB = 1.380649e-23 # J/K
32
33def ion_saturation_current(ne, Te, A_probe, mi=MP):
34 """
35 Compute ion saturation current.
36
37 Isat ≈ 0.61 * ne * e * A * sqrt(kTe/mi)
38
39 Parameters:
40 -----------
41 ne : float
42 Electron density [m^-3]
43 Te : float
44 Electron temperature [eV]
45 A_probe : float
46 Probe area [m^2]
47 mi : float
48 Ion mass [kg]
49
50 Returns:
51 --------
52 Isat : float
53 Ion saturation current [A]
54 """
55 Te_joule = Te * QE
56 cs = np.sqrt(Te_joule / mi) # Bohm velocity
57 Isat = 0.61 * ne * QE * A_probe * cs
58 return Isat
59
60def electron_saturation_current(ne, Te, A_probe):
61 """
62 Compute electron saturation current.
63
64 Iesat = ne * e * A * sqrt(kTe/(2πme))
65
66 Parameters:
67 -----------
68 ne : float
69 Electron density [m^-3]
70 Te : float
71 Electron temperature [eV]
72 A_probe : float
73 Probe area [m^2]
74
75 Returns:
76 --------
77 Iesat : float
78 Electron saturation current [A]
79 """
80 Te_joule = Te * QE
81 vth_e = np.sqrt(8 * Te_joule / (np.pi * ME))
82 Iesat = ne * QE * A_probe * vth_e / 4
83 return Iesat
84
85def langmuir_probe_iv(V, ne, Te, Vp, A_probe, mi=MP):
86 """
87 Generate idealized Langmuir probe I-V characteristic.
88
89 Parameters:
90 -----------
91 V : array
92 Probe voltage [V]
93 ne : float
94 Electron density [m^-3]
95 Te : float
96 Electron temperature [eV]
97 Vp : float
98 Plasma potential [V]
99 A_probe : float
100 Probe area [m^2]
101 mi : float
102 Ion mass [kg]
103
104 Returns:
105 --------
106 I : array
107 Probe current [A]
108 """
109 Isat = ion_saturation_current(ne, Te, A_probe, mi)
110 Iesat = electron_saturation_current(ne, Te, A_probe)
111
112 I = np.zeros_like(V)
113
114 for i, v in enumerate(V):
115 if v < Vp - 5 * Te: # Deep in ion saturation
116 I[i] = -Isat
117 elif v < Vp: # Electron retardation region
118 I[i] = -Isat + Iesat * np.exp(QE * (v - Vp) / (Te * QE))
119 else: # Electron saturation
120 I[i] = -Isat + Iesat * (1 + 0.1 * (v - Vp) / Te)
121
122 return I
123
124def add_noise(I, noise_level=0.02):
125 """Add Gaussian noise to current."""
126 noise = np.random.normal(0, noise_level * np.abs(I).max(), size=I.shape)
127 return I + noise
128
129def fit_electron_temperature(V, I, Vp_guess):
130 """
131 Fit electron temperature from electron retardation region.
132
133 ln(Ie) = ln(I0) + e·V/(kTe)
134
135 Parameters:
136 -----------
137 V : array
138 Voltage [V]
139 I : array
140 Current [A]
141 Vp_guess : float
142 Initial guess for plasma potential [V]
143
144 Returns:
145 --------
146 Te_fit : float
147 Fitted temperature [eV]
148 Vp_fit : float
149 Fitted plasma potential [V]
150 """
151 # Select retardation region (where current is positive and increasing)
152 mask = (V < Vp_guess) & (V > Vp_guess - 20) & (I > 0)
153
154 if mask.sum() < 5:
155 return None, None
156
157 V_fit = V[mask]
158 I_fit = I[mask]
159
160 # Linear fit to ln(I) vs V
161 ln_I = np.log(I_fit)
162
163 # Remove any infinities or NaNs
164 valid = np.isfinite(ln_I)
165 V_fit = V_fit[valid]
166 ln_I = ln_I[valid]
167
168 if len(V_fit) < 3:
169 return None, None
170
171 # Fit: ln(I) = a + b*V, where b = e/(kTe)
172 coeffs = np.polyfit(V_fit, ln_I, 1)
173 slope = coeffs[0]
174 intercept = coeffs[1]
175
176 # Extract Te (in eV)
177 Te_fit = 1.0 / slope # Already in eV since slope = 1/Te
178
179 # Extract Vp from where the fitted line crosses zero current
180 # Actually, find floating potential (where I=0)
181 Vf = -intercept / slope
182
183 # Plasma potential is a few Te above floating potential
184 Vp_fit = Vf + 3 * Te_fit
185
186 return Te_fit, Vp_fit
187
188def plot_langmuir_probe():
189 """
190 Create comprehensive visualization of Langmuir probe analysis.
191 """
192 # True plasma parameters
193 ne_true = 1e17 # m^-3
194 Te_true = 3.0 # eV
195 Ti_true = 0.3 # eV
196 Vp_true = 10.0 # V
197 A_probe = 1e-4 # m^2 (1 cm^2)
198
199 print("=" * 70)
200 print("Langmuir Probe I-V Characteristic Analysis")
201 print("=" * 70)
202 print("True plasma parameters:")
203 print(f" Electron density: {ne_true:.2e} m^-3")
204 print(f" Electron temperature: {Te_true:.2f} eV")
205 print(f" Plasma potential: {Vp_true:.2f} V")
206 print(f" Probe area: {A_probe*1e4:.2f} cm^2")
207 print("=" * 70)
208
209 # Generate voltage sweep
210 V = np.linspace(-30, 30, 500)
211
212 # Generate ideal I-V curve
213 I_ideal = langmuir_probe_iv(V, ne_true, Te_true, Vp_true, A_probe)
214
215 # Add noise
216 np.random.seed(42)
217 I_noisy = add_noise(I_ideal, noise_level=0.05)
218
219 # Create figure
220 fig = plt.figure(figsize=(14, 12))
221 gs = GridSpec(3, 2, figure=fig, hspace=0.4, wspace=0.35)
222
223 # Plot 1: I-V characteristic
224 ax1 = fig.add_subplot(gs[0, :])
225
226 ax1.plot(V, I_ideal * 1e3, 'b-', linewidth=2, label='Ideal')
227 ax1.plot(V, I_noisy * 1e3, 'r.', markersize=3, alpha=0.5, label='With noise')
228
229 # Mark important points
230 ax1.axvline(x=Vp_true, color='g', linestyle='--', linewidth=2,
231 label=f'Plasma potential Vp = {Vp_true:.1f} V')
232 ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
233
234 # Mark floating potential (where I=0)
235 Vf_idx = np.argmin(np.abs(I_ideal))
236 Vf = V[Vf_idx]
237 ax1.plot(Vf, 0, 'mo', markersize=10, label=f'Floating potential Vf = {Vf:.1f} V')
238
239 ax1.set_xlabel('Probe Voltage (V)', fontsize=12)
240 ax1.set_ylabel('Probe Current (mA)', fontsize=12)
241 ax1.set_title('Langmuir Probe I-V Characteristic',
242 fontsize=14, fontweight='bold')
243 ax1.grid(True, alpha=0.3)
244 ax1.legend(fontsize=10, loc='upper left')
245
246 # Annotate regions
247 ax1.text(-20, -5, 'Ion\nSaturation', fontsize=10, ha='center',
248 bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.5))
249 ax1.text(5, 10, 'Electron\nRetardation', fontsize=10, ha='center',
250 bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.5))
251 ax1.text(25, 30, 'Electron\nSaturation', fontsize=10, ha='center',
252 bbox=dict(boxstyle='round', facecolor='pink', alpha=0.5))
253
254 # Plot 2: Semi-log plot for temperature fitting
255 ax2 = fig.add_subplot(gs[1, 0])
256
257 # Plot ln(I) vs V in retardation region
258 mask_retard = (V > Vf) & (V < Vp_true) & (I_noisy > 0)
259 V_retard = V[mask_retard]
260 I_retard = I_noisy[mask_retard]
261
262 ax2.semilogy(V, np.abs(I_noisy) * 1e3, 'r.', markersize=3, alpha=0.5,
263 label='Data')
264
265 # Fit temperature
266 Te_fit, Vp_fit = fit_electron_temperature(V, I_noisy, Vp_true + 5)
267
268 if Te_fit is not None:
269 # Plot fitted line
270 V_fit_range = np.linspace(Vf, Vp_true, 100)
271 I_fit_line = np.exp((V_fit_range - (Vp_fit - 3*Te_fit)) / Te_fit) * 1e-3
272
273 ax2.semilogy(V_fit_range, I_fit_line * 1e3, 'b-', linewidth=2,
274 label=f'Fit: Te = {Te_fit:.2f} eV')
275
276 print(f"\nFitted parameters:")
277 print(f" Te (fitted): {Te_fit:.2f} eV (true: {Te_true:.2f} eV)")
278 print(f" Vp (fitted): {Vp_fit:.2f} V (true: {Vp_true:.2f} V)")
279
280 ax2.set_xlabel('Probe Voltage (V)', fontsize=12)
281 ax2.set_ylabel('|Current| (mA)', fontsize=12)
282 ax2.set_title('Semi-log Plot for Te Extraction',
283 fontsize=12, fontweight='bold')
284 ax2.grid(True, alpha=0.3, which='both')
285 ax2.legend(fontsize=10)
286 ax2.set_xlim([Vf - 5, Vp_true + 5])
287
288 # Plot 3: First derivative (dI/dV)
289 ax3 = fig.add_subplot(gs[1, 1])
290
291 dI_dV = np.gradient(I_noisy, V)
292
293 ax3.plot(V, dI_dV * 1e3, 'b-', linewidth=1.5)
294 ax3.axvline(x=Vp_true, color='g', linestyle='--', linewidth=2,
295 label='True Vp')
296
297 # Find maximum of derivative (plasma potential estimate)
298 Vp_derivative = V[np.argmax(dI_dV)]
299 ax3.axvline(x=Vp_derivative, color='r', linestyle='--', linewidth=2,
300 label=f'Vp from d²I/dV² = {Vp_derivative:.1f} V')
301
302 ax3.set_xlabel('Probe Voltage (V)', fontsize=12)
303 ax3.set_ylabel('dI/dV (mA/V)', fontsize=12)
304 ax3.set_title('First Derivative (Plasma Potential from Peak)',
305 fontsize=12, fontweight='bold')
306 ax3.grid(True, alpha=0.3)
307 ax3.legend(fontsize=10)
308
309 # Plot 4: Double probe configuration
310 ax4 = fig.add_subplot(gs[2, 0])
311
312 # Double probe: symmetric, no electron saturation
313 # Current limited by smaller of two ion saturation currents
314 I_double = np.zeros_like(V)
315 Isat = ion_saturation_current(ne_true, Te_true, A_probe)
316
317 for i, v in enumerate(V):
318 if v < -10:
319 I_double[i] = -Isat
320 elif v > 10:
321 I_double[i] = Isat
322 else:
323 # Transition region
324 I_double[i] = Isat * np.tanh(v / (2 * Te_true))
325
326 I_double_noisy = add_noise(I_double, noise_level=0.05)
327
328 ax4.plot(V, I_double_noisy * 1e3, 'b-', linewidth=1.5)
329 ax4.axhline(y=Isat * 1e3, color='r', linestyle='--', linewidth=2,
330 label=f'Ion saturation ±{Isat*1e3:.2f} mA')
331 ax4.axhline(y=-Isat * 1e3, color='r', linestyle='--', linewidth=2)
332 ax4.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
333
334 ax4.set_xlabel('Voltage Difference V1-V2 (V)', fontsize=12)
335 ax4.set_ylabel('Current (mA)', fontsize=12)
336 ax4.set_title('Double Probe Configuration',
337 fontsize=12, fontweight='bold')
338 ax4.grid(True, alpha=0.3)
339 ax4.legend(fontsize=10)
340
341 # Plot 5: Triple probe schematic and analysis
342 ax5 = fig.add_subplot(gs[2, 1])
343 ax5.axis('off')
344
345 # Triple probe explanation
346 triple_text = """
347 Triple Probe Configuration
348
349 Three electrodes:
350 • Probe 1: Floating (I₁ = 0)
351 • Probe 2: Floating (I₂ = 0)
352 • Probe 3: Draws current (I₃ = -I₁ - I₂)
353
354 Advantages:
355 ✓ No time-varying bias needed
356 ✓ Fast measurement (~μs)
357 ✓ Good for fluctuations
358
359 Formulas:
360 Te = (V₁ - V₂) / ln[(I₃+ - I₃⁺)/(I₃⁺ - I₃⁻)]
361 ne ∝ I_sat / sqrt(Te)
362
363 Limitations:
364 ✗ Assumes Maxwellian distribution
365 ✗ Requires calibration
366 ✗ More complex geometry
367 """
368
369 ax5.text(0.1, 0.95, triple_text, transform=ax5.transAxes,
370 fontsize=10, verticalalignment='top', family='monospace',
371 bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.7))
372
373 # Add simple schematic
374 from matplotlib.patches import Circle, FancyArrowPatch
375
376 # Draw three probes
377 probe_y = [0.4, 0.35, 0.3]
378 probe_x = [0.7, 0.75, 0.8]
379
380 for i, (x, y) in enumerate(zip(probe_x, probe_y)):
381 circle = Circle((x, y), 0.02, transform=ax5.transAxes,
382 facecolor='red', edgecolor='black', linewidth=2)
383 ax5.add_patch(circle)
384 ax5.text(x, y - 0.05, f'P{i+1}', transform=ax5.transAxes,
385 ha='center', fontsize=9, fontweight='bold')
386
387 plt.suptitle('Langmuir Probe Diagnostics: I-V Analysis',
388 fontsize=16, fontweight='bold', y=0.995)
389
390 plt.savefig('langmuir_probe.png', dpi=150, bbox_inches='tight')
391 print(f"\nPlot saved as 'langmuir_probe.png'")
392
393 # Error analysis
394 if Te_fit is not None:
395 Te_error = abs(Te_fit - Te_true) / Te_true * 100
396 Vp_error = abs(Vp_fit - Vp_true) / Vp_true * 100
397
398 print(f"\nFitting errors:")
399 print(f" Te error: {Te_error:.1f}%")
400 print(f" Vp error: {Vp_error:.1f}%")
401
402 print("=" * 70)
403
404 plt.show()
405
406if __name__ == "__main__":
407 plot_langmuir_probe()