fusion_core/
rf_heating.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SCPN Fusion Core — RF Heating
7//! ICRH ray tracing for RF heating deposition.
8//!
9//! Port of `rf_heating.py` (200 lines).
10//! Uses Hamiltonian ray equations in cold plasma approximation.
11
12use fusion_types::constants::MU0_SI;
13use fusion_types::error::{FusionError, FusionResult};
14
15/// Deuterium charge [C]. Python line 23.
16const Q_D: f64 = 1.602e-19;
17
18/// Deuterium mass [kg]. Python line 24.
19const M_D: f64 = 3.34e-27;
20
21/// ICRH frequency [Hz]. Python line 25.
22const FREQ_HZ: f64 = 50e6;
23
24/// Angular frequency [rad/s].
25const OMEGA_WAVE: f64 = 2.0 * std::f64::consts::PI * FREQ_HZ;
26
27/// Magnetic field on axis [T]. Python line 34.
28const B0: f64 = 5.3;
29
30/// Major radius [m]. Python line 35.
31const R0: f64 = 6.2;
32
33/// Peak electron density [m⁻³]. Python line 56.
34const N_E_PEAK: f64 = 1.0e20;
35
36/// Density Gaussian width [m²]. Python line 56.
37const DENSITY_WIDTH: f64 = 2.0;
38
39/// Finite difference step for gradient [m]. Python line 96.
40const GRAD_EPS: f64 = 1e-3;
41
42/// Initial wavenumber [m⁻¹]. Python line 139.
43const K0_INITIAL: f64 = 10.0;
44
45/// Number of rays to trace. Python line 131.
46const N_RAYS: usize = 10;
47
48/// Antenna radial position [m]. Python line 131.
49const R_ANTENNA: f64 = 9.0;
50
51/// Resonance layer intersection tolerance [m]. Python line 179.
52const RESONANCE_TOL: f64 = 0.1;
53
54/// Ray state: (R, Z, k_R, k_Z).
55pub type RayState = [f64; 4];
56
57/// Result of a single ray trace.
58#[derive(Debug, Clone)]
59pub struct RayTraceResult {
60    /// Ray trajectory points (R, Z).
61    pub trajectory: Vec<(f64, f64)>,
62    /// Resonance crossing position, if found.
63    pub resonance_point: Option<(f64, f64)>,
64    /// Resonance layer R coordinate.
65    pub r_resonance: f64,
66}
67
68/// ICRH ray tracing system.
69pub struct RFHeatingSystem {
70    pub omega: f64,
71    pub b0: f64,
72    pub r0: f64,
73}
74
75impl RFHeatingSystem {
76    /// Create with ITER-like parameters.
77    pub fn new() -> Self {
78        RFHeatingSystem {
79            omega: OMEGA_WAVE,
80            b0: B0,
81            r0: R0,
82        }
83    }
84
85    fn validate_runtime_parameters(&self) -> FusionResult<()> {
86        if !self.omega.is_finite() || self.omega <= 0.0 {
87            return Err(FusionError::ConfigError(format!(
88                "RF omega must be finite and > 0, got {}",
89                self.omega
90            )));
91        }
92        if !self.b0.is_finite() || self.b0 <= 0.0 {
93            return Err(FusionError::ConfigError(format!(
94                "RF b0 must be finite and > 0, got {}",
95                self.b0
96            )));
97        }
98        if !self.r0.is_finite() || self.r0 <= 0.0 {
99            return Err(FusionError::ConfigError(format!(
100                "RF r0 must be finite and > 0, got {}",
101                self.r0
102            )));
103        }
104        Ok(())
105    }
106
107    /// Toroidal magnetic field at radius R: B(R) = B0 · R0 / R.
108    fn b_field(&self, r: f64) -> f64 {
109        self.b0 * self.r0 / r
110    }
111
112    /// Electron density profile (Gaussian). Python line 56.
113    fn density(&self, r: f64, z: f64) -> f64 {
114        let dist_sq = (r - self.r0).powi(2) + z.powi(2);
115        N_E_PEAK * (-dist_sq / DENSITY_WIDTH).exp()
116    }
117
118    /// Alfvén speed: v_A = B / sqrt(μ₀ n_e m_D).
119    fn alfven_speed(&self, r: f64, z: f64) -> f64 {
120        let b = self.b_field(r);
121        let n = self.density(r, z).max(1e10); // floor to avoid division by zero
122        b / (MU0_SI * n * M_D).sqrt()
123    }
124
125    /// Dispersion relation: D = k²v_A² - ω².
126    fn dispersion(&self, r: f64, z: f64, kr: f64, kz: f64) -> f64 {
127        let k_sq = kr * kr + kz * kz;
128        let v_a = self.alfven_speed(r, z);
129        k_sq * v_a * v_a - self.omega * self.omega
130    }
131
132    /// Ray equations (Hamiltonian): ds/dt = (dR/dt, dZ/dt, dk_R/dt, dk_Z/dt).
133    ///
134    /// dR/dt = -∂D/∂k_R,  dZ/dt = -∂D/∂k_Z
135    /// dk_R/dt = ∂D/∂R,   dk_Z/dt = ∂D/∂Z
136    fn ray_rhs(&self, state: &RayState) -> RayState {
137        let [r, z, kr, kz] = *state;
138        let eps = GRAD_EPS;
139
140        // ∂D/∂k_R, ∂D/∂k_Z via finite differences
141        let dd_dkr = (self.dispersion(r, z, kr + eps, kz) - self.dispersion(r, z, kr - eps, kz))
142            / (2.0 * eps);
143        let dd_dkz = (self.dispersion(r, z, kr, kz + eps) - self.dispersion(r, z, kr, kz - eps))
144            / (2.0 * eps);
145
146        // ∂D/∂R, ∂D/∂Z via finite differences
147        let dd_dr = (self.dispersion(r + eps, z, kr, kz) - self.dispersion(r - eps, z, kr, kz))
148            / (2.0 * eps);
149        let dd_dz = (self.dispersion(r, z + eps, kr, kz) - self.dispersion(r, z - eps, kr, kz))
150            / (2.0 * eps);
151
152        [
153            -dd_dkr, // dR/dt
154            -dd_dkz, // dZ/dt
155            dd_dr,   // dk_R/dt
156            dd_dz,   // dk_Z/dt
157        ]
158    }
159
160    /// Adaptive RK4 integration step with sub-stepping.
161    /// The Alfvén speed can be ~1e7 m/s, so group velocity dR/dt ~ k·v_A²
162    /// can be extremely large. We use adaptive sub-stepping to keep
163    /// displacement per step < 0.1 m.
164    fn rk4_step(&self, state: &RayState, dt: f64) -> RayState {
165        let k1 = self.ray_rhs(state);
166
167        let s2 = [
168            state[0] + 0.5 * dt * k1[0],
169            state[1] + 0.5 * dt * k1[1],
170            state[2] + 0.5 * dt * k1[2],
171            state[3] + 0.5 * dt * k1[3],
172        ];
173        let k2 = self.ray_rhs(&s2);
174
175        let s3 = [
176            state[0] + 0.5 * dt * k2[0],
177            state[1] + 0.5 * dt * k2[1],
178            state[2] + 0.5 * dt * k2[2],
179            state[3] + 0.5 * dt * k2[3],
180        ];
181        let k3 = self.ray_rhs(&s3);
182
183        let s4 = [
184            state[0] + dt * k3[0],
185            state[1] + dt * k3[1],
186            state[2] + dt * k3[2],
187            state[3] + dt * k3[3],
188        ];
189        let k4 = self.ray_rhs(&s4);
190
191        [
192            state[0] + (dt / 6.0) * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]),
193            state[1] + (dt / 6.0) * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]),
194            state[2] + (dt / 6.0) * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]),
195            state[3] + (dt / 6.0) * (k1[3] + 2.0 * k2[3] + 2.0 * k3[3] + k4[3]),
196        ]
197    }
198
199    /// Ion cyclotron resonance radius.
200    ///
201    /// ω_ci = qB/m → B_res = ωm/q → R_res = B0·R0/B_res
202    pub fn resonance_radius(&self) -> FusionResult<f64> {
203        self.validate_runtime_parameters()?;
204        let b_res = self.omega * M_D / Q_D;
205        if !b_res.is_finite() || b_res <= 0.0 {
206            return Err(FusionError::ConfigError(
207                "RF resonance field must be finite and > 0".to_string(),
208            ));
209        }
210        let r_res = self.b0 * self.r0 / b_res;
211        if !r_res.is_finite() || r_res <= 0.0 {
212            return Err(FusionError::ConfigError(
213                "RF resonance radius became non-finite or non-positive".to_string(),
214            ));
215        }
216        Ok(r_res)
217    }
218
219    /// Trace a single ray from given initial conditions.
220    ///
221    /// Uses adaptive time stepping: estimates RHS magnitude and scales dt
222    /// so the spatial displacement per step is bounded.
223    pub fn trace_single_ray(
224        &self,
225        r0: f64,
226        z0: f64,
227        kr0: f64,
228        kz0: f64,
229        n_steps: usize,
230        _dt_hint: f64,
231    ) -> FusionResult<RayTraceResult> {
232        self.validate_runtime_parameters()?;
233        if !r0.is_finite() || r0 <= 0.0 || !z0.is_finite() || !kr0.is_finite() || !kz0.is_finite() {
234            return Err(FusionError::ConfigError(
235                "RF ray initial conditions must be finite with r0 > 0".to_string(),
236            ));
237        }
238        if n_steps == 0 {
239            return Err(FusionError::ConfigError(
240                "RF ray tracing requires n_steps >= 1".to_string(),
241            ));
242        }
243
244        let r_res = self.resonance_radius()?;
245        let mut state: RayState = [r0, z0, kr0, kz0];
246        let mut trajectory = Vec::with_capacity(n_steps);
247        let mut resonance_point = None;
248
249        // Maximum spatial displacement per step [m]
250        const MAX_DISP: f64 = 0.05;
251
252        for _ in 0..n_steps {
253            if state.iter().any(|v| !v.is_finite()) {
254                return Err(FusionError::ConfigError(
255                    "RF ray state became non-finite".to_string(),
256                ));
257            }
258            trajectory.push((state[0], state[1]));
259
260            // Check resonance crossing
261            if resonance_point.is_none() && (state[0] - r_res).abs() < RESONANCE_TOL {
262                resonance_point = Some((state[0], state[1]));
263            }
264
265            // Adaptive dt: estimate RHS magnitude
266            let rhs = self.ray_rhs(&state);
267            if rhs.iter().any(|v| !v.is_finite()) {
268                return Err(FusionError::ConfigError(
269                    "RF ray RHS became non-finite".to_string(),
270                ));
271            }
272            let v_mag = (rhs[0] * rhs[0] + rhs[1] * rhs[1]).sqrt();
273            if !v_mag.is_finite() {
274                return Err(FusionError::ConfigError(
275                    "RF ray velocity magnitude became non-finite".to_string(),
276                ));
277            }
278            let dt = if v_mag > 1e-10 {
279                MAX_DISP / v_mag
280            } else {
281                1e-3
282            };
283            if !dt.is_finite() || dt <= 0.0 {
284                return Err(FusionError::ConfigError(
285                    "RF adaptive timestep became invalid".to_string(),
286                ));
287            }
288
289            state = self.rk4_step(&state, dt);
290
291            // Boundary check (stay within reasonable domain)
292            if state[0] < 0.5 || state[0] > 15.0 || state[1].abs() > 10.0 {
293                break;
294            }
295        }
296
297        Ok(RayTraceResult {
298            trajectory,
299            resonance_point,
300            r_resonance: r_res,
301        })
302    }
303
304    /// Trace all rays from the antenna array. Python lines 126-157.
305    pub fn trace_rays(&self) -> FusionResult<Vec<RayTraceResult>> {
306        self.validate_runtime_parameters()?;
307        let n_steps = 500;
308        let dt = 0.0; // adaptive
309
310        let mut results = Vec::with_capacity(N_RAYS);
311        for i in 0..N_RAYS {
312            let z_antenna = -1.0 + 2.0 * (i as f64) / ((N_RAYS - 1) as f64);
313            let result =
314                self.trace_single_ray(R_ANTENNA, z_antenna, -K0_INITIAL, 0.0, n_steps, dt)?;
315            results.push(result);
316        }
317        Ok(results)
318    }
319}
320
321impl Default for RFHeatingSystem {
322    fn default() -> Self {
323        Self::new()
324    }
325}
326
327#[cfg(test)]
328mod tests {
329    use super::*;
330
331    #[test]
332    fn test_resonance_radius() {
333        let rf = RFHeatingSystem::new();
334        let r_res = rf
335            .resonance_radius()
336            .expect("valid RF resonance-radius parameters");
337        // For ITER: B0=5.3T, R0=6.2m, f=50MHz
338        // B_res = 2π·50e6·3.34e-27/1.602e-19 ≈ 6.54 T
339        // R_res = 5.3·6.2/6.54 ≈ 5.02 m
340        assert!(r_res > 3.0 && r_res < 8.0, "Resonance radius: {r_res}");
341    }
342
343    #[test]
344    fn test_trace_single_ray() {
345        let rf = RFHeatingSystem::new();
346        let result = rf
347            .trace_single_ray(R_ANTENNA, 0.0, -K0_INITIAL, 0.0, 500, 0.0)
348            .expect("valid RF single-ray trace parameters");
349
350        // Ray should propagate inward
351        assert!(result.trajectory.len() > 1, "Should have trajectory points");
352
353        // First point at antenna
354        let (r0, _z0) = result.trajectory[0];
355        assert!((r0 - R_ANTENNA).abs() < 1e-10, "Should start at antenna");
356    }
357
358    #[test]
359    fn test_trace_rays_array() {
360        let rf = RFHeatingSystem::new();
361        let results = rf
362            .trace_rays()
363            .expect("valid RF multi-ray trace parameters");
364        assert_eq!(results.len(), N_RAYS);
365        for result in &results {
366            assert!(!result.trajectory.is_empty());
367        }
368    }
369
370    #[test]
371    fn test_alfven_speed_positive() {
372        let rf = RFHeatingSystem::new();
373        let v_a = rf.alfven_speed(6.2, 0.0);
374        assert!(
375            v_a > 0.0 && v_a.is_finite(),
376            "v_A should be positive: {v_a}"
377        );
378    }
379
380    #[test]
381    fn test_rf_heating_rejects_invalid_runtime_inputs() {
382        let mut rf = RFHeatingSystem::new();
383        rf.b0 = f64::NAN;
384        assert!(rf.resonance_radius().is_err());
385        assert!(rf.trace_rays().is_err());
386
387        let rf = RFHeatingSystem::new();
388        assert!(rf
389            .trace_single_ray(0.0, 0.0, -K0_INITIAL, 0.0, 500, 0.0)
390            .is_err());
391        assert!(rf
392            .trace_single_ray(R_ANTENNA, 0.0, -K0_INITIAL, 0.0, 0, 0.0)
393            .is_err());
394    }
395}
fusion_core/rf_heating.rs

fusion_core/
rf_heating.rs
