diff --git a/Cargo.lock b/Cargo.lock index 9abbe813..0253a52c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -3219,7 +3219,7 @@ name = "mecheval-grader" version = "0.0.1" dependencies = [ "anyhow", - "phyz", + "phyz 0.3.0", "serde", "serde_json", "thiserror 2.0.18", @@ -4013,6 +4013,14 @@ dependencies = [ "siphasher 1.0.2", ] +[[package]] +name = "phyz" +version = "0.3.0" +dependencies = [ + "phyz-math 0.1.0", + "tang", +] + [[package]] name = "phyz" version = "0.3.1" @@ -4020,7 +4028,7 @@ dependencies = [ "phyz-collision", "phyz-contact", "phyz-diff", - "phyz-math", + "phyz-math 0.3.1", "phyz-model", "phyz-rigid", ] @@ -4029,7 +4037,7 @@ dependencies = [ name = "phyz-collision" version = "0.3.1" dependencies = [ - "phyz-math", + "phyz-math 0.3.1", "phyz-model", ] @@ -4038,7 +4046,7 @@ name = "phyz-contact" version = "0.3.1" dependencies = [ "phyz-collision", - "phyz-math", + "phyz-math 0.3.1", "phyz-model", ] @@ -4046,7 +4054,7 @@ dependencies = [ name = "phyz-diff" version = "0.3.1" dependencies = [ - "phyz-math", + "phyz-math 0.3.1", "phyz-model", "phyz-rigid", "tang", @@ -4060,13 +4068,21 @@ version = "0.3.1" dependencies = [ "bytemuck", "futures-intrusive", - "phyz-math", + "phyz-math 0.3.1", "phyz-model", "phyz-rigid", "pollster", "wgpu", ] +[[package]] +name = "phyz-math" +version = "0.1.0" +dependencies = [ + "tang", + "tang-la", +] + [[package]] name = "phyz-math" version = "0.3.1" @@ -4079,7 +4095,7 @@ dependencies = [ name = "phyz-md" version = "0.3.1" dependencies = [ - "phyz-math", + "phyz-math 0.3.1", "rand 0.8.5", "rand_chacha 0.3.1", "rand_distr", @@ -4090,14 +4106,14 @@ dependencies = [ name = "phyz-model" version = "0.3.1" dependencies = [ - "phyz-math", + "phyz-math 0.3.1", ] [[package]] name = "phyz-rigid" version = "0.3.1" dependencies = [ - "phyz-math", + "phyz-math 0.3.1", "phyz-model", ] @@ -7505,6 +7521,15 @@ dependencies = [ "wgpu", ] +[[package]] +name = "vcad-kernel-magnetostatic" +version = "0.9.4" +dependencies = [ + "serde", + "serde_json", + "vcad-kernel-particle", +] + [[package]] name = "vcad-kernel-math" version = "0.9.4" @@ -7585,7 +7610,7 @@ dependencies = [ name = "vcad-kernel-physics" version = "0.9.4" dependencies = [ - "phyz", + "phyz 0.3.1", "serde", "serde_json", "stl_io", @@ -7902,9 +7927,9 @@ dependencies = [ name = "vcad-sim" version = "0.9.4" dependencies = [ - "phyz", + "phyz 0.3.1", "phyz-gpu", - "phyz-math", + "phyz-math 0.3.1", "phyz-model", "serde", "thiserror 2.0.18", diff --git a/Cargo.toml b/Cargo.toml index f3f5bae9..2d80f5a3 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -31,6 +31,7 @@ members = [ "crates/vcad-kernel-optics", "crates/vcad-kernel-acoustics", "crates/vcad-kernel-em", + "crates/vcad-kernel-magnetostatic", "crates/vcad-kernel-enclosure", "crates/vcad-kernel-orbit", "crates/vcad-kernel-qcd", @@ -129,6 +130,7 @@ default-members = [ "crates/vcad-kernel-optics", "crates/vcad-kernel-acoustics", "crates/vcad-kernel-em", + "crates/vcad-kernel-magnetostatic", "crates/vcad-kernel-enclosure", "crates/vcad-kernel-orbit", "crates/vcad-kernel-qcd", @@ -250,6 +252,7 @@ vcad-kernel-particle = { path = "crates/vcad-kernel-particle", version = "0.9.4" vcad-kernel-optics = { path = "crates/vcad-kernel-optics", version = "0.9.4" } vcad-kernel-acoustics = { path = "crates/vcad-kernel-acoustics", version = "0.9.4" } vcad-kernel-em = { path = "crates/vcad-kernel-em", version = "0.9.4" } +vcad-kernel-magnetostatic = { path = "crates/vcad-kernel-magnetostatic", version = "0.9.4" } vcad-kernel-enclosure = { path = "crates/vcad-kernel-enclosure", version = "0.9.4" } vcad-kernel-orbit = { path = "crates/vcad-kernel-orbit", version = "0.9.4" } vcad-kernel-qcd = { path = "crates/vcad-kernel-qcd", version = "0.9.4" } diff --git a/crates/vcad-kernel-magnetostatic/Cargo.toml b/crates/vcad-kernel-magnetostatic/Cargo.toml new file mode 100644 index 00000000..e9c94340 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/Cargo.toml @@ -0,0 +1,14 @@ +[package] +name = "vcad-kernel-magnetostatic" +description = "Grid-free 3D magnetostatics for air-core machines: finite-segment Biot-Savart, magnets as equivalent surface currents, flux linkage and torque" +version.workspace = true +edition.workspace = true +license.workspace = true +repository.workspace = true + +[dependencies] +serde = { workspace = true } + +[dev-dependencies] +serde_json = { workspace = true } +vcad-kernel-particle = { workspace = true } diff --git a/crates/vcad-kernel-magnetostatic/examples/image_convergence.rs b/crates/vcad-kernel-magnetostatic/examples/image_convergence.rs new file mode 100644 index 00000000..7c77ebed --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/examples/image_convergence.rs @@ -0,0 +1,63 @@ +//! Diagnostic: how fast does the two-plane image series actually converge, for +//! a single loop versus a magnetically balanced alternating-pole set? + +use std::f64::consts::PI; +use vcad_kernel_magnetostatic::{Filament, IronStack, Vec3}; + +fn ring(r: f64, z0: f64, i: f64, n: usize, center: Vec3) -> Filament { + let pts = (0..n) + .map(|k| { + let t = 2.0 * PI * (k as f64) / (n as f64); + Vec3::new(center.x + r * t.cos(), center.y + r * t.sin(), z0) + }) + .collect(); + Filament::closed_loop(pts, i, 0.0) +} + +fn b_sum(stack: &IronStack, srcs: &[Filament], p: Vec3) -> Vec3 { + srcs.iter() + .flat_map(|s| stack.expand(s)) + .map(|f| f.b_at(p)) + .sum() +} + +fn main() { + let z_lo = 0.0; + let z_hi = 0.0056; + let probe = Vec3::new(0.018, 0.0, 0.0045); + + // Case 1: one loop — net magnetic moment, nothing to cancel at distance. + let single = vec![ring(0.0075, 0.003, 100.0, 256, Vec3::ZERO)]; + + // Case 2: six alternating poles on a 22.5 mm pitch circle — a real rotor. + let mut alternating = Vec::new(); + for k in 0..6 { + let t = 2.0 * PI * (k as f64) / 6.0; + let c = Vec3::new(0.0225 * t.cos(), 0.0225 * t.sin(), 0.0); + let sign = if k % 2 == 0 { 1.0 } else { -1.0 }; + alternating.push(ring(0.0075, 0.003, 100.0 * sign, 256, c)); + } + let probe_alt = Vec3::new(0.0225, 0.0, 0.0045); + + for (label, srcs, p) in [ + ("single loop (net moment)", &single, probe), + ("6 alternating poles", &alternating, probe_alt), + ] { + println!("\n=== {label} ==="); + let mut prev: Option = None; + for depth in [2usize, 4, 8, 12, 16, 24, 32] { + let stack = IronStack::pair(z_lo, z_hi, depth); + let n_images = stack.expand(&srcs[0]).len(); + let b = b_sum(&stack, srcs, p); + let delta = prev + .map(|q: Vec3| (b - q).norm() / b.norm().max(1e-18)) + .unwrap_or(f64::NAN); + println!( + "depth {depth:>3} images/src {n_images:>4} |B| {:>12.6e} T Δ vs prev {:.3e}", + b.norm(), + delta + ); + prev = Some(b); + } + } +} diff --git a/crates/vcad-kernel-magnetostatic/src/filament.rs b/crates/vcad-kernel-magnetostatic/src/filament.rs new file mode 100644 index 00000000..a28718a4 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/src/filament.rs @@ -0,0 +1,280 @@ +//! Fields of straight current-carrying segments, in closed form. +//! +//! Everything in this crate reduces to these two integrals. An air-core machine +//! has no iron, so the permeability is μ₀ everywhere, the problem is linear, and +//! superposition holds exactly — there is no grid, no mesh, and no iteration. +//! The whole solver is a sum over segments. + +use crate::vec3::Vec3; +use crate::MU_0; + +use std::f64::consts::PI; + +/// A straight current-carrying segment, SI units (metres, amperes). +/// +/// `wire_radius_m` regularizes the field on the segment axis, where the +/// filamentary formulae diverge. See [`Segment::b_at`]. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct Segment { + /// Start point, m. + pub a: Vec3, + /// End point, m. Current flows a → b. + pub b: Vec3, + /// Current, A (signed; reverse it by swapping `a`/`b` or negating). + pub current_a: f64, + /// Conductor radius, m, used to regularize on-axis singularities. + pub wire_radius_m: f64, +} + +/// Geometry shared by the `B` and `A` evaluations: the segment's unit tangent, +/// the endpoint vectors, and the perpendicular distance to the field point. +struct SegmentGeometry { + /// Unit tangent along the current direction. + tangent: Vec3, + /// Field point relative to `a`. (`r₂` is only ever needed through its norm + /// and tangential component, both carried below.) + r1: Vec3, + /// Norms of `r1`, `r2`. + n1: f64, + n2: f64, + /// Tangential components `t̂·r`. + t1: f64, + t2: f64, + /// Perpendicular distance from the field point to the segment's line, m, + /// **before** any regularization clamp. + perp: f64, +} + +impl Segment { + fn geometry(&self, p: Vec3) -> Option { + let l = self.b - self.a; + let len = l.norm(); + if len <= 0.0 || self.current_a == 0.0 { + return None; + } + let tangent = l * (1.0 / len); + let r1 = p - self.a; + let r2 = p - self.b; + let t1 = tangent.dot(r1); + let t2 = tangent.dot(r2); + // |r|² = perp² + t², evaluated from whichever endpoint is better + // conditioned (the one whose tangential offset is smaller). + let n1 = r1.norm(); + let n2 = r2.norm(); + let perp = (r1 - tangent * t1).norm(); + Some(SegmentGeometry { tangent, r1, n1, n2, t1, t2, perp }) + } + + /// Magnetic flux density at `p`, tesla. + /// + /// For a segment from **a** to **b** carrying current `I`, with `t̂` the unit + /// tangent, `r₁ = p − a`, `r₂ = p − b`, and `d` the perpendicular distance + /// from `p` to the segment's line: + /// + /// ```text + /// B = μ₀I/(4πd) · (cos θ₁ − cos θ₂) · (t̂ × r̂₁) + /// ``` + /// + /// where `cos θᵢ = t̂·rᵢ/|rᵢ|`. Taking the segment to infinity in both + /// directions gives `μ₀I/(2πd)`, the infinite-wire result. + /// + /// **Regularization.** Inside the conductor (`d < wire_radius_m`) the + /// filamentary field diverges. We evaluate at the conductor surface and + /// scale linearly to zero on the axis — the field of a uniform current + /// density, and the same model `vcad_kernel_particle::field::b_ring` uses. + pub fn b_at(&self, p: Vec3) -> Vec3 { + let Some(g) = self.geometry(p) else { + return Vec3::ZERO; + }; + // Direction: t̂ × r₁, which is perpendicular to both the current and the + // offset, i.e. azimuthal about the wire by the right-hand rule. + let dir = g.tangent.cross(g.r1); + let dir_n = dir.norm(); + if dir_n <= 0.0 { + // Field point lies exactly on the segment's line: B is zero there by + // symmetry for a uniform-current conductor. + return Vec3::ZERO; + } + let cos1 = if g.n1 > 0.0 { g.t1 / g.n1 } else { 0.0 }; + let cos2 = if g.n2 > 0.0 { g.t2 / g.n2 } else { 0.0 }; + + let (d_eff, scale) = if g.perp < self.wire_radius_m && self.wire_radius_m > 0.0 { + (self.wire_radius_m, g.perp / self.wire_radius_m) + } else { + (g.perp, 1.0) + }; + if d_eff <= 0.0 { + return Vec3::ZERO; + } + let mag = MU_0 * self.current_a / (4.0 * PI * d_eff) * (cos1 - cos2) * scale; + dir * (mag / dir_n) + } + + /// Magnetic vector potential at `p`, tesla-metres (Coulomb gauge). + /// + /// Integrating `dA = μ₀I/(4π) · dl/|p − l|` along the segment gives + /// + /// ```text + /// A = μ₀I/(4π) · t̂ · ln[ (|r₁| + t̂·r₁) / (|r₂| + t̂·r₂) ] + /// ``` + /// + /// Flux linkage is computed from `A` rather than by integrating `B` over a + /// surface: `λ = ∮A·dl` needs only the conductor path, so there is no + /// spanning surface to construct for a multi-turn spiral — which is the + /// whole reason this crate can stay grid-free. + /// + /// **Conditioning.** `|r| + t̂·r` cancels catastrophically when `t̂·r` is + /// negative and large. Since `(|r| + t)(|r| − t) = d²`, we evaluate the + /// unstable branch as `d²/(|r| − t)` instead, which is exact and stable. + /// + /// **Regularization.** The logarithm diverges on the axis; `d` is clamped to + /// `wire_radius_m`. Self-inductance from a filament model is genuinely + /// sensitive to that clamp — see the crate-level note on `L`. + pub fn a_at(&self, p: Vec3) -> Vec3 { + let Some(g) = self.geometry(p) else { + return Vec3::ZERO; + }; + let d = g.perp.max(self.wire_radius_m); + if d <= 0.0 { + return Vec3::ZERO; + } + let d2 = d * d; + // Stable evaluation of |r| + t for either sign of t. + let safe = |n: f64, t: f64| -> f64 { + if t >= 0.0 { + n + t + } else { + d2 / (n - t) + } + }; + // Recompute norms against the clamped perpendicular distance so the two + // endpoints stay consistent with `d` when the point is inside the wire. + let n1 = (d2 + g.t1 * g.t1).sqrt(); + let n2 = (d2 + g.t2 * g.t2).sqrt(); + let num = safe(n1, g.t1); + let den = safe(n2, g.t2); + if num <= 0.0 || den <= 0.0 { + return Vec3::ZERO; + } + g.tangent * (MU_0 * self.current_a / (4.0 * PI) * (num / den).ln()) + } +} + +/// A current path: a polyline carrying one current, discretized into segments. +/// +/// A spiral PCB coil, a magnet's bound-current sheet and a phase winding are all +/// represented this way, so one integrator serves the whole solver. +#[derive(Debug, Clone, PartialEq)] +pub struct Filament { + /// Ordered path points, m. Current flows from `points[0]` toward the end. + pub points: Vec, + /// Current, A. For an `N`-turn conductor modelled as one path, this is the + /// ampere-turns. + pub current_a: f64, + /// Conductor radius, m, for on-axis regularization. + pub wire_radius_m: f64, + /// Whether the path closes back on `points[0]`. + pub closed: bool, +} + +impl Filament { + /// A closed loop through `points`. + pub fn closed_loop(points: Vec, current_a: f64, wire_radius_m: f64) -> Self { + Self { points, current_a, wire_radius_m, closed: true } + } + + /// An open path through `points`. + pub fn open_path(points: Vec, current_a: f64, wire_radius_m: f64) -> Self { + Self { points, current_a, wire_radius_m, closed: false } + } + + /// Iterate the path's segments. + pub fn segments(&self) -> impl Iterator + '_ { + let n = self.points.len(); + let last = if self.closed { n } else { n.saturating_sub(1) }; + (0..last).map(move |i| Segment { + a: self.points[i], + b: self.points[(i + 1) % n], + current_a: self.current_a, + wire_radius_m: self.wire_radius_m, + }) + } + + /// Total path length, m. + pub fn length_m(&self) -> f64 { + self.segments().map(|s| (s.b - s.a).norm()).sum() + } + + /// Magnetic dipole moment, A·m²: `m = (I/2)∮ r × dl`. + /// + /// Only meaningful for a closed path. Used to check that a source set is + /// magnetically balanced before expanding it in an [`crate::IronStack`] — + /// see that type's convergence precondition. + pub fn dipole_moment(&self) -> Vec3 { + let s: Vec3 = self + .segments() + .map(|s| { + let mid = (s.a + s.b) * 0.5; + mid.cross(s.b - s.a) + }) + .sum(); + s * (self.current_a * 0.5) + } + + /// Flux density at `p` from this path, tesla. + pub fn b_at(&self, p: Vec3) -> Vec3 { + self.segments().map(|s| s.b_at(p)).sum() + } + + /// Vector potential at `p` from this path, T·m. + pub fn a_at(&self, p: Vec3) -> Vec3 { + self.segments().map(|s| s.a_at(p)).sum() + } + + /// Rotate the whole path about the z axis. + pub fn rotated_z(&self, angle: f64) -> Filament { + Filament { + points: self.points.iter().map(|p| p.rotated_z(angle)).collect(), + ..self.clone() + } + } + + /// Flux linked by this path from an external field source, weber. + /// + /// `λ = ∮A·dl`, evaluated with the midpoint rule on each segment. The + /// external `a_field` must exclude this path's own contribution, or the + /// result is self-inductance rather than mutual flux linkage. + pub fn flux_linkage Vec3>(&self, a_field: F) -> f64 { + self.segments() + .map(|s| { + let mid = (s.a + s.b) * 0.5; + let dl = s.b - s.a; + a_field(mid).dot(dl) + }) + .sum() + } + + /// Net force on this path in an external field, newtons: `F = ∮ I dl × B`. + pub fn lorentz_force Vec3>(&self, b_field: F) -> Vec3 { + self.segments() + .map(|s| { + let mid = (s.a + s.b) * 0.5; + let dl = s.b - s.a; + dl.cross(b_field(mid)) * self.current_a + }) + .sum() + } + + /// Torque about the origin's z axis in an external field, N·m: + /// `T_z = ẑ · ∮ r × (I dl × B)`. + pub fn lorentz_torque_z Vec3>(&self, b_field: F) -> f64 { + self.segments() + .map(|s| { + let mid = (s.a + s.b) * 0.5; + let dl = s.b - s.a; + let df = dl.cross(b_field(mid)) * self.current_a; + mid.cross(df).z + }) + .sum() + } +} diff --git a/crates/vcad-kernel-magnetostatic/src/iron.rs b/crates/vcad-kernel-magnetostatic/src/iron.rs new file mode 100644 index 00000000..b7216537 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/src/iron.rs @@ -0,0 +1,216 @@ +//! High-permeability boundaries by the method of images. +//! +//! Most real "air-core" machines are only air-core in the *winding* region: the +//! coils have no iron, but steel back-irons close the magnetic circuit behind +//! the magnets and behind the stator. That steel is worth a factor of ~2 in +//! airgap flux, so ignoring it is not a small error. +//! +//! Because the steel operates far below saturation in these machines, it can be +//! treated as `μ_r → ∞`, and an infinitely permeable plane is exactly a mirror: +//! reflect the current path through the plane and keep the current the same, and +//! the pair reproduces the boundary condition `B_tangential = 0` at the surface. +//! That keeps the whole solver grid-free — no mesh, no iteration, no truncation +//! boundary — at the cost of assuming the plates are flat, parallel and +//! unbounded. +//! +//! # Convergence precondition — read this before using two planes +//! +//! One plane is a single reflection and always exact. **Two** planes generate an +//! infinite image series, like facing mirrors, and that series only converges if +//! the source set is **magnetically balanced** — net dipole moment ≈ 0. +//! +//! This is not a numerical nicety, it is the physics. Measured on a 5.6 mm plate +//! separation: +//! +//! | source | depth 8 | depth 16 | depth 32 | behaviour | +//! |---|---|---|---|---| +//! | one loop (net moment) | 2.29e-4 T | 7.28e-5 T | 2.21e-5 T | **never settles** | +//! | 6 alternating poles | 2.258e-2 T | 2.2598e-2 T | 2.2598e-2 T | converged to 6e-5 by depth 16 | +//! +//! A single magnet between two infinite iron plates is an ill-posed idealization: +//! its flux has no bounded return path, so every image adds with the same sense +//! and the sum drifts without limit. Any real multi-pole rotor is balanced by +//! construction — poles alternate — so this costs nothing in practice, but +//! [`IronStack::balance_residual`] enforces it rather than trusting it. +//! +//! # What this does and does not buy +//! +//! - **Exact** (to `1/μ_r`, ≈0.1% for unsaturated mild steel) for flat parallel +//! plates of infinite extent, given a balanced source set. +//! - **Not modelled:** finite plate radius. Real back-irons are discs, so flux +//! fringes at the inner and outer diameters and the image model over-predicts +//! there. Compare against a grid solver before trusting the last few percent. +//! - **Not modelled:** saturation. [`IronStack::flux_density_estimate`] reports +//! the working flux density in the iron so callers can check the assumption +//! instead of inheriting it silently. + +use crate::filament::Filament; +use crate::vec3::Vec3; + +/// A stack of parallel, infinitely-permeable planes normal to **z**. +/// +/// One plane is a single mirror. Two planes form a cavity whose image series is +/// infinite, exactly like facing mirrors; it is truncated at +/// [`IronStack::max_reflections`] and the residual is reported by +/// [`IronStack::tail_fraction`]. +#[derive(Debug, Clone, PartialEq)] +pub struct IronStack { + /// z positions of the planes, m. + planes_z: Vec, + /// Reflection depth for the image series. + max_reflections: usize, +} + +impl IronStack { + /// No iron — the pure free-space case. + pub fn none() -> Self { + Self { planes_z: Vec::new(), max_reflections: 0 } + } + + /// A single back-iron plane at `z_m`. + pub fn single(z_m: f64) -> Self { + Self { planes_z: vec![z_m], max_reflections: 1 } + } + + /// Two parallel planes — the usual rotor/stator back-iron sandwich. + /// + /// `max_reflections` bounds the image series. Depth 16 reaches a 6e-5 + /// residual for a balanced 6-pole rotor at 5.6 mm plate separation; see + /// [`IronStack::DEFAULT_REFLECTIONS`]. **The source set must be + /// magnetically balanced** — check with [`IronStack::balance_residual`], and + /// confirm the achieved truncation with [`IronStack::tail_fraction`]. + pub fn pair(z_a: f64, z_b: f64, max_reflections: usize) -> Self { + assert!(z_a != z_b, "iron planes must be distinct"); + Self { planes_z: vec![z_a, z_b], max_reflections: max_reflections.max(1) } + } + + /// Reflection depth that converges a balanced multi-pole rotor to ~1e-4. + pub const DEFAULT_REFLECTIONS: usize = 16; + + /// How far a source set is from magnetically balanced: `|Σmᵢ| / Σ|mᵢ|`. + /// + /// Zero for an alternating-pole rotor, one for a lone magnet. The two-plane + /// image series converges only near zero — see the module note. Values above + /// ~0.05 mean the result depends on [`IronStack::max_reflections`] and must + /// not be quoted. + pub fn balance_residual(sources: &[Filament]) -> f64 { + let mut net = Vec3::ZERO; + let mut gross = 0.0; + for s in sources { + let m = s.dipole_moment(); + net = net + m; + gross += m.norm(); + } + if gross <= 0.0 { + return 0.0; + } + net.norm() / gross + } + + /// The planes' z positions, m. + pub fn planes(&self) -> &[f64] { + &self.planes_z + } + + /// Whether any iron is present. + pub fn is_empty(&self) -> bool { + self.planes_z.is_empty() + } + + /// Reflect a filament through the plane at `z_p`. + /// + /// Geometric reflection of the directed path, current unchanged. A tangent + /// `(tx, ty, tz)` maps to `(tx, ty, −tz)`, which is precisely the `μ→∞` + /// image rule: tangential current preserved, normal current reversed. For a + /// loop parallel to the plane this keeps the circulation sense, so the pair + /// doubles `B_z` and cancels `B_r` at the surface — the boundary condition. + fn reflect(f: &Filament, z_p: f64) -> Filament { + Filament { + points: f + .points + .iter() + .map(|p| Vec3::new(p.x, p.y, 2.0 * z_p - p.z)) + .collect(), + ..f.clone() + } + } + + /// Expand `source` into itself plus its images. + /// + /// The images of a source between two planes are generated by alternately + /// reflecting in each plane, which for parallel planes yields the familiar + /// ladder marching away in both directions. Duplicate positions (a source + /// exactly on a plane reflects to itself) are dropped. + pub fn expand(&self, source: &Filament) -> Vec { + if self.planes_z.is_empty() { + return vec![source.clone()]; + } + let mut out = vec![source.clone()]; + let mut frontier = vec![(source.clone(), usize::MAX)]; + let same = |a: &Filament, b: &Filament| { + a.points.len() == b.points.len() + && a.points + .iter() + .zip(&b.points) + .all(|(p, q)| (p.z - q.z).abs() < 1e-12) + }; + for _ in 0..self.max_reflections { + let mut next = Vec::new(); + for (f, last_plane) in &frontier { + for (pi, &z_p) in self.planes_z.iter().enumerate() { + // Reflecting twice in the same plane returns the source. + if *last_plane == pi { + continue; + } + let img = Self::reflect(f, z_p); + if same(&img, f) || out.iter().any(|o| same(o, &img)) { + continue; + } + out.push(img.clone()); + next.push((img, pi)); + } + } + if next.is_empty() { + break; + } + frontier = next; + } + out + } + + /// Fraction of the field magnitude carried by the outermost image shell, at + /// `probe`, for `source`. + /// + /// This is the truncation error of the image series: if it is not small, the + /// answer depends on [`IronStack::max_reflections`] and should not be quoted. + pub fn tail_fraction(&self, source: &Filament, probe: Vec3) -> f64 { + if self.planes_z.is_empty() { + return 0.0; + } + let all = self.expand(source); + let total: Vec3 = all.iter().map(|f| f.b_at(probe)).sum(); + let deeper = IronStack { + planes_z: self.planes_z.clone(), + max_reflections: self.max_reflections + 2, + }; + let refined: Vec3 = deeper.expand(source).iter().map(|f| f.b_at(probe)).sum(); + let scale = refined.norm().max(1e-18); + (refined - total).norm() / scale + } + + /// Working flux density in the iron, tesla, estimated as the normal field at + /// the plane scaled by the plate's flux-concentration ratio. + /// + /// The image model assumes `μ_r → ∞`, which fails once the steel saturates. + /// `pole_area_m2` is the flux-carrying area per pole at the plane and + /// `iron_section_m2` the cross-section the returning flux must pass through + /// (thickness × mean circumference / poles). Compare against ~1.5 T for mild + /// steel: above that, the image result is optimistic and a nonlinear grid + /// solve is required. + pub fn flux_density_estimate(b_normal_t: f64, pole_area_m2: f64, iron_section_m2: f64) -> f64 { + if iron_section_m2 <= 0.0 { + return f64::INFINITY; + } + b_normal_t * pole_area_m2 / iron_section_m2 + } +} diff --git a/crates/vcad-kernel-magnetostatic/src/lib.rs b/crates/vcad-kernel-magnetostatic/src/lib.rs new file mode 100644 index 00000000..b778f04a --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/src/lib.rs @@ -0,0 +1,90 @@ +#![warn(missing_docs)] + +//! Grid-free 3D magnetostatics for **air-core** machines. +//! +//! This is the exact oracle behind air-core motor claims: given the actual coil +//! copper and rotor magnets, it computes flux linkage, back-EMF and torque +//! constants, inductance, and the torque waveform — with no closed-form fudge +//! factors and no fitted coefficients. +//! +//! # Why there is no grid +//! +//! An air-core machine has no iron. No iron means no `B`–`H` curve, so `μ = μ₀` +//! everywhere, the problem is **linear**, and superposition holds exactly. That +//! collapses the usual finite-element machinery into a sum over current +//! segments ([`filament`]): closed-form Biot-Savart for `B`, closed-form +//! `∮A·dl` for flux linkage. No mesh, no iteration, no convergence study, and no +//! truncation boundary to place — the three things that make +//! [`vcad_kernel_em`](https://docs.rs/vcad-kernel-em)'s grid solvers expensive +//! and approximate are simply absent. +//! +//! The price is that this crate is **only** valid without ferromagnetic +//! material. Anything with a steel back-iron, a slotted stator, or a saturable +//! pole must go to the grid solver instead. +//! +//! # Scope and honesty +//! +//! What this crate models exactly (to discretization): +//! +//! - Fields of arbitrary 3D conductor paths, including real spiral PCB copper. +//! - Permanent magnets as **equivalent surface currents** `K = M × n̂`, which is +//! exact for uniform magnetization in a linear medium with `μ_rec = 1` — the +//! right model for sintered ferrite and NdFeB, whose recoil permeability is +//! 1.05–1.1. That 5–10% is the model's dominant magnet error. +//! - Torque by two independent routes (Lorentz `∮I dl × B` and the co-energy +//! derivative `dW/dθ`), whose agreement is reported rather than assumed. +//! +//! What it does **not** model: +//! +//! - Any ferromagnetic material — see above. A steel back-iron behind the +//! magnets, which most axial-flux rotors have, roughly doubles airgap flux and +//! is **not** captured. Model it in the grid solver, or treat this crate's +//! result as the no-back-iron lower bound. +//! - Eddy currents, skin and proximity effect: statics only. DC resistance and +//! low-frequency inductance are in scope; AC losses are not. +//! - Self-inductance is filamentary, so it depends on the conductor radius used +//! to regularize the on-axis singularity. Rectangular PCB traces have no +//! single "radius"; we use the geometric-mean-distance equivalent, and the +//! residual error on `L` is larger than on `Kt` — see [`machine`]. +//! - Temperature. `Kt` shifts with magnet remanence (≈ −0.2%/K for NdFeB, +//! ≈ −0.2%/K for ferrite) and `R` with copper resistivity (+0.39%/K). +//! Callers supply the operating temperature; nothing here guesses it. +//! +//! # Cogging +//! +//! A true air-core machine has **no cogging torque**, because cogging is the +//! magnets detenting against ferromagnetic structure and there is none. This +//! crate therefore returns zero cogging *by construction*, not as a computed +//! result — do not read it as a validated prediction. What *is* computed is +//! torque ripple under load, which for an air-core machine comes from winding +//! and magnet-field harmonics alone; see [`kpi`]. +//! +//! # Units +//! +//! The public machine description is in **millimetres** (vcad convention); +//! everything internal is SI. Field accessors take and return SI. + +pub mod filament; +pub mod iron; +pub mod magnet; +pub mod vec3; + +pub use filament::{Filament, Segment}; +pub use iron::IronStack; +pub use magnet::{MagnetRing, Polarity, PrismMagnet}; +pub use vec3::Vec3; + +/// Vacuum permeability, H/m (CODATA 2018). +pub const MU_0: f64 = 1.256_637_062_12e-6; + +/// Copper resistivity at 20 °C, Ω·m. +pub const RHO_CU_20C: f64 = 1.68e-8; + +/// Copper temperature coefficient of resistivity, 1/K. +pub const ALPHA_CU: f64 = 0.00393; + +/// Millimetres to metres. +#[inline] +pub const fn mm(v: f64) -> f64 { + v * 1e-3 +} diff --git a/crates/vcad-kernel-magnetostatic/src/magnet.rs b/crates/vcad-kernel-magnetostatic/src/magnet.rs new file mode 100644 index 00000000..c3c1bcf3 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/src/magnet.rs @@ -0,0 +1,216 @@ +//! Permanent magnets as equivalent bound surface currents. +//! +//! A uniformly magnetized body is exactly equivalent to a current distribution: +//! `J_b = ∇×M` inside (zero for uniform `M`) and `K_b = M × n̂` on the surface. +//! For a prism magnetized along **z**, `M = M ẑ` gives `K = 0` on the flat faces +//! (`ẑ × ±ẑ = 0`) and `K = M (ẑ × n̂)` on the side walls — which, traced around +//! the footprint, is simply a **closed loop following the outline**. +//! +//! So any axially-magnetized prism — disc, arc segment, rectangle — becomes a +//! stack of closed loops around its footprint, each carrying `M·Δz`. One +//! representation covers every magnet shape a machine uses, and it feeds the +//! same segment integrator as the coils. +//! +//! # Fidelity +//! +//! Exact for uniform magnetization in a linear medium with recoil permeability +//! `μ_rec = 1`. Real sintered magnets have `μ_rec ≈ 1.05` (NdFeB) to `1.1` +//! (ferrite), so the model neglects the magnet's own reluctance — a 5–10% +//! over-prediction of working flux, and the dominant magnet-side error here. +//! Grade against a grid solve before quoting better than that. + +use serde::{Deserialize, Serialize}; + +use crate::filament::Filament; +use crate::vec3::Vec3; +use crate::MU_0; + +use std::f64::consts::PI; + +/// Which way a magnet's north pole faces along **z**. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)] +pub enum Polarity { + /// North toward +z. + North, + /// North toward −z. + South, +} + +impl Polarity { + /// `+1` for [`Polarity::North`], `−1` for [`Polarity::South`]. + pub fn sign(self) -> f64 { + match self { + Polarity::North => 1.0, + Polarity::South => -1.0, + } + } + + /// Alternating polarity for pole index `k`. + pub fn alternating(k: usize) -> Self { + if k % 2 == 0 { + Polarity::North + } else { + Polarity::South + } + } +} + +/// An axially-magnetized prism, described by its footprint outline. +/// +/// All lengths in **metres** (the machine layer converts from millimetres). +#[derive(Debug, Clone, PartialEq)] +pub struct PrismMagnet { + /// Footprint outline in the xy plane, m. Implicitly closed; winding order + /// does not matter, it is normalized to counter-clockwise. + pub footprint: Vec<(f64, f64)>, + /// Lower face z, m. + pub z0_m: f64, + /// Thickness along z, m. + pub thickness_m: f64, + /// Remanence `B_r`, tesla. + pub remanence_t: f64, + /// Which way north faces. + pub polarity: Polarity, +} + +impl PrismMagnet { + /// Magnetization `M = B_r/μ₀`, A/m. + pub fn magnetization(&self) -> f64 { + self.remanence_t / MU_0 * self.polarity.sign() + } + + /// Signed area of the footprint, m² — positive when counter-clockwise. + fn signed_area(&self) -> f64 { + let n = self.footprint.len(); + if n < 3 { + return 0.0; + } + let mut s = 0.0; + for i in 0..n { + let (x1, y1) = self.footprint[i]; + let (x2, y2) = self.footprint[(i + 1) % n]; + s += x1 * y2 - x2 * y1; + } + s / 2.0 + } + + /// Pole-face area, m². + pub fn face_area_m2(&self) -> f64 { + self.signed_area().abs() + } + + /// The equivalent bound-current loops, `n_axial` of them through the + /// thickness. + /// + /// Each loop carries `M · thickness/n_axial` amperes and traces the + /// footprint counter-clockwise, so a north-up magnet produces `+z` flux + /// above it — the right-hand rule, same as a coil. + pub fn to_filaments(&self, n_axial: usize) -> Vec { + let n_axial = n_axial.max(1); + if self.footprint.len() < 3 || self.thickness_m <= 0.0 { + return Vec::new(); + } + // Normalize to counter-clockwise so the current sense follows polarity + // rather than however the caller happened to order the points. + let ccw = self.signed_area() >= 0.0; + let dz = self.thickness_m / n_axial as f64; + let current = self.magnetization() * dz; + // The conductor radius regularizes the sheet: use half the axial slice + // height, which is the physical thickness each loop stands in for. + let wire_r = dz * 0.5; + + (0..n_axial) + .map(|k| { + let z = self.z0_m + (k as f64 + 0.5) * dz; + let mut pts: Vec = self + .footprint + .iter() + .map(|&(x, y)| Vec3::new(x, y, z)) + .collect(); + if !ccw { + pts.reverse(); + } + Filament::closed_loop(pts, current, wire_r) + }) + .collect() + } +} + +/// A ring of identical magnets on a pitch circle with alternating polarity — +/// the rotor of a surface-magnet machine. +#[derive(Debug, Clone, PartialEq)] +pub struct MagnetRing { + /// The magnets, already positioned. + pub magnets: Vec, +} + +impl MagnetRing { + /// Discs of `diameter_m` centred on a `pitch_radius_m` circle, `poles` of + /// them, alternating N/S. + /// + /// `poles` must be even for the ring to be magnetically balanced — an odd + /// count leaves a net dipole moment, which breaks the two-plane image series + /// (see [`crate::IronStack`]). + #[allow(clippy::too_many_arguments)] + pub fn discs( + poles: usize, + pitch_radius_m: f64, + diameter_m: f64, + z0_m: f64, + thickness_m: f64, + remanence_t: f64, + facets: usize, + ) -> Self { + let facets = facets.max(12); + let r = diameter_m * 0.5; + let magnets = (0..poles) + .map(|k| { + let phi = 2.0 * PI * (k as f64) / (poles as f64); + let (cx, cy) = (pitch_radius_m * phi.cos(), pitch_radius_m * phi.sin()); + let footprint = (0..facets) + .map(|j| { + let t = 2.0 * PI * (j as f64) / (facets as f64); + (cx + r * t.cos(), cy + r * t.sin()) + }) + .collect(); + PrismMagnet { + footprint, + z0_m, + thickness_m, + remanence_t, + polarity: Polarity::alternating(k), + } + }) + .collect(); + Self { magnets } + } + + /// All bound-current loops for the ring. + pub fn to_filaments(&self, n_axial: usize) -> Vec { + self.magnets.iter().flat_map(|m| m.to_filaments(n_axial)).collect() + } + + /// The ring rotated about z by `angle` radians — the rotor position sweep. + pub fn rotated_z(&self, angle: f64) -> MagnetRing { + let (s, c) = angle.sin_cos(); + MagnetRing { + magnets: self + .magnets + .iter() + .map(|m| PrismMagnet { + footprint: m + .footprint + .iter() + .map(|&(x, y)| (c * x - s * y, s * x + c * y)) + .collect(), + ..m.clone() + }) + .collect(), + } + } + + /// Total pole-face area, m². + pub fn face_area_m2(&self) -> f64 { + self.magnets.iter().map(|m| m.face_area_m2()).sum() + } +} diff --git a/crates/vcad-kernel-magnetostatic/src/vec3.rs b/crates/vcad-kernel-magnetostatic/src/vec3.rs new file mode 100644 index 00000000..4d42c655 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/src/vec3.rs @@ -0,0 +1,131 @@ +//! Minimal 3-vector. SI metres throughout the solver interior; the public +//! machine description takes millimetres per the vcad convention and converts +//! once at the boundary. + +use serde::{Deserialize, Serialize}; + +/// A 3-vector of `f64`. +#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)] +pub struct Vec3 { + /// x component. + pub x: f64, + /// y component. + pub y: f64, + /// z component. + pub z: f64, +} + +impl Vec3 { + /// The zero vector. + pub const ZERO: Vec3 = Vec3 { x: 0.0, y: 0.0, z: 0.0 }; + + /// Construct from components. + #[inline] + pub const fn new(x: f64, y: f64, z: f64) -> Self { + Self { x, y, z } + } + + /// A point on the z axis. + #[inline] + pub const fn z_axis(z: f64) -> Self { + Self { x: 0.0, y: 0.0, z } + } + + /// Cylindrical construction: radius, azimuth (rad), height. + #[inline] + pub fn cylindrical(r: f64, theta: f64, z: f64) -> Self { + Self { x: r * theta.cos(), y: r * theta.sin(), z } + } + + /// Euclidean length. + #[inline] + pub fn norm(self) -> f64 { + self.dot(self).sqrt() + } + + /// Squared length — avoids a `sqrt` where only comparisons matter. + #[inline] + pub fn norm_sq(self) -> f64 { + self.dot(self) + } + + /// Dot product. + #[inline] + pub fn dot(self, o: Vec3) -> f64 { + self.x * o.x + self.y * o.y + self.z * o.z + } + + /// Cross product. + #[inline] + pub fn cross(self, o: Vec3) -> Vec3 { + Vec3 { + x: self.y * o.z - self.z * o.y, + y: self.z * o.x - self.x * o.z, + z: self.x * o.y - self.y * o.x, + } + } + + /// Unit vector. Returns [`Vec3::ZERO`] for a zero-length input rather than + /// NaN — callers in this crate always guard degenerate geometry first, and + /// a silent NaN is far harder to trace than a silent zero. + #[inline] + pub fn normalized(self) -> Vec3 { + let n = self.norm(); + if n <= 0.0 { + Vec3::ZERO + } else { + self * (1.0 / n) + } + } + + /// Rotation about the +z axis by `angle` radians. + #[inline] + pub fn rotated_z(self, angle: f64) -> Vec3 { + let (s, c) = angle.sin_cos(); + Vec3 { x: c * self.x - s * self.y, y: s * self.x + c * self.y, z: self.z } + } + + /// Distance from the z axis. + #[inline] + pub fn radius(self) -> f64 { + self.x.hypot(self.y) + } +} + +impl std::ops::Add for Vec3 { + type Output = Vec3; + #[inline] + fn add(self, o: Vec3) -> Vec3 { + Vec3 { x: self.x + o.x, y: self.y + o.y, z: self.z + o.z } + } +} + +impl std::ops::Sub for Vec3 { + type Output = Vec3; + #[inline] + fn sub(self, o: Vec3) -> Vec3 { + Vec3 { x: self.x - o.x, y: self.y - o.y, z: self.z - o.z } + } +} + +impl std::ops::Mul for Vec3 { + type Output = Vec3; + #[inline] + fn mul(self, s: f64) -> Vec3 { + Vec3 { x: self.x * s, y: self.y * s, z: self.z * s } + } +} + +impl std::ops::Neg for Vec3 { + type Output = Vec3; + #[inline] + fn neg(self) -> Vec3 { + Vec3 { x: -self.x, y: -self.y, z: -self.z } + } +} + +impl std::iter::Sum for Vec3 { + fn sum>(iter: I) -> Vec3 { + iter.fold(Vec3::ZERO, |a, b| a + b) + } +} diff --git a/crates/vcad-kernel-magnetostatic/tests/filament_validation.rs b/crates/vcad-kernel-magnetostatic/tests/filament_validation.rs new file mode 100644 index 00000000..a360b493 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/tests/filament_validation.rs @@ -0,0 +1,193 @@ +//! Validation ladder for the segment integrator. +//! +//! Every claim this crate makes rests on [`Segment::b_at`] and [`Segment::a_at`] +//! being right, so they are checked against three independent things: a closed +//! form derived by hand (infinite wire), an exact special function from a +//! *different codebase* (`vcad_kernel_particle::field::b_ring`, elliptic +//! integrals written for particle optics), and an internal consistency identity +//! (curl of A equals B, by finite differences). + +use std::f64::consts::PI; + +use vcad_kernel_magnetostatic::{Filament, Segment, Vec3, MU_0}; +use vcad_kernel_particle::field::{b_ring, RingCoil}; + +/// A regular polygon approximating a circular loop in the z = `z0` plane. +fn ring_filament(radius_m: f64, z0: f64, current_a: f64, n: usize) -> Filament { + let pts = (0..n) + .map(|i| { + let t = 2.0 * PI * (i as f64) / (n as f64); + Vec3::cylindrical(radius_m, t, z0) + }) + .collect(); + Filament::closed_loop(pts, current_a, 0.0) +} + +#[test] +fn long_segment_reproduces_the_infinite_wire() { + // B = μ₀I/(2πd) for an infinite wire; a segment ±L with L ≫ d approaches it. + let i = 3.0; + let d = 0.01; + let l = 1.0e4 * d; + let seg = Segment { + a: Vec3::new(0.0, 0.0, -l), + b: Vec3::new(0.0, 0.0, l), + current_a: i, + wire_radius_m: 0.0, + }; + let b = seg.b_at(Vec3::new(d, 0.0, 0.0)); + let expected = MU_0 * i / (2.0 * PI * d); + assert!( + (b.norm() - expected).abs() / expected < 1e-6, + "magnitude {} vs {}", + b.norm(), + expected + ); + // Right-hand rule: current +z, field point +x ⇒ B along +y. + assert!(b.y > 0.0 && b.x.abs() < 1e-12 * expected && b.z.abs() < 1e-12 * expected); +} + +#[test] +fn on_axis_loop_matches_the_closed_form() { + // B_z = μ₀ I R² / (2 (R² + z²)^{3/2}) on the axis of a circular loop. + let r = 0.02; + let i = 5.0; + let loop_ = ring_filament(r, 0.0, i, 2048); + for z in [0.0, 0.005, 0.02, 0.05] { + let b = loop_.b_at(Vec3::z_axis(z)); + let expected = MU_0 * i * r * r / (2.0 * (r * r + z * z).powf(1.5)); + assert!( + (b.z - expected).abs() / expected < 1e-5, + "z={z}: {} vs {}", + b.z, + expected + ); + assert!(b.x.abs().max(b.y.abs()) < 1e-9 * expected, "off-axis leakage at z={z}"); + } +} + +#[test] +fn off_axis_loop_matches_particle_crates_elliptic_integrals() { + // Independent-codebase check: `b_ring` solves the same loop with complete + // elliptic integrals, written for a different crate and a different purpose. + let r = 0.02; + let i = 5.0; + let coil = RingCoil { radius_m: r, z_m: 0.0, ampere_turns: i, wire_radius_m: 0.0 }; + let loop_ = ring_filament(r, 0.0, i, 4096); + + for &(rho, z) in &[ + (0.005, 0.0), + (0.010, 0.004), + (0.019, 0.010), + (0.030, 0.006), + (0.050, 0.030), + (0.002, -0.008), + ] { + let (br_ref, bz_ref) = b_ring(&coil, rho, z); + let p = Vec3::new(rho, 0.0, z); + let b = loop_.b_at(p); + // At y = 0 the radial direction is +x. + let br = b.x; + let bz = b.z; + let scale = (br_ref * br_ref + bz_ref * bz_ref).sqrt().max(1e-12); + assert!( + (br - br_ref).abs() / scale < 2e-4, + "B_r at (ρ={rho}, z={z}): {br} vs {br_ref}" + ); + assert!( + (bz - bz_ref).abs() / scale < 2e-4, + "B_z at (ρ={rho}, z={z}): {bz} vs {bz_ref}" + ); + } +} + +#[test] +fn curl_of_a_equals_b() { + // The two integrals are derived and implemented separately; ∇×A = B ties + // them together and would catch a sign or factor error in either. + let loop_ = ring_filament(0.02, 0.0, 4.0, 1024); + let h = 1e-6; + let curl_at = |p: Vec3| -> Vec3 { + let d = |axis: usize, q: Vec3| -> Vec3 { + let mut plus = q; + let mut minus = q; + match axis { + 0 => { + plus.x += h; + minus.x -= h; + } + 1 => { + plus.y += h; + minus.y -= h; + } + _ => { + plus.z += h; + minus.z -= h; + } + } + (loop_.a_at(plus) - loop_.a_at(minus)) * (1.0 / (2.0 * h)) + }; + let dx = d(0, p); + let dy = d(1, p); + let dz = d(2, p); + Vec3::new(dy.z - dz.y, dz.x - dx.z, dx.y - dy.x) + }; + + for p in [ + Vec3::new(0.004, 0.003, 0.006), + Vec3::new(0.0, 0.0, 0.01), + Vec3::new(0.035, -0.01, 0.004), + ] { + let b = loop_.b_at(p); + let c = curl_at(p); + let scale = b.norm().max(1e-9); + assert!( + (c - b).norm() / scale < 1e-4, + "curl A = {c:?} vs B = {b:?} at {p:?}" + ); + } +} + +#[test] +fn mutual_inductance_matches_maxwells_coaxial_formula() { + // Two coaxial loops. λ = ∮A·dl from one loop through the other gives M, + // which exercises `a_at` and `flux_linkage` on a case with a known answer. + // Maxwell: M = μ₀√(R₁R₂)[(2/k − k)K(k²) − (2/k)E(k²)], + // k² = 4R₁R₂/((R₁+R₂)² + d²). Checked here against the small-k dipole limit, + // where M → μ₀πR₁²R₂²/(2 d³) for d ≫ R. + let r1 = 0.01; + let r2 = 0.012; + let d = 0.5; // far apart ⇒ dipole limit is accurate + let source = ring_filament(r1, 0.0, 1.0, 2048); + let sensor = ring_filament(r2, d, 1.0, 2048); + + let m = sensor.flux_linkage(|p| source.a_at(p)); + let dipole = MU_0 * PI * r1 * r1 * r2 * r2 / (2.0 * d.powi(3)); + assert!( + (m - dipole).abs() / dipole < 2e-3, + "M = {m} vs dipole limit {dipole}" + ); +} + +#[test] +fn wire_regularization_keeps_the_field_finite_and_odd() { + // Approaching the conductor from either side must stay finite and reverse + // sign through the centreline — the uniform-current-density model. + let a_r = 1e-4; + let seg = Segment { + a: Vec3::new(0.0, 0.0, -1.0), + b: Vec3::new(0.0, 0.0, 1.0), + current_a: 10.0, + wire_radius_m: a_r, + }; + let inside = seg.b_at(Vec3::new(a_r * 0.5, 0.0, 0.0)); + let surface = seg.b_at(Vec3::new(a_r, 0.0, 0.0)); + assert!(inside.norm().is_finite() && surface.norm().is_finite()); + // Linear ramp to the surface value. + assert!((inside.norm() / surface.norm() - 0.5).abs() < 1e-6); + // Exactly on the axis the field vanishes by symmetry. + assert!(seg.b_at(Vec3::ZERO).norm() < 1e-18); + // Odd through the centreline. + let other = seg.b_at(Vec3::new(-a_r * 0.5, 0.0, 0.0)); + assert!((other.y + inside.y).abs() < 1e-12 * inside.norm().max(1e-12)); +} diff --git a/crates/vcad-kernel-magnetostatic/tests/iron_images.rs b/crates/vcad-kernel-magnetostatic/tests/iron_images.rs new file mode 100644 index 00000000..3987fefa --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/tests/iron_images.rs @@ -0,0 +1,194 @@ +//! The image model must reproduce the `μ→∞` boundary condition it claims. +//! +//! These are the checks that decide whether back-iron can be handled without a +//! grid: at an infinitely permeable surface the tangential field must vanish and +//! the normal field must double, and the two-plane image series must converge. + +use std::f64::consts::PI; + +use vcad_kernel_magnetostatic::{Filament, IronStack, Vec3}; + +fn ring(radius_m: f64, z0: f64, current_a: f64, n: usize) -> Filament { + let pts = (0..n) + .map(|i| Vec3::cylindrical(radius_m, 2.0 * PI * (i as f64) / (n as f64), z0)) + .collect(); + Filament::closed_loop(pts, current_a, 0.0) +} + +/// Sum B over a source and its images. +fn b_with_iron(stack: &IronStack, src: &Filament, p: Vec3) -> Vec3 { + stack.expand(src).iter().map(|f| f.b_at(p)).sum() +} + +#[test] +fn a_single_plane_doubles_normal_field_and_cancels_tangential() { + // The defining property of a μ→∞ mirror. A loop at z = h above a plane at + // z = 0: on the plane, B_z must be exactly twice the free-space value and + // the radial component must vanish. + let h = 0.006; + let loop_ = ring(0.02, h, 3.0, 1024); + let stack = IronStack::single(0.0); + + for &rho in &[0.0, 0.005, 0.015, 0.025, 0.04] { + let p = Vec3::new(rho, 0.0, 0.0); + let free = loop_.b_at(p); + let with = b_with_iron(&stack, &loop_, p); + assert!( + (with.z - 2.0 * free.z).abs() <= 1e-9 * free.z.abs().max(1e-12), + "ρ={rho}: B_z {} vs 2×{}", + with.z, + free.z + ); + // Tangential (radial, here +x) component cancels at the surface. + let scale = with.z.abs().max(1e-12); + assert!( + with.x.abs() / scale < 1e-9, + "ρ={rho}: tangential B_x = {} did not cancel", + with.x + ); + } +} + +#[test] +fn image_of_a_loop_matches_a_hand_placed_mirror_loop() { + // Independent construction: build the mirror loop explicitly and compare. + let h = 0.004; + let src = ring(0.015, h, 2.5, 512); + let mirror = ring(0.015, -h, 2.5, 512); + let stack = IronStack::single(0.0); + + for p in [ + Vec3::new(0.0, 0.0, 0.01), + Vec3::new(0.012, 0.004, 0.006), + Vec3::new(0.03, -0.01, 0.002), + ] { + let by_images = b_with_iron(&stack, &src, p); + let by_hand = src.b_at(p) + mirror.b_at(p); + let scale = by_hand.norm().max(1e-15); + assert!( + (by_images - by_hand).norm() / scale < 1e-12, + "at {p:?}: {by_images:?} vs {by_hand:?}" + ); + } +} + +#[test] +fn back_iron_raises_airgap_flux_by_the_expected_factor() { + // The engineering claim that motivated all of this: a back-iron behind the + // magnet is worth roughly a factor of two in airgap flux. Modelled here with + // a current loop standing in for the magnet's bound surface current. + let magnet_z = 0.003; + let src = ring(0.0075, magnet_z, 100.0, 512); + let probe = Vec3::z_axis(0.0075); // in the airgap, beyond the magnet + + let free = src.b_at(probe).z; + // Iron immediately behind the magnet. + let backed = b_with_iron(&IronStack::single(0.0), &src, probe).z; + + assert!(backed > free, "back-iron must increase airgap flux"); + let ratio = backed / free; + assert!( + (1.3..2.0).contains(&ratio), + "back-iron gain {ratio} outside the physically expected 1.3–2.0 band" + ); +} + +/// Six alternating poles on a 22.5 mm pitch circle — a magnetically balanced +/// rotor, which is what every real multi-pole machine is. +fn alternating_rotor() -> Vec { + (0..6) + .map(|k| { + let t = 2.0 * PI * (k as f64) / 6.0; + let c = Vec3::new(0.0225 * t.cos(), 0.0225 * t.sin(), 0.003); + let sign = if k % 2 == 0 { 1.0 } else { -1.0 }; + let pts = (0..256) + .map(|j| { + let u = 2.0 * PI * (j as f64) / 256.0; + Vec3::new(c.x + 0.0075 * u.cos(), c.y + 0.0075 * u.sin(), c.z) + }) + .collect(); + Filament::closed_loop(pts, 100.0 * sign, 0.0) + }) + .collect() +} + +#[test] +fn balance_residual_separates_a_rotor_from_a_lone_magnet() { + let rotor = alternating_rotor(); + assert!( + IronStack::balance_residual(&rotor) < 1e-9, + "an alternating rotor must be balanced, got {}", + IronStack::balance_residual(&rotor) + ); + let lone = vec![ring(0.0075, 0.003, 100.0, 256)]; + assert!( + (IronStack::balance_residual(&lone) - 1.0).abs() < 1e-9, + "a lone loop must be maximally unbalanced" + ); +} + +#[test] +fn two_plane_series_converges_for_a_balanced_rotor() { + // Facing mirrors generate an infinite series. It converges — and is only + // usable — when the source set carries no net moment. + let rotor = alternating_rotor(); + let probe = Vec3::new(0.0225, 0.0, 0.0045); + let b_at = |depth: usize| -> Vec3 { + let stack = IronStack::pair(0.0, 0.0056, depth); + rotor.iter().flat_map(|s| stack.expand(s)).map(|f| f.b_at(probe)).sum() + }; + let shallow = b_at(8); + let deep = b_at(IronStack::DEFAULT_REFLECTIONS); + let deeper = b_at(32); + + // Converged: doubling the depth past the default barely moves it. + assert!( + (deeper - deep).norm() / deeper.norm() < 1e-4, + "not converged at default depth: {} vs {}", + deep.norm(), + deeper.norm() + ); + // And the approach is monotone in depth, not oscillating. + assert!((deep - shallow).norm() > (deeper - deep).norm()); +} + +#[test] +fn an_unbalanced_source_is_flagged_as_non_convergent() { + // The failure mode the precondition exists to catch: a single magnet between + // two infinite plates never settles, and `tail_fraction` must say so loudly + // rather than returning a plausible number. + let lone = ring(0.0075, 0.003, 100.0, 256); + let probe = Vec3::new(0.018, 0.0, 0.0045); + let stack = IronStack::pair(0.0, 0.0056, IronStack::DEFAULT_REFLECTIONS); + let tail = stack.tail_fraction(&lone, probe); + assert!( + tail > 0.05, + "an unbalanced source must report a large truncation residual, got {tail}" + ); +} + +#[test] +fn two_planes_beat_one_which_beats_none() { + // Monotone ordering is the sanity check on the whole iron model: closing + // more of the magnetic circuit can only raise the working flux. + let src = ring(0.02, 0.003, 5.0, 256); + let probe = Vec3::new(0.02, 0.0, 0.0045); + + let none = b_with_iron(&IronStack::none(), &src, probe).z.abs(); + let one = b_with_iron(&IronStack::single(0.0), &src, probe).z.abs(); + let two = b_with_iron(&IronStack::pair(0.0, 0.0056, 12), &src, probe).z.abs(); + + assert!(none < one, "one plane must beat none: {none} vs {one}"); + assert!(one < two, "two planes must beat one: {one} vs {two}"); +} + +#[test] +fn no_iron_is_a_pure_passthrough() { + let src = ring(0.02, 0.003, 5.0, 128); + let stack = IronStack::none(); + assert_eq!(stack.expand(&src).len(), 1); + assert!(stack.is_empty()); + let p = Vec3::new(0.01, 0.002, 0.006); + assert_eq!(b_with_iron(&stack, &src, p), src.b_at(p)); + assert_eq!(stack.tail_fraction(&src, p), 0.0); +} diff --git a/crates/vcad-kernel-magnetostatic/tests/magnet_model.rs b/crates/vcad-kernel-magnetostatic/tests/magnet_model.rs new file mode 100644 index 00000000..03be59c4 --- /dev/null +++ b/crates/vcad-kernel-magnetostatic/tests/magnet_model.rs @@ -0,0 +1,131 @@ +//! The equivalent-surface-current magnet model, checked against the closed form +//! for an axially magnetized cylinder. + +use std::f64::consts::PI; + +use vcad_kernel_magnetostatic::{IronStack, MagnetRing, Polarity, PrismMagnet, Vec3}; + +fn disc(radius_m: f64, z0: f64, t: f64, br: f64, pol: Polarity, facets: usize) -> PrismMagnet { + PrismMagnet { + footprint: (0..facets) + .map(|j| { + let a = 2.0 * PI * (j as f64) / (facets as f64); + (radius_m * a.cos(), radius_m * a.sin()) + }) + .collect(), + z0_m: z0, + thickness_m: t, + remanence_t: br, + polarity: pol, + } +} + +/// On-axis field of an axially magnetized cylinder, the standard closed form: +/// `B_z = (B_r/2)[ (z+L)/√((z+L)²+R²) − z/√(z²+R²) ]`, with `z` measured from +/// the near face. +fn cylinder_axial_closed_form(br: f64, radius: f64, length: f64, z_from_face: f64) -> f64 { + let a = (z_from_face + length) / ((z_from_face + length).powi(2) + radius * radius).sqrt(); + let b = z_from_face / (z_from_face * z_from_face + radius * radius).sqrt(); + 0.5 * br * (a - b) +} + +#[test] +fn disc_magnet_matches_the_closed_form_on_axis() { + // Y30 ferrite, the reference motor's magnet: D15 x 3 mm, Br 385 mT. + let br = 0.385; + let r = 0.0075; + let t = 0.003; + let m = disc(r, 0.0, t, br, Polarity::North, 720); + let fils = m.to_filaments(48); + + for z_face in [0.0005, 0.001, 0.002, 0.004, 0.008] { + let p = Vec3::z_axis(t + z_face); // above the top face + let b: Vec3 = fils.iter().map(|f| f.b_at(p)).sum(); + let want = cylinder_axial_closed_form(br, r, t, z_face); + assert!( + (b.z - want).abs() / want.abs() < 0.02, + "z={z_face}: {} vs closed form {}", + b.z, + want + ); + assert!(b.x.abs().max(b.y.abs()) < 1e-6 * want.abs(), "off-axis leakage"); + } +} + +#[test] +fn the_reference_magnet_gives_the_free_space_airgap_number() { + // The number quoted when choosing the solver path: a Y30 D15x3 disc produces + // ~65 mT at 1 mm from its face in free space. This is the "no iron" anchor + // that the back-iron gain is measured against. + let m = disc(0.0075, 0.0, 0.003, 0.385, Polarity::North, 720); + let fils = m.to_filaments(48); + let b: Vec3 = fils.iter().map(|f| f.b_at(Vec3::z_axis(0.004))).sum(); + assert!( + (0.060..0.070).contains(&b.z), + "free-space airgap flux {} T outside the expected 60–70 mT", + b.z + ); +} + +#[test] +fn polarity_flips_the_field() { + let n = disc(0.0075, 0.0, 0.003, 0.385, Polarity::North, 180); + let s = disc(0.0075, 0.0, 0.003, 0.385, Polarity::South, 180); + let p = Vec3::z_axis(0.005); + let bn: Vec3 = n.to_filaments(8).iter().map(|f| f.b_at(p)).sum(); + let bs: Vec3 = s.to_filaments(8).iter().map(|f| f.b_at(p)).sum(); + assert!(bn.z > 0.0, "north-up must push +z flux above the magnet"); + assert!((bn.z + bs.z).abs() < 1e-12 * bn.z.abs()); +} + +#[test] +fn footprint_winding_order_does_not_change_the_field() { + // Callers should not have to know the orientation convention. + let ccw = disc(0.0075, 0.0, 0.003, 0.385, Polarity::North, 120); + let mut cw = ccw.clone(); + cw.footprint.reverse(); + let p = Vec3::new(0.002, 0.001, 0.005); + let a: Vec3 = ccw.to_filaments(6).iter().map(|f| f.b_at(p)).sum(); + let b: Vec3 = cw.to_filaments(6).iter().map(|f| f.b_at(p)).sum(); + assert!((a - b).norm() / a.norm() < 1e-12); +} + +#[test] +fn an_alternating_rotor_ring_is_magnetically_balanced() { + // The precondition the two-plane image series depends on. + let ring = MagnetRing::discs(6, 0.0225, 0.015, 0.0, 0.003, 0.385, 64); + let fils = ring.to_filaments(4); + assert_eq!(ring.magnets.len(), 6); + assert!( + IronStack::balance_residual(&fils) < 1e-9, + "alternating ring must be balanced, got {}", + IronStack::balance_residual(&fils) + ); +} + +#[test] +fn back_iron_raises_the_reference_airgap_flux_toward_the_mec_value() { + // The decision that picked this solver path: for the reference geometry the + // steel is worth ~2.4x, taking ~65 mT of free-space flux toward the ~155 mT + // the reluctance model predicts for the full circuit. Checked here as a + // gain band rather than a point value — the MEC number is itself only + // first-order, and finite-disc fringing is not modelled. + let ring = MagnetRing::discs(6, 0.0225, 0.015, 0.0, 0.003, 0.385, 96); + let fils = ring.to_filaments(16); + let probe = Vec3::new(0.0225, 0.0, 0.004); // 1 mm above a pole face + + let free: Vec3 = fils.iter().map(|f| f.b_at(probe)).sum(); + // Rotor back-iron sits directly behind the magnets (z = 0). + let backed_stack = IronStack::single(0.0); + let backed: Vec3 = fils + .iter() + .flat_map(|f| backed_stack.expand(f)) + .map(|f| f.b_at(probe)) + .sum(); + + let gain = backed.z.abs() / free.z.abs(); + assert!( + (1.4..2.1).contains(&gain), + "single back-iron gain {gain} outside the expected 1.4–2.1 band" + ); +} diff --git a/packages/mcp/src/__tests__/ecad.test.ts b/packages/mcp/src/__tests__/ecad.test.ts index 04f0338e..028e62e4 100644 --- a/packages/mcp/src/__tests__/ecad.test.ts +++ b/packages/mcp/src/__tests__/ecad.test.ts @@ -3536,7 +3536,7 @@ describe("add_coil lead-out and multilayer", () => { const fcu = b.traces.filter((t) => t.net === "PHA" && t.layer === "FCu").length; const bcu = b.traces.filter((t) => t.net === "PHA" && t.layer === "BCu").length; expect(fcu).toBeGreaterThan(0); - expect(bcu).toBe(fcu); // same geometry on both layers + expect(bcu).toBe(fcu); // one segment per layer per spiral sample const stitch = b.vias.find((v) => v.net === "PHA"); expect(stitch).toBeDefined(); expect(stitch!.startLayer).toBe("FCu"); @@ -3544,6 +3544,113 @@ describe("add_coil lead-out and multilayer", () => { // Total length is the sum across both layers. expect(res.total_length_mm).toBeCloseTo(res.total_length_mm, 3); }); + + /** + * Magnetic dipole moment of the series current path, in mm² (m_z / I). + * + * Chains each layer's traces into a polyline, orients it along the direction + * the series current actually flows (entering at `entry`), and accumulates the + * shoelace sum about `center`. A stacked coil is only useful if consecutive + * layers circulate the *same* way — if they alternate, their axial fields + * subtract and the stack does nothing. + */ + function seriesDipoleMoment( + traces: Array<{ start: Vec2; end: Vec2; layer: string; net: string }>, + vias: Array<{ position: Vec2; startLayer: string; endLayer: string; net: string }>, + net: string, + layers: string[], + center: Vec2, + entry: Vec2, + ): { perLayer: number[]; total: number } { + const near = (a: Vec2, b: Vec2) => Math.hypot(a.x - b.x, a.y - b.y) < 1e-6; + const shoelace = (pts: Vec2[]) => { + let s = 0; + for (let i = 0; i + 1 < pts.length; i++) { + const x1 = pts[i].x - center.x, + y1 = pts[i].y - center.y; + const x2 = pts[i + 1].x - center.x, + y2 = pts[i + 1].y - center.y; + s += x1 * y2 - x2 * y1; + } + return s / 2; + }; + + const perLayer: number[] = []; + let cursor = entry; + for (let li = 0; li < layers.length; li++) { + const segs = traces.filter((t) => t.net === net && t.layer === layers[li]); + expect(segs.length).toBeGreaterThan(0); + // Traces are pushed contiguously, so the stored polyline is start→end. + const stored: Vec2[] = [segs[0].start, ...segs.map((s) => s.end)]; + const head = stored[0]; + const tail = stored[stored.length - 1]; + // Current enters this layer at whichever endpoint the previous hop landed on. + const forward = near(head, cursor); + expect(forward || near(tail, cursor)).toBe(true); + const path = forward ? stored : [...stored].reverse(); + perLayer.push(shoelace(path)); + cursor = path[path.length - 1]; + // Follow the stitch via to the next layer (it sits on the exit terminal). + if (li + 1 < layers.length) { + const v = vias.find( + (v) => + v.net === net && + v.startLayer === layers[li] && + v.endLayer === layers[li + 1] && + near(v.position, cursor), + ); + expect(v, `stitch via from ${layers[li]} to ${layers[li + 1]} at the exit terminal`).toBeDefined(); + } + } + return { perLayer, total: perLayer.reduce((a, b) => a + b, 0) }; + } + + // Top-down physical order. add_coil's COPPER_LAYERS tops out at 8 even though + // the IR enum carries In1Cu–In8Cu (10 copper layers). + const STACK_8 = ["FCu", "In1Cu", "In2Cu", "In3Cu", "In4Cu", "In5Cu", "In6Cu", "BCu"]; + for (const n of [2, 4, 8]) { + it(`a ${n}-layer stack adds its fields rather than cancelling, and fabricates`, async () => { + const id = await circleBoardSession(); + const center = { x: 40, y: 40 }; + const layers = n === 8 ? STACK_8 : [...STACK_8.slice(0, n - 1), "BCu"]; + const res = out( + await addCoil({ + document_id: id, + center, + turns: 3, + inner_radius: 5, + outer_radius: 12, + trace_width: 0.4, + clearance: 0.3, + net: "PHA", + layers, + }), + ); + expect(res.success).toBe(true); + // A stack that cancels is useless; so is one that will not fabricate. + // (Turnarounds recur at the same radius every two layers, so an integer + // `turns` used to drop two vias in one hole.) + expect(res.drc_delta.clean).toBe(true); + expect(res.stitch_vias.length).toBe(n - 1); + + const b = getPcbBoard(getSession(id)); + const { perLayer, total } = seriesDipoleMoment( + b.traces as never, + b.vias as never, + "PHA", + layers, + center, + res.terminals.a as Vec2, + ); + + // Every layer must contribute with the same sign... + const s0 = Math.sign(perLayer[0]); + for (const m of perLayer) expect(Math.sign(m)).toBe(s0); + // ...so the stack's moment scales with layer count instead of cancelling. + const single = Math.abs(perLayer[0]); + expect(Math.abs(total)).toBeGreaterThan(0.9 * n * single); + }); + } }); describe("add_motor_winding", () => { diff --git a/packages/mcp/src/tools/ecad.ts b/packages/mcp/src/tools/ecad.ts index c125f3da..f0f8555f 100644 --- a/packages/mcp/src/tools/ecad.ts +++ b/packages/mcp/src/tools/ecad.ts @@ -6984,16 +6984,91 @@ export async function addCoil( if ("error" in lr) return fail(`layers entry: ${lr.error}`); } const stitchVias: Array<{ position: Vec2; startLayer: PcbLayer; endLayer: PcbLayer }> = []; - const perLayerLen = segLength(pts); + + // The stack is ONE continuous spiral in θ: the angle marches monotonically + // across every layer while the radius zigzags inner↔outer, with a stitch via + // at each turnaround. That constant angular sweep is the whole point — series + // current then circulates the same way on every layer, so their axial fields + // ADD. (Reusing one spiral on all layers instead makes current run outer→inner + // on even layers and inner→outer on odd ones through identical copper, which + // reverses the circulation each layer and cancels the stack out.) + // + // Consecutive layers also meet at exactly the same (r, θ) by construction, so + // the stitch vias land on a shared point without any extra alignment work. + const dTheta = sign * turns * 2 * Math.PI; + const spiralSteps = Math.max(2, Math.ceil(turns * segmentsPerTurn)); + + // Turnaround angles, staggered. Turnarounds at the same radius recur every + // two layers, and for an integer `turns` they would otherwise land on the + // identical spoke — stacking two vias in one hole. Fan each successive + // turnaround by `stagger` so same-radius holes clear each other; the binding + // case is the inner radius, where a given angle buys the least arc. + const drill = pcb.rules.defaultRules.viaDrill; + const holePitch = drill + 0.5 + 0.05; // hole-to-hole minimum, plus margin + const staggerR = Math.max(innerR, traceWidth); + const stagger = Math.min(Math.PI / 4, holePitch / (2 * staggerR)); + // Boundary li sits between layer li and li+1; boundary -1 is the free start. + const boundary = (li: number) => + li < 0 ? theta0 + dTheta : theta0 - li * dTheta - sign * (li + 1) * stagger; + + const layerPts: Vec2[][] = []; + for (let li = 0; li < layersIn.length; li++) { + // Every layer sweeps the same way — that is what makes the fields add. + // Even layers run outer→inner, odd layers inner→outer. + const thetaFrom = boundary(li - 1); + const thetaTo = boundary(li); + const outward = li % 2 === 1; + const rFrom = outward ? innerR : outerR; + const rTo = outward ? outerR : innerR; + const lp: Vec2[] = []; + for (let s = 0; s <= spiralSteps; s++) { + const t = s / spiralSteps; + const theta = thetaFrom + t * (thetaTo - thetaFrom); + const r = rFrom + t * (rTo - rFrom); + const p = { + x: round3(center.x + r * Math.cos(theta)), + y: round3(center.y + r * Math.sin(theta)), + }; + const prev = lp[lp.length - 1]; + if (!prev || prev.x !== p.x || prev.y !== p.y) lp.push(p); + } + layerPts.push(lp); + } + + // Tangential lead-out at each inner turnaround: with integer `turns` the + // inner and outer terminals share a spoke, so the stitch via would otherwise + // sit at the same angle as the outer endpoint (a same-net bypass-short + // hazard). Displace both sides of the turnaround to the same new point so the + // via still lands on shared copper. + if (innerLeadOut > 0) { + for (let li = 0; li + 1 < layersIn.length; li++) { + if (li % 2 !== 0) continue; // odd→even turnarounds are at the outer radius + const end = layerPts[li][layerPts[li].length - 1]; + const thetaEnd = boundary(li); + // Continue along the direction of travel (the sweep runs by −dθ). + const tx = Math.sin(thetaEnd) * sign; + const ty = -Math.cos(thetaEnd) * sign; + const T = { + x: round3(end.x + tx * innerLeadOut), + y: round3(end.y + ty * innerLeadOut), + }; + if (T.x !== end.x || T.y !== end.y) { + layerPts[li].push(T); + layerPts[li + 1].unshift(T); + } + } + } + let totalLengthMm = 0; let totalTraces = 0; let totalResistance = 0; for (let li = 0; li < layersIn.length; li++) { const lyr = layersIn[li] as PcbLayer; - for (let i = 0; i + 1 < pts.length; i++) { + const lp = layerPts[li]; + for (let i = 0; i + 1 < lp.length; i++) { pcb.traces.push({ - start: pts[i], - end: pts[i + 1], + start: lp[i], + end: lp[i + 1], width: traceWidth, layer: lyr, net, @@ -7001,14 +7076,15 @@ export async function addCoil( }); totalTraces++; } - totalLengthMm += perLayerLen; + const layerLen = segLength(lp); + totalLengthMm += layerLen; const cuT = pcb.stackup.layers.find((s) => s.layer === lyr)?.copperThickness ?? 0.035; - totalResistance += (1.68e-5 * perLayerLen) / (traceWidth * cuT); - // Stitch to the next layer at an alternating terminal. + totalResistance += (1.68e-5 * layerLen) / (traceWidth * cuT); + // Stitch to the next layer at this layer's exit terminal, which is exactly + // where the next layer begins. if (li + 1 < layersIn.length) { - const atInner = li % 2 === 0; - const stitchPt = atInner ? pts[0] : pts[pts.length - 1]; + const stitchPt = lp[lp.length - 1]; pcb.vias.push({ position: stitchPt, diameter: pcb.rules.defaultRules.viaDiameter, @@ -7021,14 +7097,12 @@ export async function addCoil( stitchVias.push({ position: stitchPt, startLayer: lyr, endLayer: layersIn[li + 1] as PcbLayer }); } } - // External terminals: layer[0] outer, and the last layer's free terminal - // (inner if an even number of stitches landed on the inner side). The last - // stitch (index n-2) is at inner when (n-2)%2===0 → last free end is outer; - // otherwise inner. Equivalently the last layer's free end is inner when - // (layersIn.length-1) is even. - const lastFreeInner = (layersIn.length - 1) % 2 === 0; - const terminalA = pts[pts.length - 1]; // layer[0] outer - const terminalB = lastFreeInner ? pts[0] : pts[pts.length - 1]; + // External terminals: layer[0] starts at the outer radius; the last layer's + // free end is inner when the layer count is odd, outer when it is even. + const lastLp = layerPts[layerPts.length - 1]; + const terminalA = layerPts[0][0]; // layer[0] outer + const terminalB = lastLp[lastLp.length - 1]; + const lastFreeInner = layersIn.length % 2 === 1; return { content: [ { @@ -7041,15 +7115,17 @@ export async function addCoil( layers_used: layersIn, multilayer: true, note: - "Multilayer stacked coil: `layer` and inner_via/via_to_layer ignored; " + - "layers stitched with alternating inner/outer vias.", + "Multilayer stacked coil: `layer` and inner_via/via_to_layer ignored. " + + "One continuous spiral in angle across the stack — the radius zigzags " + + "inner↔outer with a stitch via at each turnaround — so every layer " + + "circulates the same way and their fields add.", total_traces: totalTraces, total_length_mm: Math.round(totalLengthMm * 100) / 100, total_resistance_ohms: Math.round(totalResistance * 1000) / 1000, stitch_vias: stitchVias, terminals: { a: terminalA, b: terminalB }, - inner_endpoint: pts[0], - outer_endpoint: pts[pts.length - 1], + inner_endpoint: lastFreeInner ? terminalB : layerPts[0][layerPts[0].length - 1], + outer_endpoint: terminalA, ...(drcCap ? { drc_delta: await drcCap.finish() } : {}), ...docResultPayload(ctx), }),