1#!/usr/bin/env python3
2"""
3Kruskal-Shafranov Stability Criterion
4
5Analyzes the Kruskal-Shafranov criterion for kink instability:
6 q(a) > m/n
7
8where q(a) is the safety factor at the plasma edge, and (m,n) are
9the poloidal and toroidal mode numbers.
10
11Computes stability boundaries and maximum stable current vs aspect ratio.
12
13Author: Claude
14Date: 2026-02-14
15"""
16
17import numpy as np
18import matplotlib.pyplot as plt
19from scipy.integrate import simpson
20
21
22def current_profile(r, a, I_total, profile='parabolic', alpha=2.0):
23 """
24 Define current density profile.
25
26 Parameters
27 ----------
28 r : array_like
29 Radial positions
30 a : float
31 Plasma radius
32 I_total : float
33 Total plasma current
34 profile : str
35 Profile type: 'uniform', 'parabolic', 'peaked'
36 alpha : float
37 Profile peaking factor
38
39 Returns
40 -------
41 J : ndarray
42 Current density
43 """
44 if profile == 'uniform':
45 J_0 = I_total / (np.pi * a**2)
46 J = np.where(r <= a, J_0, 0.0)
47
48 elif profile == 'parabolic':
49 J_0 = (alpha + 1) * I_total / (np.pi * a**2)
50 J = np.where(r <= a, J_0 * (1 - (r/a)**2)**alpha, 0.0)
51
52 elif profile == 'peaked':
53 # Highly peaked on axis
54 J_0 = 3 * I_total / (np.pi * a**2)
55 J = np.where(r <= a, J_0 * np.exp(-3*(r/a)**2), 0.0)
56
57 else:
58 raise ValueError(f"Unknown profile: {profile}")
59
60 return J
61
62
63def safety_factor_from_current(r, J, R0, B0, a):
64 """
65 Compute safety factor q(r) from current density.
66
67 For circular cross-section:
68 q(r) = (r * B_φ) / (R0 * B_θ)
69
70 where B_θ(r) = μ₀ I(r) / (2π r)
71
72 Parameters
73 ----------
74 r : ndarray
75 Radial positions
76 J : ndarray
77 Current density
78 R0 : float
79 Major radius
80 B0 : float
81 Toroidal field on axis
82 a : float
83 Minor radius
84
85 Returns
86 -------
87 q : ndarray
88 Safety factor
89 """
90 mu_0 = 4 * np.pi * 1e-7
91
92 q = np.zeros_like(r)
93
94 for i, ri in enumerate(r):
95 if ri < 1e-10:
96 q[i] = 0.0
97 else:
98 # Enclosed current
99 r_int = r[r <= ri]
100 J_int = J[r <= ri]
101 if len(r_int) > 1:
102 I_enclosed = 2 * np.pi * simpson(r_int * J_int, x=r_int)
103 else:
104 I_enclosed = 0.0
105
106 # Poloidal field
107 B_theta = mu_0 * I_enclosed / (2 * np.pi * ri)
108
109 # Toroidal field (assume ~ 1/R)
110 B_phi = B0 * R0 / (R0 + ri * np.cos(0)) # At outboard midplane
111
112 # Safety factor
113 q[i] = ri * B_phi / (R0 * B_theta + 1e-10)
114
115 return q
116
117
118def kruskal_shafranov_criterion(q_edge, m, n):
119 """
120 Check Kruskal-Shafranov criterion.
121
122 Criterion: q(a) > m/n for stability
123
124 Parameters
125 ----------
126 q_edge : float
127 Safety factor at edge
128 m : int
129 Poloidal mode number
130 n : int
131 Toroidal mode number
132
133 Returns
134 -------
135 bool
136 True if stable
137 """
138 return q_edge > m / n
139
140
141def maximum_stable_current(R0, a, B0, m, n, profile='parabolic'):
142 """
143 Compute maximum stable current for given mode (m,n).
144
145 Parameters
146 ----------
147 R0 : float
148 Major radius
149 a : float
150 Minor radius
151 B0 : float
152 Toroidal field
153 m, n : int
154 Mode numbers
155 profile : str
156 Current profile type
157
158 Returns
159 -------
160 I_max : float
161 Maximum stable current
162 """
163 # Binary search for maximum current
164 I_min, I_max = 1e3, 10e6 # Search range: 1 kA to 10 MA
165 tolerance = 1e3 # 1 kA
166
167 r = np.linspace(0, a, 200)
168
169 while I_max - I_min > tolerance:
170 I_test = (I_min + I_max) / 2
171
172 J = current_profile(r, a, I_test, profile=profile)
173 q = safety_factor_from_current(r, J, R0, B0, a)
174 q_edge = q[-1]
175
176 if kruskal_shafranov_criterion(q_edge, m, n):
177 # Stable, can increase current
178 I_min = I_test
179 else:
180 # Unstable, decrease current
181 I_max = I_test
182
183 return I_min
184
185
186def plot_q_profile(r, q, a, m_values=[1, 2, 3], n=1):
187 """
188 Plot q profile with stability boundaries.
189
190 Parameters
191 ----------
192 r : ndarray
193 Radial grid
194 q : ndarray
195 Safety factor
196 a : float
197 Minor radius
198 m_values : list
199 Mode numbers to show
200 n : int
201 Toroidal mode number
202 """
203 fig, ax = plt.subplots(figsize=(10, 7))
204
205 # Normalized radius
206 r_norm = r / a
207
208 # Plot q profile
209 ax.plot(r_norm, q, 'b-', linewidth=3, label='q(r)')
210
211 # Plot stability boundaries
212 colors = ['red', 'orange', 'green', 'purple', 'brown']
213 for i, m in enumerate(m_values):
214 q_crit = m / n
215 ax.axhline(y=q_crit, color=colors[i % len(colors)],
216 linestyle='--', linewidth=2,
217 label=f'q = {m}/{n} (m={m}, n={n})')
218
219 # Mark q at edge
220 q_edge = q[-1]
221 ax.plot(1.0, q_edge, 'ro', markersize=12, label=f'q(a) = {q_edge:.2f}')
222
223 ax.set_xlabel('Normalized radius r/a', fontsize=13)
224 ax.set_ylabel('Safety factor q', fontsize=13)
225 ax.set_title('Safety Factor Profile with Stability Boundaries',
226 fontsize=14, fontweight='bold')
227 ax.legend(fontsize=11, loc='best')
228 ax.grid(True, alpha=0.3)
229 ax.set_xlim([0, 1])
230 ax.set_ylim([0, np.max(q) * 1.2])
231
232 plt.tight_layout()
233 plt.savefig('q_profile.png', dpi=150, bbox_inches='tight')
234 plt.show()
235
236
237def plot_stability_diagram(R0, B0, profile='parabolic'):
238 """
239 Plot stability diagram: I_max vs aspect ratio.
240
241 Parameters
242 ----------
243 R0 : float
244 Major radius
245 B0 : float
246 Toroidal field
247 profile : str
248 Current profile type
249 """
250 # Scan over aspect ratio
251 epsilon_values = np.linspace(0.1, 0.5, 20) # a/R0
252 a_values = epsilon_values * R0
253
254 # Compute max current for different modes
255 modes = [(1, 1), (2, 1), (3, 1), (3, 2)]
256 I_max_data = {mode: [] for mode in modes}
257
258 print("Computing stability boundaries...")
259 for a in a_values:
260 for mode in modes:
261 m, n = mode
262 I_max = maximum_stable_current(R0, a, B0, m, n, profile=profile)
263 I_max_data[mode].append(I_max / 1e6) # Convert to MA
264
265 # Plot
266 fig, ax = plt.subplots(figsize=(11, 7))
267
268 colors = ['red', 'blue', 'green', 'purple']
269 markers = ['o', 's', '^', 'D']
270
271 for i, mode in enumerate(modes):
272 m, n = mode
273 ax.plot(epsilon_values, I_max_data[mode],
274 color=colors[i], marker=markers[i], markersize=7,
275 linewidth=2, label=f'm/n = {m}/{n}')
276
277 ax.set_xlabel('Inverse aspect ratio ε = a/R₀', fontsize=13)
278 ax.set_ylabel('Maximum stable current I_p (MA)', fontsize=13)
279 ax.set_title('Kruskal-Shafranov Stability Limit',
280 fontsize=14, fontweight='bold')
281 ax.legend(fontsize=11, loc='best')
282 ax.grid(True, alpha=0.3)
283
284 # Add text box with parameters
285 textstr = f'R₀ = {R0:.1f} m\nB₀ = {B0:.1f} T\nProfile: {profile}'
286 props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
287 ax.text(0.05, 0.95, textstr, transform=ax.transAxes,
288 fontsize=11, verticalalignment='top', bbox=props)
289
290 plt.tight_layout()
291 plt.savefig('stability_boundary.png', dpi=150, bbox_inches='tight')
292 plt.show()
293
294
295def parameter_scan(R0, a, B0, profile='parabolic'):
296 """
297 Scan current and plot marginal stability curves.
298
299 Parameters
300 ----------
301 R0 : float
302 Major radius
303 a : float
304 Minor radius
305 B0 : float
306 Toroidal field
307 profile : str
308 Current profile type
309 """
310 # Current scan
311 I_values = np.linspace(0.1e6, 5e6, 30) # 0.1 to 5 MA
312 r = np.linspace(0, a, 200)
313
314 q_edge_values = []
315
316 for I in I_values:
317 J = current_profile(r, a, I, profile=profile)
318 q = safety_factor_from_current(r, J, R0, B0, a)
319 q_edge_values.append(q[-1])
320
321 # Plot
322 fig, ax = plt.subplots(figsize=(10, 7))
323
324 ax.plot(I_values / 1e6, q_edge_values, 'b-', linewidth=3,
325 label='q(a) vs I_p')
326
327 # Mark stability boundaries
328 modes = [(1, 1), (2, 1), (3, 1), (3, 2)]
329 colors = ['red', 'orange', 'green', 'purple']
330
331 for i, (m, n) in enumerate(modes):
332 q_crit = m / n
333 ax.axhline(y=q_crit, color=colors[i], linestyle='--',
334 linewidth=2, label=f'm/n = {m}/{n}')
335
336 ax.set_xlabel('Plasma current I_p (MA)', fontsize=13)
337 ax.set_ylabel('Edge safety factor q(a)', fontsize=13)
338 ax.set_title('Marginal Stability: q(a) vs Plasma Current',
339 fontsize=14, fontweight='bold')
340 ax.legend(fontsize=11, loc='best')
341 ax.grid(True, alpha=0.3)
342
343 # Shade unstable regions
344 for i, (m, n) in enumerate(modes):
345 q_crit = m / n
346 if i == 0: # Only shade below q=1 for clarity
347 ax.fill_between(I_values / 1e6, 0, q_crit,
348 where=(np.array(q_edge_values) < q_crit),
349 alpha=0.2, color='red',
350 label='Unstable region')
351
352 plt.tight_layout()
353 plt.savefig('marginal_stability.png', dpi=150, bbox_inches='tight')
354 plt.show()
355
356
357def main():
358 """Main execution function."""
359 # Tokamak parameters
360 R0 = 1.5 # Major radius (m)
361 a = 0.5 # Minor radius (m)
362 B0 = 2.5 # Toroidal field (T)
363 I_p = 1.0e6 # Plasma current (A) = 1 MA
364 profile = 'parabolic'
365
366 print("Kruskal-Shafranov Stability Analysis")
367 print("=" * 60)
368 print(f"Tokamak parameters:")
369 print(f" Major radius R₀: {R0:.2f} m")
370 print(f" Minor radius a: {a:.2f} m")
371 print(f" Aspect ratio A = R₀/a: {R0/a:.2f}")
372 print(f" Toroidal field B₀: {B0:.2f} T")
373 print(f" Plasma current I_p: {I_p/1e6:.2f} MA")
374 print(f" Current profile: {profile}")
375 print()
376
377 # Compute q profile
378 r = np.linspace(0, a, 200)
379 J = current_profile(r, a, I_p, profile=profile)
380 q = safety_factor_from_current(r, J, R0, B0, a)
381
382 q_0 = q[10] # Near axis (avoiding r=0 singularity)
383 q_edge = q[-1]
384
385 print(f"Safety factor:")
386 print(f" q(0) ≈ {q_0:.2f}")
387 print(f" q(a) = {q_edge:.2f}")
388 print()
389
390 # Check stability for different modes
391 modes = [(1, 1), (2, 1), (3, 1), (3, 2), (4, 1)]
392 print("Stability check (Kruskal-Shafranov criterion):")
393 print(f" {'Mode (m,n)':<12} {'q_crit':<8} {'Stable?'}")
394 print(" " + "-" * 35)
395
396 for m, n in modes:
397 stable = kruskal_shafranov_criterion(q_edge, m, n)
398 q_crit = m / n
399 status = "✓ STABLE" if stable else "✗ UNSTABLE"
400 print(f" ({m},{n}){' '*(10-len(f'({m},{n})'))} {q_crit:<8.2f} {status}")
401
402 print()
403
404 # Plot q profile
405 plot_q_profile(r, q, a, m_values=[1, 2, 3], n=1)
406 print("q profile plot saved as 'q_profile.png'")
407
408 # Plot stability boundary vs aspect ratio
409 plot_stability_diagram(R0, B0, profile=profile)
410 print("Stability boundary plot saved as 'stability_boundary.png'")
411
412 # Parameter scan
413 parameter_scan(R0, a, B0, profile=profile)
414 print("Marginal stability plot saved as 'marginal_stability.png'")
415
416
417if __name__ == '__main__':
418 main()