1#!/usr/bin/env python3
2"""
3Spectral Line Analysis for Plasma Temperature and Density
4
5This script demonstrates spectral line diagnostics including:
6- Doppler broadening (temperature measurement)
7- Stark broadening (density measurement)
8- Voigt profile fitting
9
10Key Physics:
11- Doppler width: ΔλD/λ = sqrt(2kTi/mic²)
12- Stark width: proportional to electron density
13- Line shape: convolution of Gaussian (Doppler) + Lorentzian (Stark)
14
15Author: Plasma Physics Examples
16License: MIT
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from matplotlib.gridspec import GridSpec
22from scipy.special import wofz
23from scipy.optimize import curve_fit
24
25# Physical constants
26QE = 1.602176634e-19 # C
27ME = 9.10938356e-31 # kg
28MP = 1.672621898e-27 # kg
29C = 2.99792458e8 # m/s
30KB = 1.380649e-23 # J/K
31
32def gaussian_profile(wavelength, lambda0, Ti, mi):
33 """
34 Gaussian (Doppler) line profile.
35
36 I(λ) ∝ exp(-(λ-λ0)²/(2σD²))
37
38 where σD = (λ0/c) * sqrt(2kTi/mi)
39
40 Parameters:
41 -----------
42 wavelength : array
43 Wavelength [m]
44 lambda0 : float
45 Central wavelength [m]
46 Ti : float
47 Ion temperature [eV]
48 mi : float
49 Ion mass [kg]
50
51 Returns:
52 --------
53 profile : array
54 Normalized intensity
55 """
56 Ti_joule = Ti * QE
57
58 # Doppler width
59 sigma_D = (lambda0 / C) * np.sqrt(2 * Ti_joule / mi)
60
61 # Gaussian profile
62 profile = np.exp(-(wavelength - lambda0)**2 / (2 * sigma_D**2))
63 profile /= profile.max() # Normalize
64
65 return profile
66
67def lorentzian_profile(wavelength, lambda0, ne):
68 """
69 Lorentzian (Stark) line profile.
70
71 I(λ) ∝ γ² / ((λ-λ0)² + γ²)
72
73 where γ is the Stark width, proportional to ne.
74
75 Parameters:
76 -----------
77 wavelength : array
78 Wavelength [m]
79 lambda0 : float
80 Central wavelength [m]
81 ne : float
82 Electron density [m^-3]
83
84 Returns:
85 --------
86 profile : array
87 Normalized intensity
88 """
89 # Stark width (empirical for hydrogen Balmer-α at 656.3 nm)
90 # γ ≈ 1.0e-26 * ne^(2/3) [m]
91 gamma = 1.0e-26 * ne**(2/3)
92
93 # Lorentzian profile
94 profile = gamma**2 / ((wavelength - lambda0)**2 + gamma**2)
95 profile /= profile.max() # Normalize
96
97 return profile
98
99def voigt_profile(wavelength, lambda0, Ti, mi, ne):
100 """
101 Voigt profile: convolution of Gaussian and Lorentzian.
102
103 Uses the Faddeeva function for efficient computation.
104
105 Parameters:
106 -----------
107 wavelength : array
108 Wavelength [m]
109 lambda0 : float
110 Central wavelength [m]
111 Ti : float
112 Ion temperature [eV]
113 mi : float
114 Ion mass [kg]
115 ne : float
116 Electron density [m^-3]
117
118 Returns:
119 --------
120 profile : array
121 Normalized intensity
122 """
123 Ti_joule = Ti * QE
124
125 # Doppler width
126 sigma_D = (lambda0 / C) * np.sqrt(2 * Ti_joule / mi)
127
128 # Stark width
129 gamma = 1.0e-26 * ne**(2/3)
130
131 # Voigt profile using Faddeeva function
132 z = ((wavelength - lambda0) + 1j * gamma) / (sigma_D * np.sqrt(2))
133 profile = np.real(wofz(z)) / (sigma_D * np.sqrt(2 * np.pi))
134
135 profile /= profile.max() # Normalize
136
137 return profile
138
139def instrument_broadening(wavelength, profile, sigma_inst):
140 """
141 Apply instrumental broadening (Gaussian convolution).
142
143 Parameters:
144 -----------
145 wavelength : array
146 Wavelength [m]
147 profile : array
148 Original profile
149 sigma_inst : float
150 Instrumental width [m]
151
152 Returns:
153 --------
154 broadened : array
155 Broadened profile
156 """
157 # Gaussian kernel
158 kernel = np.exp(-wavelength**2 / (2 * sigma_inst**2))
159 kernel /= kernel.sum()
160
161 # Convolve
162 broadened = np.convolve(profile, kernel, mode='same')
163 broadened /= broadened.max()
164
165 return broadened
166
167def fit_voigt_profile(wavelength, intensity, lambda0_guess, Ti_guess, mi, ne_guess):
168 """
169 Fit Voigt profile to extract Ti and ne.
170
171 Parameters:
172 -----------
173 wavelength : array
174 Wavelength [m]
175 intensity : array
176 Measured intensity
177 lambda0_guess : float
178 Initial guess for central wavelength [m]
179 Ti_guess : float
180 Initial guess for temperature [eV]
181 mi : float
182 Ion mass [kg]
183 ne_guess : float
184 Initial guess for density [m^-3]
185
186 Returns:
187 --------
188 Ti_fit, ne_fit, lambda0_fit : fitted parameters
189 """
190 # Define fitting function
191 def voigt_fit_func(lam, lambda0, Ti, ne):
192 return voigt_profile(lam, lambda0, Ti, mi, ne)
193
194 # Initial guess
195 p0 = [lambda0_guess, Ti_guess, ne_guess]
196
197 # Bounds
198 bounds = ([lambda0_guess - 1e-10, 0.1, 1e16],
199 [lambda0_guess + 1e-10, 1000, 1e21])
200
201 # Fit
202 try:
203 popt, _ = curve_fit(voigt_fit_func, wavelength, intensity,
204 p0=p0, bounds=bounds, maxfev=10000)
205 lambda0_fit, Ti_fit, ne_fit = popt
206 return Ti_fit, ne_fit, lambda0_fit
207 except:
208 return None, None, None
209
210def plot_spectral_line_analysis():
211 """
212 Create comprehensive visualization of spectral line analysis.
213 """
214 # Plasma parameters (true values)
215 Ti_true = 10.0 # eV
216 ne_true = 5e17 # m^-3
217 mi = MP # Hydrogen
218 lambda0 = 656.3e-9 # m (H-alpha line)
219
220 # Instrumental width
221 sigma_inst = 0.05e-9 # m (50 pm resolution)
222
223 print("=" * 70)
224 print("Spectral Line Analysis: Doppler and Stark Broadening")
225 print("=" * 70)
226 print(f"True parameters:")
227 print(f" Ion temperature: {Ti_true:.1f} eV")
228 print(f" Electron density: {ne_true:.2e} m^-3")
229 print(f" Line: H-alpha at {lambda0*1e9:.1f} nm")
230 print(f" Instrumental resolution: {sigma_inst*1e12:.0f} pm")
231 print("=" * 70)
232
233 # Create wavelength array
234 dlambda = 2e-9 # ±2 nm around line center
235 wavelength = np.linspace(lambda0 - dlambda, lambda0 + dlambda, 2000)
236
237 # Create figure
238 fig = plt.figure(figsize=(14, 12))
239 gs = GridSpec(3, 2, figure=fig, hspace=0.4, wspace=0.35)
240
241 # Plot 1: Individual broadening mechanisms
242 ax1 = fig.add_subplot(gs[0, :])
243
244 # Gaussian (Doppler)
245 profile_gaussian = gaussian_profile(wavelength, lambda0, Ti_true, mi)
246
247 # Lorentzian (Stark)
248 profile_lorentzian = lorentzian_profile(wavelength, lambda0, ne_true)
249
250 # Voigt (combined)
251 profile_voigt = voigt_profile(wavelength, lambda0, Ti_true, mi, ne_true)
252
253 ax1.plot((wavelength - lambda0) * 1e12, profile_gaussian, 'b-', linewidth=2,
254 label='Doppler (Gaussian)')
255 ax1.plot((wavelength - lambda0) * 1e12, profile_lorentzian, 'r-', linewidth=2,
256 label='Stark (Lorentzian)')
257 ax1.plot((wavelength - lambda0) * 1e12, profile_voigt, 'g-', linewidth=2.5,
258 label='Voigt (convolution)')
259
260 ax1.set_xlabel('Wavelength - λ₀ (pm)', fontsize=12)
261 ax1.set_ylabel('Normalized Intensity', fontsize=12)
262 ax1.set_title('Line Broadening Mechanisms', fontsize=13, fontweight='bold')
263 ax1.grid(True, alpha=0.3)
264 ax1.legend(fontsize=11)
265 ax1.set_xlim([-2000, 2000])
266
267 # Plot 2: Effect of temperature on Doppler width
268 ax2 = fig.add_subplot(gs[1, 0])
269
270 Ti_values = [1, 5, 10, 50] # eV
271 colors_T = plt.cm.Reds(np.linspace(0.3, 1, len(Ti_values)))
272
273 for Ti, color in zip(Ti_values, colors_T):
274 profile = gaussian_profile(wavelength, lambda0, Ti, mi)
275 ax2.plot((wavelength - lambda0) * 1e12, profile, color=color,
276 linewidth=2, label=f'Ti = {Ti} eV')
277
278 ax2.set_xlabel('Wavelength - λ₀ (pm)', fontsize=11)
279 ax2.set_ylabel('Normalized Intensity', fontsize=11)
280 ax2.set_title('Doppler Broadening vs Temperature',
281 fontsize=12, fontweight='bold')
282 ax2.grid(True, alpha=0.3)
283 ax2.legend(fontsize=10)
284 ax2.set_xlim([-1000, 1000])
285
286 # Plot 3: Effect of density on Stark width
287 ax3 = fig.add_subplot(gs[1, 1])
288
289 ne_values = [1e17, 5e17, 1e18, 5e18] # m^-3
290 colors_n = plt.cm.Blues(np.linspace(0.3, 1, len(ne_values)))
291
292 for ne, color in zip(ne_values, colors_n):
293 profile = lorentzian_profile(wavelength, lambda0, ne)
294 ax3.plot((wavelength - lambda0) * 1e12, profile, color=color,
295 linewidth=2, label=f'ne = {ne:.1e} m⁻³')
296
297 ax3.set_xlabel('Wavelength - λ₀ (pm)', fontsize=11)
298 ax3.set_ylabel('Normalized Intensity', fontsize=11)
299 ax3.set_title('Stark Broadening vs Density', fontsize=12, fontweight='bold')
300 ax3.grid(True, alpha=0.3)
301 ax3.legend(fontsize=10)
302 ax3.set_xlim([-500, 500])
303
304 # Plot 4: Measured spectrum with noise and fitting
305 ax4 = fig.add_subplot(gs[2, 0])
306
307 # Generate synthetic measured spectrum
308 profile_measured = voigt_profile(wavelength, lambda0, Ti_true, mi, ne_true)
309
310 # Apply instrumental broadening
311 profile_measured = instrument_broadening(wavelength, profile_measured, sigma_inst)
312
313 # Add noise
314 np.random.seed(42)
315 noise = np.random.normal(0, 0.02, size=profile_measured.shape)
316 profile_measured_noisy = profile_measured + noise
317
318 # Plot measured data
319 ax4.plot((wavelength - lambda0) * 1e12, profile_measured_noisy, 'k.',
320 markersize=2, alpha=0.5, label='Measured (with noise)')
321
322 # Fit Voigt profile
323 Ti_fit, ne_fit, lambda0_fit = fit_voigt_profile(
324 wavelength, profile_measured_noisy, lambda0, Ti_true, mi, ne_true)
325
326 if Ti_fit is not None:
327 profile_fit = voigt_profile(wavelength, lambda0_fit, Ti_fit, mi, ne_fit)
328
329 ax4.plot((wavelength - lambda0) * 1e12, profile_fit, 'r-', linewidth=2,
330 label=f'Fit: Ti={Ti_fit:.1f} eV, ne={ne_fit:.2e}')
331
332 print(f"\nFitting results:")
333 print(f" Ti (fitted): {Ti_fit:.2f} eV (true: {Ti_true:.1f} eV)")
334 print(f" ne (fitted): {ne_fit:.2e} m^-3 (true: {ne_true:.2e} m^-3)")
335 print(f" λ0 (fitted): {lambda0_fit*1e9:.4f} nm (true: {lambda0*1e9:.4f} nm)")
336
337 Ti_error = abs(Ti_fit - Ti_true) / Ti_true * 100
338 ne_error = abs(ne_fit - ne_true) / ne_true * 100
339
340 print(f"\nFitting errors:")
341 print(f" Ti error: {Ti_error:.1f}%")
342 print(f" ne error: {ne_error:.1f}%")
343
344 ax4.set_xlabel('Wavelength - λ₀ (pm)', fontsize=11)
345 ax4.set_ylabel('Normalized Intensity', fontsize=11)
346 ax4.set_title('Voigt Profile Fitting', fontsize=12, fontweight='bold')
347 ax4.grid(True, alpha=0.3)
348 ax4.legend(fontsize=10)
349 ax4.set_xlim([-1000, 1000])
350
351 # Plot 5: Width extraction method
352 ax5 = fig.add_subplot(gs[2, 1])
353
354 # Plot FWHM measurement
355 profile_clean = voigt_profile(wavelength, lambda0, Ti_true, mi, ne_true)
356
357 # Find FWHM
358 half_max = profile_clean.max() / 2
359 idx_half = np.where(profile_clean >= half_max)[0]
360 lambda_left = wavelength[idx_half[0]]
361 lambda_right = wavelength[idx_half[-1]]
362 fwhm = lambda_right - lambda_left
363
364 ax5.plot((wavelength - lambda0) * 1e12, profile_clean, 'b-', linewidth=2)
365 ax5.axhline(y=half_max, color='r', linestyle='--', linewidth=1,
366 label=f'FWHM = {fwhm*1e12:.1f} pm')
367 ax5.axvline(x=(lambda_left - lambda0) * 1e12, color='r', linestyle=':', linewidth=1)
368 ax5.axvline(x=(lambda_right - lambda0) * 1e12, color='r', linestyle=':', linewidth=1)
369
370 # Mark FWHM points
371 ax5.plot([(lambda_left - lambda0) * 1e12, (lambda_right - lambda0) * 1e12],
372 [half_max, half_max], 'ro', markersize=8)
373
374 # Theoretical Doppler FWHM
375 Ti_joule = Ti_true * QE
376 sigma_D = (lambda0 / C) * np.sqrt(2 * Ti_joule / mi)
377 fwhm_doppler = 2 * np.sqrt(2 * np.log(2)) * sigma_D
378
379 ax5.text(0.05, 0.95, f'FWHM (measured): {fwhm*1e12:.1f} pm\n'
380 f'FWHM (Doppler only): {fwhm_doppler*1e12:.1f} pm\n'
381 f'Broadening ratio: {fwhm/fwhm_doppler:.2f}',
382 transform=ax5.transAxes, fontsize=9, verticalalignment='top',
383 bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.6))
384
385 ax5.set_xlabel('Wavelength - λ₀ (pm)', fontsize=11)
386 ax5.set_ylabel('Normalized Intensity', fontsize=11)
387 ax5.set_title('FWHM Measurement', fontsize=12, fontweight='bold')
388 ax5.grid(True, alpha=0.3)
389 ax5.legend(fontsize=10)
390 ax5.set_xlim([-800, 800])
391
392 plt.suptitle('Spectral Line Analysis: Temperature and Density Diagnostics',
393 fontsize=16, fontweight='bold', y=0.995)
394
395 plt.savefig('spectral_line_analysis.png', dpi=150, bbox_inches='tight')
396 print(f"\nPlot saved as 'spectral_line_analysis.png'")
397 print("=" * 70)
398
399 plt.show()
400
401if __name__ == "__main__":
402 plot_spectral_line_analysis()