1#!/usr/bin/env python3
2"""
3Tokamak Profile Analysis
4
5Generates and analyzes standard tokamak radial profiles including:
6- Pressure p(ψ)
7- Safety factor q(ψ)
8- Current density J(ψ)
9- Density n(ψ)
10- Temperature T(ψ)
11
12Compares L-mode and H-mode (with pedestal) profiles.
13Computes key dimensionless parameters: βN, li, βp.
14
15Author: Claude
16Date: 2026-02-14
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from scipy.integrate import simpson
22from scipy.constants import mu_0, e, m_p
23
24# Physical constants
25MU_0 = mu_0
26E_CHARGE = e
27M_PROTON = m_p
28
29
30class TokamakProfiles:
31 """
32 Class to generate and analyze tokamak equilibrium profiles.
33
34 Attributes
35 ----------
36 psi_norm : ndarray
37 Normalized flux coordinate (0 at axis, 1 at edge)
38 """
39
40 def __init__(self, n_points=200):
41 """
42 Initialize profile generator.
43
44 Parameters
45 ----------
46 n_points : int
47 Number of radial points
48 """
49 self.psi_norm = np.linspace(0, 1, n_points)
50 self.n_points = n_points
51
52 def temperature_profile(self, T0, T_edge, mode='L', alpha_T=2.0):
53 """
54 Generate temperature profile.
55
56 Parameters
57 ----------
58 T0 : float
59 Central temperature (keV)
60 T_edge : float
61 Edge temperature (keV)
62 mode : str
63 'L' for L-mode, 'H' for H-mode with pedestal
64 alpha_T : float
65 Profile shape parameter
66
67 Returns
68 -------
69 T : ndarray
70 Temperature profile (keV)
71 """
72 psi = self.psi_norm
73
74 if mode == 'L':
75 # L-mode: smooth parabolic profile
76 T = T_edge + (T0 - T_edge) * (1 - psi**alpha_T)
77
78 elif mode == 'H':
79 # H-mode: core + pedestal
80 T_ped = 0.2 * T0 # Pedestal temperature
81 ped_width = 0.05
82 ped_position = 0.9
83
84 # Core profile
85 T_core = T_ped + (T0 - T_ped) * (1 - (psi/ped_position)**alpha_T)
86 T_core = np.where(psi < ped_position, T_core, T_ped)
87
88 # Pedestal drop
89 pedestal_mask = (psi >= ped_position) & (psi <= ped_position + ped_width)
90 T_drop = T_ped + (T_edge - T_ped) * ((psi - ped_position) / ped_width)
91 T = np.where(pedestal_mask, T_drop, T_core)
92 T = np.where(psi > ped_position + ped_width, T_edge, T)
93
94 else:
95 raise ValueError(f"Unknown mode: {mode}")
96
97 return T
98
99 def density_profile(self, n0, n_edge, mode='L', alpha_n=0.5):
100 """
101 Generate density profile.
102
103 Parameters
104 ----------
105 n0 : float
106 Central density (10^19 m^-3)
107 n_edge : float
108 Edge density (10^19 m^-3)
109 mode : str
110 'L' or 'H' mode
111 alpha_n : float
112 Profile shape parameter
113
114 Returns
115 -------
116 n : ndarray
117 Density profile (10^19 m^-3)
118 """
119 psi = self.psi_norm
120
121 if mode == 'L':
122 n = n_edge + (n0 - n_edge) * (1 - psi**alpha_n)
123
124 elif mode == 'H':
125 n_ped = 0.5 * n0
126 ped_width = 0.05
127 ped_position = 0.9
128
129 n_core = n_ped + (n0 - n_ped) * (1 - (psi/ped_position)**alpha_n)
130 n_core = np.where(psi < ped_position, n_core, n_ped)
131
132 pedestal_mask = (psi >= ped_position) & (psi <= ped_position + ped_width)
133 n_drop = n_ped + (n_edge - n_ped) * ((psi - ped_position) / ped_width)
134 n = np.where(pedestal_mask, n_drop, n_core)
135 n = np.where(psi > ped_position + ped_width, n_edge, n)
136
137 return n
138
139 def pressure_profile(self, T, n):
140 """
141 Compute pressure from temperature and density.
142
143 p = 2 * n * k_B * T (factor 2 for electrons + ions)
144
145 Parameters
146 ----------
147 T : ndarray
148 Temperature (keV)
149 n : ndarray
150 Density (10^19 m^-3)
151
152 Returns
153 -------
154 p : ndarray
155 Pressure (kPa)
156 """
157 # Convert to SI units
158 T_J = T * 1e3 * E_CHARGE # keV -> J
159 n_SI = n * 1e19 # 10^19 m^-3 -> m^-3
160
161 p = 2 * n_SI * T_J # Pa
162 return p / 1e3 # Convert to kPa
163
164 def safety_factor_profile(self, q0, q_edge, mode='monotonic'):
165 """
166 Generate safety factor profile.
167
168 Parameters
169 ----------
170 q0 : float
171 Central safety factor
172 q_edge : float
173 Edge safety factor
174 mode : str
175 'monotonic' or 'reversed' (for advanced scenarios)
176
177 Returns
178 -------
179 q : ndarray
180 Safety factor profile
181 """
182 psi = self.psi_norm
183
184 if mode == 'monotonic':
185 # Standard monotonic q profile
186 q = q0 + (q_edge - q0) * psi**2
187
188 elif mode == 'reversed':
189 # Reversed shear (for ITB scenarios)
190 q_min = 1.5
191 psi_min = 0.4
192 q = np.where(psi < psi_min,
193 q0 + (q_min - q0) * (psi / psi_min)**2,
194 q_min + (q_edge - q_min) * ((psi - psi_min) / (1 - psi_min))**2)
195
196 return q
197
198 def current_density_profile(self, q, B0, R0, a):
199 """
200 Compute current density from safety factor.
201
202 J ~ dq/dψ (simplified relation)
203
204 Parameters
205 ----------
206 q : ndarray
207 Safety factor
208 B0 : float
209 Toroidal field (T)
210 R0 : float
211 Major radius (m)
212 a : float
213 Minor radius (m)
214
215 Returns
216 -------
217 J : ndarray
218 Current density (MA/m^2)
219 """
220 psi = self.psi_norm
221 r = psi * a # Approximate minor radius
222
223 # J_phi ~ (r*B0)/(R0*q) for circular geometry
224 J = np.zeros_like(psi)
225 J[1:] = (r[1:] * B0) / (R0 * q[1:] + 1e-10)
226 J = J / (a**2) # Normalize
227
228 # Convert to MA/m^2
229 J = J * 10 # Scaling factor
230
231 return J
232
233
234def compute_beta_parameters(p, B0, R0, a, I_p):
235 """
236 Compute dimensionless beta parameters.
237
238 Parameters
239 ----------
240 p : ndarray
241 Pressure profile (Pa)
242 B0 : float
243 Toroidal field (T)
244 R0 : float
245 Major radius (m)
246 a : float
247 Minor radius (m)
248 I_p : float
249 Plasma current (MA)
250
251 Returns
252 -------
253 dict
254 Dictionary with beta_N, beta_p, beta_toroidal, li
255 """
256 # Volume-averaged pressure
257 p_SI = p * 1e3 # kPa -> Pa
258 p_avg = np.mean(p_SI)
259
260 # Toroidal beta
261 beta_toroidal = 2 * MU_0 * p_avg / B0**2
262
263 # Poloidal beta (simplified)
264 B_p = MU_0 * I_p * 1e6 / (2 * np.pi * a) # Poloidal field estimate
265 beta_p = 2 * MU_0 * p_avg / B_p**2
266
267 # Normalized beta
268 beta_N = beta_toroidal / (I_p / (a * B0))
269
270 # Internal inductance (simplified, assuming circular cross-section)
271 li = 1.0 # Typical value
272
273 return {
274 'beta_toroidal': beta_toroidal * 100, # Convert to %
275 'beta_p': beta_p,
276 'beta_N': beta_N,
277 'li': li
278 }
279
280
281def plot_profiles(profiles_L, profiles_H, psi_norm):
282 """
283 Plot comparison of L-mode and H-mode profiles.
284
285 Parameters
286 ----------
287 profiles_L : dict
288 L-mode profiles
289 profiles_H : dict
290 H-mode profiles
291 psi_norm : ndarray
292 Normalized flux coordinate
293 """
294 fig, axes = plt.subplots(2, 3, figsize=(15, 10))
295
296 # Temperature
297 axes[0, 0].plot(psi_norm, profiles_L['T'], 'b-', linewidth=2, label='L-mode')
298 axes[0, 0].plot(psi_norm, profiles_H['T'], 'r-', linewidth=2, label='H-mode')
299 axes[0, 0].set_xlabel('Normalized flux ψ_N', fontsize=11)
300 axes[0, 0].set_ylabel('Temperature (keV)', fontsize=11)
301 axes[0, 0].set_title('Temperature Profile', fontsize=12, fontweight='bold')
302 axes[0, 0].legend(fontsize=10)
303 axes[0, 0].grid(True, alpha=0.3)
304
305 # Density
306 axes[0, 1].plot(psi_norm, profiles_L['n'], 'b-', linewidth=2, label='L-mode')
307 axes[0, 1].plot(psi_norm, profiles_H['n'], 'r-', linewidth=2, label='H-mode')
308 axes[0, 1].set_xlabel('Normalized flux ψ_N', fontsize=11)
309 axes[0, 1].set_ylabel('Density (10¹⁹ m⁻³)', fontsize=11)
310 axes[0, 1].set_title('Density Profile', fontsize=12, fontweight='bold')
311 axes[0, 1].legend(fontsize=10)
312 axes[0, 1].grid(True, alpha=0.3)
313
314 # Pressure
315 axes[0, 2].plot(psi_norm, profiles_L['p'], 'b-', linewidth=2, label='L-mode')
316 axes[0, 2].plot(psi_norm, profiles_H['p'], 'r-', linewidth=2, label='H-mode')
317 axes[0, 2].set_xlabel('Normalized flux ψ_N', fontsize=11)
318 axes[0, 2].set_ylabel('Pressure (kPa)', fontsize=11)
319 axes[0, 2].set_title('Pressure Profile', fontsize=12, fontweight='bold')
320 axes[0, 2].legend(fontsize=10)
321 axes[0, 2].grid(True, alpha=0.3)
322
323 # Safety factor
324 axes[1, 0].plot(psi_norm, profiles_L['q'], 'b-', linewidth=2, label='L-mode')
325 axes[1, 0].plot(psi_norm, profiles_H['q'], 'r-', linewidth=2, label='H-mode')
326 axes[1, 0].axhline(y=1, color='k', linestyle='--', alpha=0.3)
327 axes[1, 0].axhline(y=2, color='k', linestyle='--', alpha=0.3)
328 axes[1, 0].set_xlabel('Normalized flux ψ_N', fontsize=11)
329 axes[1, 0].set_ylabel('Safety factor q', fontsize=11)
330 axes[1, 0].set_title('Safety Factor', fontsize=12, fontweight='bold')
331 axes[1, 0].legend(fontsize=10)
332 axes[1, 0].grid(True, alpha=0.3)
333
334 # Current density
335 axes[1, 1].plot(psi_norm, profiles_L['J'], 'b-', linewidth=2, label='L-mode')
336 axes[1, 1].plot(psi_norm, profiles_H['J'], 'r-', linewidth=2, label='H-mode')
337 axes[1, 1].set_xlabel('Normalized flux ψ_N', fontsize=11)
338 axes[1, 1].set_ylabel('J (MA/m²)', fontsize=11)
339 axes[1, 1].set_title('Current Density', fontsize=12, fontweight='bold')
340 axes[1, 1].legend(fontsize=10)
341 axes[1, 1].grid(True, alpha=0.3)
342
343 # Pressure gradient (related to stability)
344 dp_L = np.gradient(profiles_L['p'], psi_norm)
345 dp_H = np.gradient(profiles_H['p'], psi_norm)
346 axes[1, 2].plot(psi_norm, np.abs(dp_L), 'b-', linewidth=2, label='L-mode')
347 axes[1, 2].plot(psi_norm, np.abs(dp_H), 'r-', linewidth=2, label='H-mode')
348 axes[1, 2].set_xlabel('Normalized flux ψ_N', fontsize=11)
349 axes[1, 2].set_ylabel('|dp/dψ| (kPa)', fontsize=11)
350 axes[1, 2].set_title('Pressure Gradient', fontsize=12, fontweight='bold')
351 axes[1, 2].legend(fontsize=10)
352 axes[1, 2].grid(True, alpha=0.3)
353
354 plt.suptitle('Tokamak Radial Profiles: L-mode vs H-mode',
355 fontsize=14, fontweight='bold')
356 plt.tight_layout()
357 plt.savefig('tokamak_profiles.png', dpi=150, bbox_inches='tight')
358 plt.show()
359
360
361def main():
362 """Main execution function."""
363 # Initialize profile generator
364 profiles_gen = TokamakProfiles(n_points=200)
365
366 # Tokamak parameters
367 R0 = 1.8 # Major radius (m)
368 a = 0.5 # Minor radius (m)
369 B0 = 3.0 # Toroidal field (T)
370 I_p = 1.0 # Plasma current (MA)
371
372 # Profile parameters
373 T0 = 10.0 # Central temperature (keV)
374 T_edge = 0.1 # Edge temperature (keV)
375 n0 = 8.0 # Central density (10^19 m^-3)
376 n_edge = 1.0 # Edge density (10^19 m^-3)
377 q0 = 1.0 # Central q
378 q_edge = 4.0 # Edge q
379
380 print("Tokamak Profile Analysis")
381 print("=" * 60)
382 print(f"Device parameters:")
383 print(f" Major radius R0: {R0:.2f} m")
384 print(f" Minor radius a: {a:.2f} m")
385 print(f" Toroidal field B0: {B0:.1f} T")
386 print(f" Plasma current I_p: {I_p:.1f} MA")
387 print()
388
389 # Generate L-mode profiles
390 T_L = profiles_gen.temperature_profile(T0, T_edge, mode='L', alpha_T=2.0)
391 n_L = profiles_gen.density_profile(n0, n_edge, mode='L', alpha_n=0.5)
392 p_L = profiles_gen.pressure_profile(T_L, n_L)
393 q_L = profiles_gen.safety_factor_profile(q0, q_edge, mode='monotonic')
394 J_L = profiles_gen.current_density_profile(q_L, B0, R0, a)
395
396 # Generate H-mode profiles
397 T_H = profiles_gen.temperature_profile(T0, T_edge, mode='H', alpha_T=2.0)
398 n_H = profiles_gen.density_profile(n0, n_edge, mode='H', alpha_n=0.5)
399 p_H = profiles_gen.pressure_profile(T_H, n_H)
400 q_H = profiles_gen.safety_factor_profile(q0, q_edge, mode='monotonic')
401 J_H = profiles_gen.current_density_profile(q_H, B0, R0, a)
402
403 profiles_L = {'T': T_L, 'n': n_L, 'p': p_L, 'q': q_L, 'J': J_L}
404 profiles_H = {'T': T_H, 'n': n_H, 'p': p_H, 'q': q_H, 'J': J_H}
405
406 # Compute beta parameters
407 beta_L = compute_beta_parameters(p_L, B0, R0, a, I_p)
408 beta_H = compute_beta_parameters(p_H, B0, R0, a, I_p)
409
410 print("L-mode parameters:")
411 print(f" β_toroidal: {beta_L['beta_toroidal']:.3f}%")
412 print(f" β_poloidal: {beta_L['beta_p']:.3f}")
413 print(f" β_N: {beta_L['beta_N']:.3f}")
414 print(f" li: {beta_L['li']:.2f}")
415 print()
416
417 print("H-mode parameters:")
418 print(f" β_toroidal: {beta_H['beta_toroidal']:.3f}%")
419 print(f" β_poloidal: {beta_H['beta_p']:.3f}")
420 print(f" β_N: {beta_H['beta_N']:.3f}")
421 print(f" li: {beta_H['li']:.2f}")
422 print()
423
424 # Plot profiles
425 plot_profiles(profiles_L, profiles_H, profiles_gen.psi_norm)
426 print("Profile comparison plot saved as 'tokamak_profiles.png'")
427
428
429if __name__ == '__main__':
430 main()