fusion_core/
vacuum.rs

1// ─────────────────────────────────────────────────────────────────────
2// SCPN Fusion Core — Vacuum
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//! Vacuum magnetic field from coils via toroidal Green's function.
9//!
10//! Port of fusion_kernel.py `calculate_vacuum_field()` (lines 59-93).
11//! Uses complete elliptic integrals K(m), E(m) for the exact toroidal
12//! single-coil flux function (Smythe/Jackson form).
13
14use fusion_math::elliptic::{ellipe, ellipk};
15use fusion_types::config::CoilConfig;
16use fusion_types::error::{FusionError, FusionResult};
17use fusion_types::state::Grid2D;
18use ndarray::Array2;
19
20/// Calculate vacuum magnetic flux from coils using the toroidal Green's function.
21///
22/// For each coil at `(Rc, Zc)` with current `I`, the flux contribution is:
23///
24///   Ψ = (μ₀ I / 2π) · √((R + Rc)² + (Z - Zc)²) · ((2 - k²)K(k²) - 2E(k²)) / k²
25///
26/// where k² = 4 R Rc / ((R + Rc)² + (Z - Zc)²), clipped to [1e-9, 0.999999].
27///
28/// Returns Ψ_vacuum `[nz, nr]`.
29pub fn calculate_vacuum_field(
30    grid: &Grid2D,
31    coils: &[CoilConfig],
32    mu0: f64,
33) -> FusionResult<Array2<f64>> {
34    if !mu0.is_finite() || mu0 <= 0.0 {
35        return Err(FusionError::ConfigError(format!(
36            "vacuum permeability mu0 must be finite and > 0, got {mu0}"
37        )));
38    }
39    if grid.nz == 0 || grid.nr == 0 {
40        return Err(FusionError::ConfigError(
41            "vacuum field grid must have nz,nr >= 1".to_string(),
42        ));
43    }
44    if grid.rr.nrows() != grid.nz || grid.rr.ncols() != grid.nr {
45        return Err(FusionError::ConfigError(format!(
46            "grid.rr shape mismatch: expected ({}, {}), got ({}, {})",
47            grid.nz,
48            grid.nr,
49            grid.rr.nrows(),
50            grid.rr.ncols()
51        )));
52    }
53    if grid.zz.nrows() != grid.nz || grid.zz.ncols() != grid.nr {
54        return Err(FusionError::ConfigError(format!(
55            "grid.zz shape mismatch: expected ({}, {}), got ({}, {})",
56            grid.nz,
57            grid.nr,
58            grid.zz.nrows(),
59            grid.zz.ncols()
60        )));
61    }
62    if grid.rr.iter().any(|v| !v.is_finite()) || grid.zz.iter().any(|v| !v.is_finite()) {
63        return Err(FusionError::ConfigError(
64            "vacuum field grid coordinates must be finite".to_string(),
65        ));
66    }
67
68    let nz = grid.nz;
69    let nr = grid.nr;
70    let mut psi_vac: Array2<f64> = Array2::zeros((nz, nr));
71
72    for coil in coils {
73        let rc = coil.r;
74        let zc = coil.z;
75        let current = coil.current;
76        if !rc.is_finite() || rc <= 0.0 {
77            return Err(FusionError::ConfigError(format!(
78                "coil radius must be finite and > 0, got {}",
79                coil.r
80            )));
81        }
82        if !zc.is_finite() || !current.is_finite() {
83            return Err(FusionError::ConfigError(
84                "coil z/current must be finite".to_string(),
85            ));
86        }
87
88        let prefactor = (mu0 * current) / (2.0 * std::f64::consts::PI);
89        if !prefactor.is_finite() {
90            return Err(FusionError::ConfigError(
91                "coil prefactor became non-finite".to_string(),
92            ));
93        }
94
95        for iz in 0..nz {
96            for ir in 0..nr {
97                let r = grid.rr[[iz, ir]];
98                let z = grid.zz[[iz, ir]];
99
100                let dz = z - zc;
101                let r_plus_rc = r + rc;
102                let r_plus_rc_sq = r_plus_rc * r_plus_rc;
103                let denom = r_plus_rc_sq + dz * dz;
104                if !denom.is_finite() || denom <= 0.0 {
105                    return Err(FusionError::ConfigError(
106                        "vacuum field denominator became non-finite or non-positive".to_string(),
107                    ));
108                }
109
110                // k² parameter
111                let k2_raw = (4.0 * r * rc) / denom;
112                if !k2_raw.is_finite() {
113                    return Err(FusionError::ConfigError(
114                        "vacuum field k^2 became non-finite".to_string(),
115                    ));
116                }
117                let mut k2 = k2_raw;
118                // CRITICAL: clip to avoid singularity at ellipk(1.0) and division by k2=0
119                k2 = k2.clamp(1e-9, 0.999999);
120
121                let k_val = ellipk(k2);
122                let e_val = ellipe(k2);
123                if !k_val.is_finite() || !e_val.is_finite() {
124                    return Err(FusionError::ConfigError(
125                        "vacuum field elliptic integrals became non-finite".to_string(),
126                    ));
127                }
128
129                let sqrt_term = (r * rc).sqrt();
130                let k = k2.sqrt();
131                let term = ((2.0 - k2) * k_val - 2.0 * e_val) / k;
132                let coil_flux = prefactor * sqrt_term * term;
133                if !coil_flux.is_finite() {
134                    return Err(FusionError::ConfigError(
135                        "vacuum field coil flux contribution became non-finite".to_string(),
136                    ));
137                }
138
139                psi_vac[[iz, ir]] += coil_flux;
140            }
141        }
142    }
143
144    if psi_vac.iter().any(|v| !v.is_finite()) {
145        return Err(FusionError::ConfigError(
146            "vacuum field result contains non-finite values".to_string(),
147        ));
148    }
149
150    Ok(psi_vac)
151}
152
153#[cfg(test)]
154mod tests {
155    use super::*;
156    use fusion_types::config::ReactorConfig;
157    use std::path::PathBuf;
158
159    fn project_root() -> PathBuf {
160        PathBuf::from(env!("CARGO_MANIFEST_DIR"))
161            .join("..")
162            .join("..")
163            .join("..")
164    }
165
166    fn config_path(relative: &str) -> String {
167        project_root().join(relative).to_string_lossy().to_string()
168    }
169
170    #[test]
171    fn test_vacuum_field_not_nan() {
172        let cfg = ReactorConfig::from_file(&config_path("iter_config.json")).unwrap();
173        let grid = cfg.create_grid();
174        let psi_vac = calculate_vacuum_field(&grid, &cfg.coils, cfg.physics.vacuum_permeability)
175            .expect("valid vacuum-field inputs");
176
177        // No NaN values
178        assert!(
179            !psi_vac.iter().any(|v| v.is_nan()),
180            "Vacuum field contains NaN"
181        );
182
183        // Should have non-zero values (coils produce flux)
184        let max_val = psi_vac.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
185        assert!(max_val.abs() > 1e-10, "Vacuum field is all zeros");
186    }
187
188    #[test]
189    fn test_vacuum_field_symmetry() {
190        // Single coil at Z=0 should produce field symmetric about Z=0
191        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
192        let coil = CoilConfig {
193            name: "test".to_string(),
194            r: 5.0,
195            z: 0.0,
196            current: 1.0,
197        };
198        let psi = calculate_vacuum_field(&grid, &[coil], 1.0).expect("valid vacuum-field inputs");
199
200        // Check symmetry: psi[iz, ir] ≈ psi[nz-1-iz, ir]
201        for ir in 0..33 {
202            for iz in 0..16 {
203                let diff = (psi[[iz, ir]] - psi[[32 - iz, ir]]).abs();
204                assert!(
205                    diff < 1e-10,
206                    "Symmetry broken at iz={iz}, ir={ir}: {} vs {}",
207                    psi[[iz, ir]],
208                    psi[[32 - iz, ir]]
209                );
210            }
211        }
212    }
213
214    #[test]
215    fn test_vacuum_field_single_coil_shape() {
216        // Flux from a single coil should peak near the coil position
217        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
218        let coil = CoilConfig {
219            name: "test".to_string(),
220            r: 5.0,
221            z: 0.0,
222            current: 1.0,
223        };
224        let psi = calculate_vacuum_field(&grid, &[coil], 1.0).expect("valid vacuum-field inputs");
225
226        // Find maximum
227        let mut max_val = f64::NEG_INFINITY;
228        let mut max_iz = 0;
229        let mut max_ir = 0;
230        for iz in 0..33 {
231            for ir in 0..33 {
232                if psi[[iz, ir]] > max_val {
233                    max_val = psi[[iz, ir]];
234                    max_iz = iz;
235                    max_ir = ir;
236                }
237            }
238        }
239
240        // Maximum should be near Z=0 (middle row)
241        assert!(
242            (max_iz as f64 - 16.0).abs() < 3.0,
243            "Peak not near Z=0: iz={max_iz}"
244        );
245        // Maximum should be near R=5.0 (middle of grid)
246        assert!(
247            (max_ir as f64 - 16.0).abs() < 3.0,
248            "Peak not near R=5: ir={max_ir}"
249        );
250    }
251
252    #[test]
253    fn test_vacuum_field_rejects_invalid_runtime_inputs() {
254        let grid = Grid2D::new(17, 17, 1.0, 9.0, -5.0, 5.0);
255        let bad_coil = CoilConfig {
256            name: "bad".to_string(),
257            r: 0.0,
258            z: 0.0,
259            current: 1.0,
260        };
261        assert!(calculate_vacuum_field(&grid, &[bad_coil], 1.0).is_err());
262        assert!(calculate_vacuum_field(&grid, &[], f64::NAN).is_err());
263    }
264}
fusion_core/vacuum.rs

fusion_core/
vacuum.rs
