1#!/usr/bin/env python3
2"""
3Ballooning Mode Instability
4
5Solves the ballooning equation in s-α coordinates to find stability boundaries.
6
7s = r q'/q (magnetic shear)
8α = -2μ₀ Rq² p'/B² (normalized pressure gradient)
9
10Identifies first and second stability regions in the s-α diagram.
11
12Author: Claude
13Date: 2026-02-14
14"""
15
16import numpy as np
17import matplotlib.pyplot as plt
18from scipy.integrate import solve_ivp
19
20
21def ballooning_equation(theta, y, s, alpha):
22 """
23 Ballooning equation in ballooning representation.
24
25 d²X/dθ² + [s²θ² + α sin(θ)]X = 0
26
27 Parameters
28 ----------
29 theta : float
30 Poloidal angle
31 y : array
32 [X, dX/dtheta]
33 s : float
34 Magnetic shear
35 alpha : float
36 Pressure gradient parameter
37
38 Returns
39 -------
40 dydt theta : array
41 [dX/dθ, d²X/dθ²]
42 """
43 X, dX = y
44
45 # Ballooning equation coefficient
46 coeff = s**2 * theta**2 + alpha * np.sin(theta)
47
48 d2X = -coeff * X
49
50 return [dX, d2X]
51
52
53def solve_ballooning(s, alpha, theta_max=20):
54 """
55 Solve ballooning equation and determine stability.
56
57 Parameters
58 ----------
59 s : float
60 Magnetic shear
61 alpha : float
62 Pressure gradient
63 theta_max : float
64 Integration range
65
66 Returns
67 -------
68 stable : bool
69 True if stable
70 growth_indicator : float
71 Positive if unstable
72 """
73 # Initial conditions: X(0) = 1, dX/dθ(0) = 0
74 y0 = [1.0, 0.0]
75
76 # Integrate
77 theta_span = (0, theta_max)
78 theta_eval = np.linspace(0, theta_max, 500)
79
80 try:
81 sol = solve_ivp(ballooning_equation, theta_span, y0,
82 args=(s, alpha), t_eval=theta_eval,
83 method='RK45', max_step=0.1)
84
85 X = sol.y[0]
86
87 # Check for exponential growth
88 # If |X| grows exponentially, mode is unstable
89 growth_indicator = np.log(np.abs(X[-1]) / np.abs(X[0]) + 1e-10)
90
91 # Threshold for instability
92 stable = growth_indicator < 1.0
93
94 return stable, growth_indicator
95
96 except:
97 # Integration failed, likely unstable
98 return False, 10.0
99
100
101def create_s_alpha_diagram(s_range, alpha_range, resolution=50):
102 """
103 Create stability diagram in s-α space.
104
105 Parameters
106 ----------
107 s_range : tuple
108 (s_min, s_max)
109 alpha_range : tuple
110 (alpha_min, alpha_max)
111 resolution : int
112 Grid resolution
113
114 Returns
115 -------
116 S, Alpha : ndarray
117 Meshgrids
118 stability_map : ndarray
119 Stability map (1 = stable, 0 = unstable)
120 """
121 s_arr = np.linspace(s_range[0], s_range[1], resolution)
122 alpha_arr = np.linspace(alpha_range[0], alpha_range[1], resolution)
123
124 S, Alpha = np.meshgrid(s_arr, alpha_arr)
125 stability_map = np.zeros_like(S)
126
127 print("Computing s-α stability diagram...")
128 for i in range(len(alpha_arr)):
129 if i % 10 == 0:
130 print(f" Progress: {100*i/len(alpha_arr):.0f}%")
131
132 for j in range(len(s_arr)):
133 stable, _ = solve_ballooning(s_arr[j], alpha_arr[i])
134 stability_map[i, j] = 1.0 if stable else 0.0
135
136 print(" Complete!")
137 return S, Alpha, stability_map
138
139
140def plot_s_alpha_diagram(S, Alpha, stability_map):
141 """
142 Plot s-α stability diagram.
143
144 Parameters
145 ----------
146 S, Alpha : ndarray
147 Coordinate grids
148 stability_map : ndarray
149 Stability (1 = stable, 0 = unstable)
150 """
151 fig, ax = plt.subplots(figsize=(11, 8))
152
153 # Plot stability regions
154 im = ax.contourf(S, Alpha, stability_map, levels=[0, 0.5, 1.0],
155 colors=['red', 'lightgreen'], alpha=0.7)
156
157 # Stability boundary
158 cs = ax.contour(S, Alpha, stability_map, levels=[0.5],
159 colors='black', linewidths=3)
160
161 ax.set_xlabel('Magnetic shear s = r q\'/q', fontsize=13)
162 ax.set_ylabel('Pressure parameter α = -2μ₀Rq²p\'/B²', fontsize=13)
163 ax.set_title('Ballooning Stability Diagram',
164 fontsize=14, fontweight='bold')
165 ax.grid(True, alpha=0.3)
166
167 # Labels for regions
168 ax.text(0.5, 0.3, 'FIRST STABLE\nREGION',
169 fontsize=14, fontweight='bold', ha='center',
170 bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
171
172 # Check for second stability region at high s
173 if np.any((S > 2) & (stability_map > 0.5)):
174 ax.text(3.5, 2.0, 'SECOND STABLE\nREGION',
175 fontsize=14, fontweight='bold', ha='center',
176 bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
177
178 ax.text(1.0, 1.5, 'UNSTABLE',
179 fontsize=14, fontweight='bold', ha='center',
180 bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.8))
181
182 # Reference lines
183 ax.axhline(y=0, color='gray', linestyle='--', alpha=0.5)
184 ax.axvline(x=0, color='gray', linestyle='--', alpha=0.5)
185
186 plt.tight_layout()
187 plt.savefig('ballooning_s_alpha_diagram.png', dpi=150, bbox_inches='tight')
188 plt.show()
189
190
191def plot_eigenfunction_examples(s_values, alpha_test):
192 """
193 Plot eigenfunctions for different shear values.
194
195 Parameters
196 ----------
197 s_values : list
198 Shear values to plot
199 alpha_test : float
200 Pressure gradient parameter
201 """
202 theta_max = 15
203 theta = np.linspace(0, theta_max, 500)
204
205 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
206 axes = axes.flatten()
207
208 for i, s in enumerate(s_values[:4]):
209 ax = axes[i]
210
211 y0 = [1.0, 0.0]
212 sol = solve_ivp(ballooning_equation, (0, theta_max), y0,
213 args=(s, alpha_test), t_eval=theta,
214 method='RK45', max_step=0.1)
215
216 X = sol.y[0]
217 stable, growth = solve_ballooning(s, alpha_test)
218
219 ax.plot(theta, X, 'b-', linewidth=2)
220 ax.axhline(y=0, color='k', linestyle='--', alpha=0.3)
221 ax.set_xlabel('Poloidal angle θ', fontsize=11)
222 ax.set_ylabel('Eigenfunction X(θ)', fontsize=11)
223
224 status = "STABLE" if stable else "UNSTABLE"
225 color = "green" if stable else "red"
226
227 ax.set_title(f's={s:.2f}, α={alpha_test:.2f}\n{status}',
228 fontsize=12, fontweight='bold', color=color)
229 ax.grid(True, alpha=0.3)
230
231 plt.suptitle('Ballooning Mode Eigenfunctions',
232 fontsize=14, fontweight='bold')
233 plt.tight_layout()
234 plt.savefig('ballooning_eigenfunctions.png', dpi=150, bbox_inches='tight')
235 plt.show()
236
237
238def plot_stability_boundaries():
239 """
240 Plot stability boundaries for different alpha values.
241 """
242 s_arr = np.linspace(0, 5, 100)
243 alpha_values = [0.5, 1.0, 1.5, 2.0, 2.5]
244
245 fig, ax = plt.subplots(figsize=(10, 7))
246
247 colors = plt.cm.viridis(np.linspace(0, 1, len(alpha_values)))
248
249 for i, alpha in enumerate(alpha_values):
250 growth_arr = []
251
252 for s in s_arr:
253 stable, growth = solve_ballooning(s, alpha)
254 growth_arr.append(growth)
255
256 ax.plot(s_arr, growth_arr, color=colors[i], linewidth=2.5,
257 label=f'α = {alpha:.1f}')
258
259 ax.axhline(y=1.0, color='red', linestyle='--', linewidth=2,
260 label='Stability threshold')
261 ax.set_xlabel('Magnetic shear s', fontsize=13)
262 ax.set_ylabel('Growth indicator', fontsize=13)
263 ax.set_title('Ballooning Stability vs Magnetic Shear',
264 fontsize=14, fontweight='bold')
265 ax.legend(fontsize=10)
266 ax.grid(True, alpha=0.3)
267 ax.set_ylim([-0.5, 5])
268
269 # Shade stable region
270 ax.fill_between(s_arr, -0.5, 1.0, alpha=0.2, color='green',
271 label='Stable region')
272
273 plt.tight_layout()
274 plt.savefig('ballooning_vs_shear.png', dpi=150, bbox_inches='tight')
275 plt.show()
276
277
278def main():
279 """Main execution function."""
280 print("Ballooning Mode Stability Analysis")
281 print("=" * 60)
282
283 # Create s-α diagram
284 print("\nCreating stability diagram in s-α space...")
285 s_range = (0, 5)
286 alpha_range = (0, 3)
287 S, Alpha, stability_map = create_s_alpha_diagram(s_range, alpha_range,
288 resolution=40)
289
290 # Find stability boundaries
291 n_stable = np.sum(stability_map > 0.5)
292 n_total = stability_map.size
293 print(f"\nStable points: {n_stable}/{n_total} ({100*n_stable/n_total:.1f}%)")
294
295 # Plot s-α diagram
296 plot_s_alpha_diagram(S, Alpha, stability_map)
297 print("s-α diagram saved as 'ballooning_s_alpha_diagram.png'")
298
299 # Plot eigenfunctions
300 print("\nComputing example eigenfunctions...")
301 s_test = [0.5, 1.0, 2.0, 3.0]
302 alpha_test = 1.0
303 plot_eigenfunction_examples(s_test, alpha_test)
304 print("Eigenfunctions saved as 'ballooning_eigenfunctions.png'")
305
306 # Plot stability vs shear
307 print("\nPlotting stability vs magnetic shear...")
308 plot_stability_boundaries()
309 print("Shear dependence saved as 'ballooning_vs_shear.png'")
310
311 print("\nAnalysis complete!")
312
313
314if __name__ == '__main__':
315 main()