1#!/usr/bin/env python3
2"""
3Banana Orbit Simulation in Tokamak
4
5This script simulates banana orbits in a simplified tokamak geometry.
6Demonstrates the difference between passing and trapped particles in
7toroidal magnetic field configuration.
8
9Author: Plasma Physics Examples
10License: MIT
11"""
12
13import numpy as np
14import matplotlib.pyplot as plt
15from scipy.constants import e, m_e
16
17
18def boris_push(x, v, q, m, B, dt):
19 """Boris algorithm for particle pushing (no E field)."""
20 t = (q * dt / (2 * m)) * B
21 t_mag2 = np.dot(t, t)
22 s = 2 * t / (1 + t_mag2)
23
24 v_prime = v + np.cross(v, t)
25 v_new = v + np.cross(v_prime, s)
26 x_new = x + v_new * dt
27
28 return x_new, v_new
29
30
31def tokamak_field(R, Z, phi, R0, B0, epsilon):
32 """
33 Simplified tokamak magnetic field using large aspect ratio approximation.
34
35 B_toroidal = B0 * R0 / R (1/R dependence)
36 B_poloidal ≈ ε * B0 (simplified, constant)
37
38 Parameters:
39 -----------
40 R, Z, phi : float
41 Cylindrical coordinates [m, m, rad]
42 R0 : float
43 Major radius [m]
44 B0 : float
45 Magnetic field at R0 [T]
46 epsilon : float
47 Inverse aspect ratio (a/R0), determines poloidal field strength
48
49 Returns:
50 --------
51 B : array (3,)
52 Magnetic field in cylindrical coordinates [B_R, B_Z, B_phi]
53 """
54 # Toroidal field (dominant)
55 B_phi = B0 * R0 / R
56
57 # Poloidal field (simplified as constant, pointing in Z direction)
58 # In reality, B_poloidal varies with position
59 B_R = 0.0
60 B_Z = epsilon * B0
61
62 return np.array([B_R, B_Z, B_phi])
63
64
65def cylindrical_to_cartesian(R, Z, phi):
66 """Convert cylindrical to Cartesian coordinates."""
67 x = R * np.cos(phi)
68 y = R * np.sin(phi)
69 z = Z
70 return np.array([x, y, z])
71
72
73def cartesian_to_cylindrical(x, y, z):
74 """Convert Cartesian to cylindrical coordinates."""
75 R = np.sqrt(x**2 + y**2)
76 Z = z
77 phi = np.arctan2(y, x)
78 return R, Z, phi
79
80
81def cyl_field_to_cart(B_cyl, phi):
82 """
83 Convert magnetic field from cylindrical to Cartesian coordinates.
84
85 B_cyl = [B_R, B_Z, B_phi]
86 """
87 B_R, B_Z, B_phi = B_cyl
88
89 B_x = B_R * np.cos(phi) - B_phi * np.sin(phi)
90 B_y = B_R * np.sin(phi) + B_phi * np.cos(phi)
91 B_z = B_Z
92
93 return np.array([B_x, B_y, B_z])
94
95
96def simulate_tokamak_orbit(q, m, R0, B0, epsilon, v0_cyl, x0_cyl, dt, n_steps):
97 """
98 Simulate particle orbit in tokamak.
99
100 Parameters:
101 -----------
102 q, m : float
103 Charge and mass
104 R0, B0 : float
105 Tokamak major radius and field
106 epsilon : float
107 Inverse aspect ratio
108 v0_cyl : array (3,)
109 Initial velocity in cylindrical coords [v_R, v_Z, v_phi]
110 x0_cyl : array (3,)
111 Initial position in cylindrical coords [R, Z, phi]
112 dt : float
113 Timestep
114 n_steps : int
115 Number of steps
116
117 Returns:
118 --------
119 trajectory : dict
120 """
121 # Initialize in Cartesian coordinates
122 x = np.zeros((n_steps, 3))
123 v = np.zeros((n_steps, 3))
124 t = np.zeros(n_steps)
125
126 R_traj = np.zeros(n_steps)
127 Z_traj = np.zeros(n_steps)
128 phi_traj = np.zeros(n_steps)
129
130 # Initial conditions
131 R0_pos, Z0_pos, phi0 = x0_cyl
132 x[0] = cylindrical_to_cartesian(R0_pos, Z0_pos, phi0)
133
134 # Convert initial velocity to Cartesian
135 v_R, v_Z, v_phi = v0_cyl
136 v[0] = np.array([
137 v_R * np.cos(phi0) - v_phi * np.sin(phi0),
138 v_R * np.sin(phi0) + v_phi * np.cos(phi0),
139 v_Z
140 ])
141
142 R_traj[0] = R0_pos
143 Z_traj[0] = Z0_pos
144 phi_traj[0] = phi0
145
146 for i in range(n_steps - 1):
147 # Get current cylindrical coordinates
148 R, Z, phi = cartesian_to_cylindrical(x[i, 0], x[i, 1], x[i, 2])
149
150 # Get magnetic field in cylindrical coords
151 B_cyl = tokamak_field(R, Z, phi, R0, B0, epsilon)
152
153 # Convert to Cartesian
154 B_cart = cyl_field_to_cart(B_cyl, phi)
155
156 # Push particle
157 x[i+1], v[i+1] = boris_push(x[i], v[i], q, m, B_cart, dt)
158 t[i+1] = t[i] + dt
159
160 # Store cylindrical coordinates
161 R_traj[i+1], Z_traj[i+1], phi_traj[i+1] = cartesian_to_cylindrical(
162 x[i+1, 0], x[i+1, 1], x[i+1, 2]
163 )
164
165 return {
166 'x': x[:, 0], 'y': x[:, 1], 'z': x[:, 2],
167 'vx': v[:, 0], 'vy': v[:, 1], 'vz': v[:, 2],
168 'R': R_traj, 'Z': Z_traj, 'phi': phi_traj,
169 't': t
170 }
171
172
173def plot_banana_orbit():
174 """Demonstrate banana orbit (trapped particle)."""
175
176 # Tokamak parameters
177 R0 = 1.0 # Major radius [m]
178 B0 = 2.0 # Magnetic field [T]
179 epsilon = 0.3 # Inverse aspect ratio (a/R0)
180 a = epsilon * R0 # Minor radius
181
182 print(f"\nTokamak Parameters:")
183 print(f" Major radius R0 = {R0} m")
184 print(f" Minor radius a = {a} m")
185 print(f" Magnetic field B0 = {B0} T")
186 print(f" Inverse aspect ratio ε = {epsilon}")
187
188 # Particle parameters
189 q = -e
190 m = m_e
191
192 # Trapped particle: start at outer midplane with velocity mostly poloidal
193 R_start = R0 + 0.5 * a # Start at outer part of torus
194 Z_start = 0.0
195 phi_start = 0.0
196
197 # Velocity: small parallel, large perpendicular (trapped condition)
198 v_thermal = 1e7 # m/s
199 v_parallel = 0.3 * v_thermal # Small parallel velocity
200 v_perp = np.sqrt(v_thermal**2 - v_parallel**2)
201
202 # Initial velocity in cylindrical coords
203 v0_cyl = np.array([0.0, v_perp, v_parallel])
204 x0_cyl = np.array([R_start, Z_start, phi_start])
205
206 # Time parameters
207 omega_c = abs(q) * B0 / m
208 T_c = 2 * np.pi / omega_c
209 dt = T_c / 50
210 n_steps = 5000
211
212 # Simulate
213 traj = simulate_tokamak_orbit(q, m, R0, B0, epsilon, v0_cyl, x0_cyl, dt, n_steps)
214
215 # Create figure
216 fig = plt.figure(figsize=(16, 10))
217
218 # Plot 1: 3D trajectory
219 ax1 = fig.add_subplot(2, 3, 1, projection='3d')
220 ax1.plot(traj['x'], traj['y'], traj['z'], 'b-', linewidth=1, alpha=0.7)
221 ax1.plot([traj['x'][0]], [traj['y'][0]], [traj['z'][0]],
222 'go', markersize=10, label='Start')
223
224 # Draw torus outline
225 theta_tor = np.linspace(0, 2*np.pi, 100)
226 for phi in [0, np.pi/2, np.pi, 3*np.pi/2]:
227 x_tor = (R0 + a * np.cos(theta_tor)) * np.cos(phi)
228 y_tor = (R0 + a * np.cos(theta_tor)) * np.sin(phi)
229 z_tor = a * np.sin(theta_tor)
230 ax1.plot(x_tor, y_tor, z_tor, 'k--', alpha=0.3, linewidth=1)
231
232 ax1.set_xlabel('x [m]', fontsize=11)
233 ax1.set_ylabel('y [m]', fontsize=11)
234 ax1.set_zlabel('z [m]', fontsize=11)
235 ax1.set_title('3D Banana Orbit', fontsize=13, fontweight='bold')
236 ax1.legend(loc='best', fontsize=9)
237
238 # Plot 2: Poloidal cross-section (R-Z plane)
239 ax2 = fig.add_subplot(2, 3, 2)
240 ax2.plot(traj['R'], traj['Z'], 'b-', linewidth=2, alpha=0.7)
241 ax2.plot(traj['R'][0], traj['Z'][0], 'go', markersize=10, label='Start')
242
243 # Draw plasma boundary (circular cross-section)
244 theta = np.linspace(0, 2*np.pi, 100)
245 R_boundary = R0 + a * np.cos(theta)
246 Z_boundary = a * np.sin(theta)
247 ax2.plot(R_boundary, Z_boundary, 'k--', linewidth=2, label='Plasma boundary')
248 ax2.plot(R0, 0, 'r*', markersize=15, label='Magnetic axis')
249
250 ax2.set_xlabel('R [m]', fontsize=12)
251 ax2.set_ylabel('Z [m]', fontsize=12)
252 ax2.set_title('Poloidal Cross-Section (Banana Shape)', fontsize=13, fontweight='bold')
253 ax2.legend(loc='best', fontsize=10)
254 ax2.grid(True, alpha=0.3)
255 ax2.set_aspect('equal')
256
257 # Plot 3: Toroidal angle evolution
258 ax3 = fig.add_subplot(2, 3, 3)
259 # Unwrap phi to avoid discontinuities
260 phi_unwrapped = np.unwrap(traj['phi'])
261 ax3.plot(traj['t'] * 1e6, phi_unwrapped * 180/np.pi, 'b-', linewidth=2)
262
263 ax3.set_xlabel('Time [μs]', fontsize=12)
264 ax3.set_ylabel('Toroidal Angle φ [degrees]', fontsize=12)
265 ax3.set_title('Toroidal Motion', fontsize=13, fontweight='bold')
266 ax3.grid(True, alpha=0.3)
267
268 # Plot 4: R vs time (shows banana tips)
269 ax4 = fig.add_subplot(2, 3, 4)
270 ax4.plot(traj['t'] * 1e6, traj['R'], 'b-', linewidth=2)
271 ax4.axhline(R0 + a, color='red', linestyle='--', linewidth=2, label='Outer boundary')
272 ax4.axhline(R0 - a, color='red', linestyle='--', linewidth=2, label='Inner boundary')
273 ax4.axhline(R0, color='green', linestyle=':', linewidth=2, label='Magnetic axis')
274
275 ax4.set_xlabel('Time [μs]', fontsize=12)
276 ax4.set_ylabel('Major Radius R [m]', fontsize=12)
277 ax4.set_title('Radial Excursion', fontsize=13, fontweight='bold')
278 ax4.legend(loc='best', fontsize=9)
279 ax4.grid(True, alpha=0.3)
280
281 # Plot 5: Z vs time (poloidal oscillation)
282 ax5 = fig.add_subplot(2, 3, 5)
283 ax5.plot(traj['t'] * 1e6, traj['Z'], 'b-', linewidth=2)
284 ax5.axhline(a, color='red', linestyle='--', linewidth=2, label='Upper boundary')
285 ax5.axhline(-a, color='red', linestyle='--', linewidth=2, label='Lower boundary')
286
287 ax5.set_xlabel('Time [μs]', fontsize=12)
288 ax5.set_ylabel('Z [m]', fontsize=12)
289 ax5.set_title('Vertical Motion', fontsize=13, fontweight='bold')
290 ax5.legend(loc='best', fontsize=9)
291 ax5.grid(True, alpha=0.3)
292
293 # Plot 6: Banana width
294 ax6 = fig.add_subplot(2, 3, 6)
295 # Calculate banana width from R excursion
296 R_max = np.max(traj['R'])
297 R_min = np.min(traj['R'])
298 banana_width = R_max - R_min
299
300 ax6.hist(traj['R'], bins=50, color='blue', alpha=0.7, edgecolor='black')
301 ax6.axvline(R_max, color='red', linestyle='--', linewidth=2,
302 label=f'Width = {banana_width*100:.2f} cm')
303 ax6.axvline(R_min, color='red', linestyle='--', linewidth=2)
304
305 ax6.set_xlabel('R [m]', fontsize=12)
306 ax6.set_ylabel('Frequency', fontsize=12)
307 ax6.set_title('Banana Width Distribution', fontsize=13, fontweight='bold')
308 ax6.legend(loc='best', fontsize=10)
309 ax6.grid(True, alpha=0.3)
310
311 print(f"\nBanana Orbit Characteristics:")
312 print(f" Banana width: {banana_width*100:.2f} cm")
313 print(f" R range: [{R_min:.3f}, {R_max:.3f}] m")
314 print(f" Z range: [{np.min(traj['Z']):.3f}, {np.max(traj['Z']):.3f}] m")
315
316 plt.tight_layout()
317 plt.savefig('banana_orbit.png', dpi=300, bbox_inches='tight')
318 plt.show()
319
320
321def plot_passing_vs_trapped():
322 """Compare passing and trapped orbits."""
323
324 # Tokamak parameters
325 R0 = 1.0
326 B0 = 2.0
327 epsilon = 0.3
328 a = epsilon * R0
329
330 # Particle parameters
331 q = -e
332 m = m_e
333
334 # Initial position
335 R_start = R0 + 0.5 * a
336 Z_start = 0.0
337 phi_start = 0.0
338
339 v_thermal = 1e7 # m/s
340
341 # Two cases: trapped vs passing
342 # Trapped: small v_parallel (< critical)
343 # Passing: large v_parallel (> critical)
344
345 cases = [
346 ('Trapped', 0.2 * v_thermal), # Small parallel velocity
347 ('Passing', 0.8 * v_thermal), # Large parallel velocity
348 ]
349
350 omega_c = abs(q) * B0 / m
351 T_c = 2 * np.pi / omega_c
352 dt = T_c / 50
353 n_steps = 5000
354
355 fig, axes = plt.subplots(1, 2, figsize=(14, 6))
356
357 colors = ['blue', 'red']
358
359 for i, (label, v_para) in enumerate(cases):
360 v_perp = np.sqrt(v_thermal**2 - v_para**2)
361 v0_cyl = np.array([0.0, v_perp, v_para])
362 x0_cyl = np.array([R_start, Z_start, phi_start])
363
364 traj = simulate_tokamak_orbit(q, m, R0, B0, epsilon, v0_cyl, x0_cyl, dt, n_steps)
365
366 # Plot in R-Z plane
367 ax = axes[0]
368 ax.plot(traj['R'], traj['Z'], color=colors[i], linewidth=2,
369 alpha=0.7, label=f'{label} (v_∥/v_th = {v_para/v_thermal:.1f})')
370
371 # Draw plasma boundary
372 ax = axes[0]
373 theta = np.linspace(0, 2*np.pi, 100)
374 R_boundary = R0 + a * np.cos(theta)
375 Z_boundary = a * np.sin(theta)
376 ax.plot(R_boundary, Z_boundary, 'k--', linewidth=2, label='Plasma boundary')
377 ax.plot(R0, 0, 'g*', markersize=15, label='Magnetic axis')
378
379 ax.set_xlabel('R [m]', fontsize=12)
380 ax.set_ylabel('Z [m]', fontsize=12)
381 ax.set_title('Passing vs Trapped Orbits', fontsize=13, fontweight='bold')
382 ax.legend(loc='best', fontsize=10)
383 ax.grid(True, alpha=0.3)
384 ax.set_aspect('equal')
385
386 # Second plot: phi vs time
387 ax = axes[1]
388 for i, (label, v_para) in enumerate(cases):
389 v_perp = np.sqrt(v_thermal**2 - v_para**2)
390 v0_cyl = np.array([0.0, v_perp, v_para])
391 x0_cyl = np.array([R_start, Z_start, phi_start])
392
393 traj = simulate_tokamak_orbit(q, m, R0, B0, epsilon, v0_cyl, x0_cyl, dt, n_steps)
394
395 phi_unwrapped = np.unwrap(traj['phi'])
396 ax.plot(traj['t'] * 1e6, phi_unwrapped, color=colors[i],
397 linewidth=2, label=label)
398
399 ax.set_xlabel('Time [μs]', fontsize=12)
400 ax.set_ylabel('Toroidal Angle φ [rad]', fontsize=12)
401 ax.set_title('Toroidal Circulation', fontsize=13, fontweight='bold')
402 ax.legend(loc='best', fontsize=10)
403 ax.grid(True, alpha=0.3)
404
405 plt.tight_layout()
406 plt.savefig('passing_vs_trapped.png', dpi=300, bbox_inches='tight')
407 plt.show()
408
409
410if __name__ == '__main__':
411 print("\n" + "="*80)
412 print("BANANA ORBIT SIMULATION IN TOKAMAK")
413 print("="*80)
414
415 print("\nPart 1: Banana Orbit (Trapped Particle)")
416 print("-" * 80)
417 plot_banana_orbit()
418
419 print("\nPart 2: Passing vs Trapped Comparison")
420 print("-" * 80)
421 plot_passing_vs_trapped()
422
423 print("\nKey Points:")
424 print(" - Trapped particles: small v_∥, banana-shaped orbits")
425 print(" - Passing particles: large v_∥, circulate around torus")
426 print(" - Banana width depends on v_∥ and radial position")
427
428 print("\nDone! Generated 2 figures:")
429 print(" - banana_orbit.png")
430 print(" - passing_vs_trapped.png")