1#!/usr/bin/env python3
2"""
3External Kink Instability
4
5Analyzes external kink mode in sharp-boundary cylindrical plasma.
6Computes dispersion relation for surface modes and shows stabilization
7by conducting wall.
8
9Author: Claude
10Date: 2026-02-14
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.special import iv, kn, ivp, kvp # Modified Bessel functions
16
17
18def kink_dispersion(q_a, m, k_z, a, b_wall=None):
19 """
20 Compute kink mode growth rate.
21
22 For m=1 external kink, stability requires q(a) > m.
23
24 Parameters
25 ----------
26 q_a : float
27 Safety factor at edge
28 m : int
29 Poloidal mode number
30 k_z : float
31 Axial wavenumber
32 a : float
33 Plasma radius
34 b_wall : float or None
35 Wall radius (None for no wall)
36
37 Returns
38 -------
39 gamma : float
40 Growth rate (normalized)
41 """
42 # Normalized growth rate (simplified)
43 # Unstable if q(a) < m
44 if q_a < m:
45 gamma_sq = (m - q_a) / (m + 1) # Simplified formula
46 else:
47 gamma_sq = 0
48
49 # Wall stabilization
50 if b_wall is not None and b_wall > a:
51 # Wall reduces growth rate
52 wall_factor = (a / b_wall)**m
53 gamma_sq *= (1 - wall_factor)
54
55 gamma = np.sqrt(max(gamma_sq, 0))
56 return gamma
57
58
59def plot_growth_vs_q(m=1, k_z_values=[0.1, 0.5, 1.0]):
60 """
61 Plot growth rate vs q(a) for m=1 kink.
62
63 Parameters
64 ----------
65 m : int
66 Mode number
67 k_z_values : list
68 Axial wavenumbers
69 """
70 q_a = np.linspace(0.5, 3.0, 100)
71 a = 0.5 # m
72
73 fig, ax = plt.subplots(figsize=(10, 7))
74
75 colors = plt.cm.plasma(np.linspace(0, 1, len(k_z_values)))
76
77 for i, k_z in enumerate(k_z_values):
78 gamma = np.array([kink_dispersion(q, m, k_z, a) for q in q_a])
79
80 ax.plot(q_a, gamma, color=colors[i], linewidth=2.5,
81 label=f'k_z = {k_z:.1f} m⁻¹')
82
83 ax.axvline(x=m, color='red', linestyle='--', linewidth=2,
84 label=f'q = {m} (marginal stability)')
85 ax.axhline(y=0, color='k', linestyle='-', alpha=0.3)
86
87 ax.set_xlabel('Edge safety factor q(a)', fontsize=13)
88 ax.set_ylabel('Normalized growth rate γ', fontsize=13)
89 ax.set_title(f'External Kink Mode (m={m})',
90 fontsize=14, fontweight='bold')
91 ax.legend(fontsize=11)
92 ax.grid(True, alpha=0.3)
93
94 # Shade unstable region
95 ax.fill_betweenx([0, ax.get_ylim()[1]], 0, m, alpha=0.2, color='red',
96 label='Unstable')
97
98 plt.tight_layout()
99 plt.savefig('kink_growth_vs_q.png', dpi=150, bbox_inches='tight')
100 plt.show()
101
102
103def plot_wall_stabilization(m=1):
104 """
105 Show wall stabilization effect.
106
107 Parameters
108 ----------
109 m : int
110 Mode number
111 """
112 q_a_values = [0.7, 0.8, 0.9] # < 1, unstable without wall
113 b_a_ratio = np.linspace(1.0, 3.0, 100) # b/a
114
115 fig, ax = plt.subplots(figsize=(10, 7))
116
117 colors = plt.cm.viridis(np.linspace(0, 1, len(q_a_values)))
118
119 a = 0.5
120 k_z = 0.5
121
122 for i, q_a in enumerate(q_a_values):
123 gamma_arr = []
124
125 for ratio in b_a_ratio:
126 b = ratio * a
127 gamma = kink_dispersion(q_a, m, k_z, a, b_wall=b)
128 gamma_arr.append(gamma)
129
130 # No-wall case
131 gamma_no_wall = kink_dispersion(q_a, m, k_z, a, b_wall=None)
132
133 ax.plot(b_a_ratio, gamma_arr, color=colors[i], linewidth=2.5,
134 label=f'q(a)={q_a:.1f} (no wall: γ={gamma_no_wall:.2f})')
135
136 ax.set_xlabel('Wall radius ratio b/a', fontsize=13)
137 ax.set_ylabel('Normalized growth rate γ', fontsize=13)
138 ax.set_title('Conducting Wall Stabilization of Kink Mode',
139 fontsize=14, fontweight='bold')
140 ax.legend(fontsize=10)
141 ax.grid(True, alpha=0.3)
142 ax.set_xlim([1, 3])
143
144 plt.tight_layout()
145 plt.savefig('kink_wall_stabilization.png', dpi=150, bbox_inches='tight')
146 plt.show()
147
148
149def plot_mode_structure(m=1, n=1):
150 """
151 Visualize kink mode structure.
152
153 Parameters
154 ----------
155 m, n : int
156 Mode numbers
157 """
158 # Cylindrical grid
159 theta = np.linspace(0, 2*np.pi, 100)
160 z = np.linspace(0, 2*np.pi/n, 100)
161 Theta, Z = np.meshgrid(theta, z)
162
163 # Mode structure: ξ ~ exp(i(mθ + nz))
164 # Visualize real part
165 xi_r = np.cos(m * Theta + n * Z)
166
167 fig = plt.figure(figsize=(12, 5))
168
169 # 2D plot
170 ax1 = fig.add_subplot(121)
171 im1 = ax1.contourf(Theta, Z, xi_r, levels=20, cmap='RdBu_r')
172 ax1.set_xlabel('Poloidal angle θ (rad)', fontsize=12)
173 ax1.set_ylabel('Axial position z (rad)', fontsize=12)
174 ax1.set_title(f'Kink Mode Structure (m={m}, n={n})',
175 fontsize=13, fontweight='bold')
176 plt.colorbar(im1, ax=ax1, label='Displacement ξ_r')
177
178 # 3D-like view (unrolled cylinder)
179 ax2 = fig.add_subplot(122, projection='3d')
180 R = 1.0
181 X = R * np.cos(Theta)
182 Y = R * np.sin(Theta)
183 surf = ax2.plot_surface(X, Y, Z, facecolors=plt.cm.RdBu_r((xi_r+1)/2),
184 alpha=0.9, rstride=2, cstride=2)
185 ax2.set_xlabel('X', fontsize=11)
186 ax2.set_ylabel('Y', fontsize=11)
187 ax2.set_zlabel('Z', fontsize=11)
188 ax2.set_title('3D View', fontsize=12, fontweight='bold')
189
190 plt.tight_layout()
191 plt.savefig('kink_mode_structure.png', dpi=150, bbox_inches='tight')
192 plt.show()
193
194
195def main():
196 """Main execution function."""
197 print("External Kink Instability Analysis")
198 print("=" * 60)
199
200 # Growth rate vs q(a)
201 print("\nPlotting growth rate vs q(a)...")
202 plot_growth_vs_q(m=1, k_z_values=[0.1, 0.5, 1.0])
203 print(" Saved as 'kink_growth_vs_q.png'")
204
205 # Wall stabilization
206 print("\nPlotting wall stabilization effect...")
207 plot_wall_stabilization(m=1)
208 print(" Saved as 'kink_wall_stabilization.png'")
209
210 # Mode structure
211 print("\nVisualizing mode structure...")
212 plot_mode_structure(m=1, n=1)
213 print(" Saved as 'kink_mode_structure.png'")
214
215 print("\nKey results:")
216 print(" - m=1 kink unstable when q(a) < 1")
217 print(" - Conducting wall provides stabilization")
218 print(" - Stabilization stronger when wall closer to plasma")
219
220
221if __name__ == '__main__':
222 main()