fusion_core/
rf_heating.rs

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

fusion_core/
rf_heating.rs
