1#!/usr/bin/env python3
2"""
3Microwave Interferometry for Plasma Density Measurement
4
5This script demonstrates microwave interferometry diagnostics including:
6- Phase shift from line-integrated density
7- Abel inversion for radial profile reconstruction
8- Multi-chord interferometer simulation
9
10Key Physics:
11- Phase shift: Δφ = (e²/2mcε₀ω) ∫ ne dl
12- Abel inversion: reconstruct n(r) from line integrals
13
14Author: Plasma Physics Examples
15License: MIT
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20from matplotlib.gridspec import GridSpec
21from scipy.integrate import simpson
22from scipy.interpolate import interp1d
23
24# Physical constants
25QE = 1.602176634e-19 # C
26ME = 9.10938356e-31 # kg
27C = 2.99792458e8 # m/s
28EPS0 = 8.854187817e-12 # F/m
29
30def phase_shift_coefficient(freq):
31 """
32 Compute phase shift coefficient K.
33
34 Δφ = K * ∫ ne dl
35
36 where K = e²/(2·me·c·ε₀·ω) = e²·λ/(4π·me·c²·ε₀)
37
38 Parameters:
39 -----------
40 freq : float
41 Microwave frequency [Hz]
42
43 Returns:
44 --------
45 K : float
46 Phase shift coefficient [rad·m²]
47 """
48 omega = 2 * np.pi * freq
49 K = QE**2 / (2 * ME * C * EPS0 * omega)
50 return K
51
52def parabolic_density_profile(r, ne0, a):
53 """
54 Parabolic density profile (typical for tokamak).
55
56 ne(r) = ne0 * (1 - (r/a)²)^α
57
58 Parameters:
59 -----------
60 r : array
61 Radial position [m]
62 ne0 : float
63 Central density [m^-3]
64 a : float
65 Minor radius [m]
66
67 Returns:
68 --------
69 ne : array
70 Electron density [m^-3]
71 """
72 alpha = 2.0 # Profile shape parameter
73 ne = ne0 * np.maximum(0, 1 - (r / a)**2)**alpha
74 return ne
75
76def line_integrated_density(impact_param, ne0, a, num_points=500):
77 """
78 Compute line-integrated density for a given impact parameter.
79
80 ∫ ne(r) dl along chord at distance y from axis
81
82 Parameters:
83 -----------
84 impact_param : float
85 Impact parameter (distance from axis) [m]
86 ne0 : float
87 Central density [m^-3]
88 a : float
89 Minor radius [m]
90 num_points : int
91 Number of integration points
92
93 Returns:
94 --------
95 ne_line : float
96 Line-integrated density [m^-2]
97 """
98 if abs(impact_param) >= a:
99 return 0.0
100
101 # Integration along chord
102 # For impact parameter y, chord goes from -sqrt(a²-y²) to +sqrt(a²-y²)
103 x_max = np.sqrt(a**2 - impact_param**2)
104 x = np.linspace(-x_max, x_max, num_points)
105
106 # Radial distance at each point along chord
107 r = np.sqrt(x**2 + impact_param**2)
108
109 # Density at each point
110 ne = parabolic_density_profile(r, ne0, a)
111
112 # Integrate
113 ne_line = simpson(ne, x=x)
114
115 return ne_line
116
117def compute_phase_shifts(impact_params, ne0, a, freq):
118 """
119 Compute phase shifts for multiple chords.
120
121 Parameters:
122 -----------
123 impact_params : array
124 Impact parameters for each chord [m]
125 ne0 : float
126 Central density [m^-3]
127 a : float
128 Minor radius [m]
129 freq : float
130 Microwave frequency [Hz]
131
132 Returns:
133 --------
134 phase_shifts : array
135 Phase shift for each chord [rad]
136 ne_lines : array
137 Line-integrated densities [m^-2]
138 """
139 K = phase_shift_coefficient(freq)
140
141 ne_lines = np.array([line_integrated_density(y, ne0, a)
142 for y in impact_params])
143
144 phase_shifts = K * ne_lines
145
146 return phase_shifts, ne_lines
147
148def abel_inversion_matrix(impact_params, r_grid):
149 """
150 Construct Abel inversion matrix using matrix method.
151
152 For discrete measurements, we solve:
153 ne_line = A · ne_radial
154
155 where A is the Abel transform matrix.
156
157 Parameters:
158 -----------
159 impact_params : array
160 Impact parameters (sorted) [m]
161 r_grid : array
162 Radial grid points [m]
163
164 Returns:
165 --------
166 A_matrix : 2D array
167 Abel transform matrix
168 """
169 n_chords = len(impact_params)
170 n_radial = len(r_grid)
171
172 A = np.zeros((n_chords, n_radial))
173
174 # For each chord i and radial shell j
175 for i, y in enumerate(impact_params):
176 for j in range(n_radial):
177 if j == 0:
178 r_inner = 0
179 r_outer = (r_grid[0] + r_grid[1]) / 2 if n_radial > 1 else r_grid[0]
180 elif j == n_radial - 1:
181 r_inner = (r_grid[j-1] + r_grid[j]) / 2
182 r_outer = r_grid[j]
183 else:
184 r_inner = (r_grid[j-1] + r_grid[j]) / 2
185 r_outer = (r_grid[j] + r_grid[j+1]) / 2
186
187 # Chord length through shell j
188 if y < r_inner:
189 # Chord passes through entire shell
190 dl = 2 * (np.sqrt(r_outer**2 - y**2) - np.sqrt(r_inner**2 - y**2))
191 elif y < r_outer:
192 # Chord passes through outer part of shell
193 dl = 2 * np.sqrt(r_outer**2 - y**2)
194 else:
195 # Chord misses shell
196 dl = 0
197
198 A[i, j] = dl
199
200 return A
201
202def reconstruct_density_profile(phase_shifts, impact_params, a, freq, n_radial=20):
203 """
204 Reconstruct radial density profile from phase shift measurements.
205
206 Parameters:
207 -----------
208 phase_shifts : array
209 Measured phase shifts [rad]
210 impact_params : array
211 Impact parameters [m]
212 a : float
213 Minor radius [m]
214 freq : float
215 Microwave frequency [Hz]
216 n_radial : int
217 Number of radial grid points
218
219 Returns:
220 --------
221 r_grid : array
222 Radial positions [m]
223 ne_reconstructed : array
224 Reconstructed density [m^-3]
225 """
226 K = phase_shift_coefficient(freq)
227
228 # Convert phase shifts to line-integrated densities
229 ne_lines = phase_shifts / K
230
231 # Radial grid
232 r_grid = np.linspace(0, a, n_radial)
233
234 # Construct Abel matrix
235 A = abel_inversion_matrix(impact_params, r_grid)
236
237 # Solve via least squares (pseudo-inverse)
238 ne_reconstructed = np.linalg.lstsq(A, ne_lines, rcond=None)[0]
239
240 # Ensure non-negative
241 ne_reconstructed = np.maximum(ne_reconstructed, 0)
242
243 return r_grid, ne_reconstructed
244
245def plot_interferometry():
246 """
247 Create comprehensive visualization of interferometry diagnostic.
248 """
249 # Plasma parameters (tokamak-like)
250 ne0 = 5e19 # m^-3 (central density)
251 a = 0.5 # m (minor radius)
252 freq = 140e9 # Hz (140 GHz, typical for tokamak)
253
254 # Wavelength
255 wavelength = C / freq
256
257 print("=" * 70)
258 print("Microwave Interferometry Diagnostic")
259 print("=" * 70)
260 print(f"Central density: {ne0:.2e} m^-3")
261 print(f"Minor radius: {a:.2f} m")
262 print(f"Microwave frequency: {freq/1e9:.0f} GHz")
263 print(f"Wavelength: {wavelength*1e3:.2f} mm")
264 print(f"Phase shift coefficient: {phase_shift_coefficient(freq):.2e} rad·m²")
265 print("=" * 70)
266
267 # Create figure
268 fig = plt.figure(figsize=(14, 12))
269 gs = GridSpec(3, 2, figure=fig, hspace=0.4, wspace=0.35)
270
271 # True density profile
272 r_fine = np.linspace(0, a, 500)
273 ne_true = parabolic_density_profile(r_fine, ne0, a)
274
275 # Plot 1: True density profile
276 ax1 = fig.add_subplot(gs[0, 0])
277 ax1.plot(r_fine * 100, ne_true / 1e19, 'b-', linewidth=2.5)
278 ax1.set_xlabel('Radius (cm)', fontsize=11)
279 ax1.set_ylabel(r'Density ($10^{19}$ m$^{-3}$)', fontsize=11)
280 ax1.set_title('True Radial Density Profile', fontsize=12, fontweight='bold')
281 ax1.grid(True, alpha=0.3)
282 ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
283
284 # Simulate multi-chord interferometer
285 n_chords = 8
286 impact_params = np.linspace(-a * 0.9, a * 0.9, n_chords)
287
288 # Compute phase shifts
289 phase_shifts, ne_lines = compute_phase_shifts(impact_params, ne0, a, freq)
290
291 # Add some noise
292 np.random.seed(42)
293 noise_level = 0.02
294 phase_shifts_noisy = phase_shifts + np.random.normal(0, noise_level * phase_shifts.max(),
295 size=phase_shifts.shape)
296
297 # Plot 2: Chord geometry
298 ax2 = fig.add_subplot(gs[0, 1])
299
300 # Draw plasma boundary
301 theta = np.linspace(0, 2 * np.pi, 100)
302 ax2.plot(a * 100 * np.cos(theta), a * 100 * np.sin(theta),
303 'k-', linewidth=2, label='Plasma boundary')
304
305 # Draw chords
306 colors = plt.cm.viridis(np.linspace(0, 1, n_chords))
307
308 for i, (y, color) in enumerate(zip(impact_params, colors)):
309 if abs(y) < a:
310 x_max = np.sqrt(a**2 - y**2)
311 ax2.plot([-x_max * 100, x_max * 100], [y * 100, y * 100],
312 color=color, linewidth=1.5, alpha=0.7, label=f'Chord {i+1}')
313
314 ax2.set_xlabel('x (cm)', fontsize=11)
315 ax2.set_ylabel('y (cm)', fontsize=11)
316 ax2.set_title('Multi-Chord Interferometer Geometry',
317 fontsize=12, fontweight='bold')
318 ax2.set_aspect('equal')
319 ax2.grid(True, alpha=0.3)
320 ax2.legend(fontsize=7, ncol=2, loc='upper right')
321
322 # Plot 3: Phase shifts vs impact parameter
323 ax3 = fig.add_subplot(gs[1, 0])
324
325 ax3.plot(impact_params * 100, phase_shifts, 'b-', linewidth=2,
326 label='True')
327 ax3.plot(impact_params * 100, phase_shifts_noisy, 'ro', markersize=8,
328 label='Measured (with noise)')
329
330 ax3.set_xlabel('Impact Parameter (cm)', fontsize=11)
331 ax3.set_ylabel('Phase Shift (rad)', fontsize=11)
332 ax3.set_title('Measured Phase Shifts', fontsize=12, fontweight='bold')
333 ax3.grid(True, alpha=0.3)
334 ax3.legend(fontsize=10)
335
336 # Plot 4: Line-integrated density
337 ax4 = fig.add_subplot(gs[1, 1])
338
339 ax4.plot(impact_params * 100, ne_lines / 1e19, 'b-', linewidth=2,
340 label='True')
341 ax4.plot(impact_params * 100, phase_shifts_noisy / phase_shift_coefficient(freq) / 1e19,
342 'ro', markersize=8, label='From measured Δφ')
343
344 ax4.set_xlabel('Impact Parameter (cm)', fontsize=11)
345 ax4.set_ylabel(r'$\int n_e \, dl$ ($10^{19}$ m$^{-2}$)', fontsize=11)
346 ax4.set_title('Line-Integrated Density', fontsize=12, fontweight='bold')
347 ax4.grid(True, alpha=0.3)
348 ax4.legend(fontsize=10)
349
350 # Plot 5: Abel inversion reconstruction
351 ax5 = fig.add_subplot(gs[2, :])
352
353 # Reconstruct density profile
354 r_recon, ne_recon = reconstruct_density_profile(phase_shifts_noisy,
355 np.abs(impact_params),
356 a, freq, n_radial=15)
357
358 # Plot true vs reconstructed
359 ax5.plot(r_fine * 100, ne_true / 1e19, 'b-', linewidth=2.5,
360 label='True profile')
361 ax5.plot(r_recon * 100, ne_recon / 1e19, 'ro-', linewidth=2,
362 markersize=8, label=f'Reconstructed ({n_chords} chords)')
363
364 # Try with fewer chords
365 n_chords_few = 4
366 impact_params_few = np.linspace(0, a * 0.9, n_chords_few)
367 phase_shifts_few, _ = compute_phase_shifts(impact_params_few, ne0, a, freq)
368 phase_shifts_few_noisy = phase_shifts_few + np.random.normal(
369 0, noise_level * phase_shifts_few.max(), size=phase_shifts_few.shape)
370
371 r_recon_few, ne_recon_few = reconstruct_density_profile(
372 phase_shifts_few_noisy, impact_params_few, a, freq, n_radial=10)
373
374 ax5.plot(r_recon_few * 100, ne_recon_few / 1e19, 'gs-', linewidth=2,
375 markersize=8, alpha=0.7, label=f'Reconstructed ({n_chords_few} chords)')
376
377 ax5.set_xlabel('Radius (cm)', fontsize=11)
378 ax5.set_ylabel(r'Density ($10^{19}$ m$^{-3}$)', fontsize=11)
379 ax5.set_title('Abel Inversion: Reconstructed Density Profile',
380 fontsize=12, fontweight='bold')
381 ax5.grid(True, alpha=0.3)
382 ax5.legend(fontsize=10, loc='upper right')
383 ax5.set_xlim([0, a * 100])
384
385 # Add text box with reconstruction quality
386 rms_error_8 = np.sqrt(np.mean((np.interp(r_recon, r_fine, ne_true) - ne_recon)**2))
387 rms_error_4 = np.sqrt(np.mean((np.interp(r_recon_few, r_fine, ne_true) - ne_recon_few)**2))
388
389 textstr = '\n'.join([
390 'Reconstruction Quality:',
391 f'8 chords: RMS error = {rms_error_8/1e19:.2f}×10¹⁹ m⁻³',
392 f'4 chords: RMS error = {rms_error_4/1e19:.2f}×10¹⁹ m⁻³',
393 '',
394 'More chords → better reconstruction'
395 ])
396
397 ax5.text(0.98, 0.97, textstr, transform=ax5.transAxes,
398 fontsize=9, verticalalignment='top', horizontalalignment='right',
399 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
400
401 plt.suptitle('Microwave Interferometry: Density Measurement & Abel Inversion',
402 fontsize=16, fontweight='bold', y=0.995)
403
404 plt.savefig('interferometry.png', dpi=150, bbox_inches='tight')
405 print(f"\nPlot saved as 'interferometry.png'")
406
407 print(f"\nReconstruction results:")
408 print(f" 8 chords: RMS error = {rms_error_8/ne0*100:.1f}% of peak density")
409 print(f" 4 chords: RMS error = {rms_error_4/ne0*100:.1f}% of peak density")
410 print(f"\nCentral density (reconstructed, 8 chords): {ne_recon[0]:.2e} m^-3")
411 print(f"Central density (true): {ne0:.2e} m^-3")
412 print(f"Error: {abs(ne_recon[0] - ne0)/ne0*100:.1f}%")
413
414 print("=" * 70)
415
416 plt.show()
417
418if __name__ == "__main__":
419 plot_interferometry()