proof_engine/electromagnetic/
charged_particles.rs

1//! Charged particle dynamics — Lorentz force, Boris pusher, cyclotron motion,
2//! drift velocities, magnetic mirrors, and particle tracing.
3
4use glam::{Vec3, Vec4};
5use std::f32::consts::PI;
6
7// ── Charged Particle ──────────────────────────────────────────────────────
8
9#[derive(Clone, Debug)]
10pub struct ChargedParticle {
11    pub position: Vec3,
12    pub velocity: Vec3,
13    pub mass: f32,
14    pub charge: f32,
15}
16
17impl ChargedParticle {
18    pub fn new(position: Vec3, velocity: Vec3, mass: f32, charge: f32) -> Self {
19        Self { position, velocity, mass, charge }
20    }
21
22    /// Kinetic energy: 0.5 * m * v^2
23    pub fn kinetic_energy(&self) -> f32 {
24        0.5 * self.mass * self.velocity.length_squared()
25    }
26
27    /// Momentum: m * v
28    pub fn momentum(&self) -> Vec3 {
29        self.mass * self.velocity
30    }
31}
32
33// ── Lorentz Force ─────────────────────────────────────────────────────────
34
35/// Lorentz force: F = q(E + v × B)
36pub fn lorentz_force(charge: f32, velocity: Vec3, e: Vec3, b: Vec3) -> Vec3 {
37    charge * (e + velocity.cross(b))
38}
39
40// ── Boris Pusher ──────────────────────────────────────────────────────────
41
42/// Boris algorithm for advancing a charged particle in EM fields.
43/// This is the standard velocity-Verlet-like scheme that exactly conserves
44/// energy in a pure magnetic field.
45pub fn boris_push(particle: &mut ChargedParticle, e: Vec3, b: Vec3, dt: f32) {
46    let q_over_m = particle.charge / particle.mass;
47    let half_dt = dt * 0.5;
48
49    // Half acceleration from E field
50    let v_minus = particle.velocity + e * q_over_m * half_dt;
51
52    // Rotation from B field (Boris rotation)
53    let t_vec = b * q_over_m * half_dt;
54    let t_mag_sq = t_vec.length_squared();
55    let s_vec = 2.0 * t_vec / (1.0 + t_mag_sq);
56
57    let v_prime = v_minus + v_minus.cross(t_vec);
58    let v_plus = v_minus + v_prime.cross(s_vec);
59
60    // Second half acceleration from E field
61    particle.velocity = v_plus + e * q_over_m * half_dt;
62
63    // Update position
64    particle.position += particle.velocity * dt;
65}
66
67/// Advance a particle one step in given E and B fields.
68pub fn step_particle(particle: &mut ChargedParticle, e_field: Vec3, b_field: Vec3, dt: f32) {
69    boris_push(particle, e_field, b_field, dt);
70}
71
72// ── Cyclotron Motion ──────────────────────────────────────────────────────
73
74/// Cyclotron (Larmor) radius: r_L = m*v_perp / (|q|*B)
75pub fn cyclotron_radius(particle: &ChargedParticle, b_magnitude: f32) -> f32 {
76    if b_magnitude < 1e-10 || particle.charge.abs() < 1e-10 {
77        return f32::INFINITY;
78    }
79    let v_perp = particle.velocity.length(); // Assume all velocity is perpendicular for simplicity
80    particle.mass * v_perp / (particle.charge.abs() * b_magnitude)
81}
82
83/// Cyclotron frequency: omega_c = |q|*B / m
84pub fn cyclotron_frequency(particle: &ChargedParticle, b_magnitude: f32) -> f32 {
85    if particle.mass < 1e-10 {
86        return 0.0;
87    }
88    particle.charge.abs() * b_magnitude / particle.mass
89}
90
91// ── Drift Velocities ─────────────────────────────────────────────────────
92
93/// E×B drift velocity: v_drift = (E × B) / B^2
94/// Independent of charge and mass.
95#[derive(Clone, Debug)]
96pub struct ExBDrift {
97    pub e: Vec3,
98    pub b: Vec3,
99}
100
101impl ExBDrift {
102    pub fn new(e: Vec3, b: Vec3) -> Self {
103        Self { e, b }
104    }
105
106    pub fn drift_velocity(&self) -> Vec3 {
107        let b2 = self.b.length_squared();
108        if b2 < 1e-10 {
109            return Vec3::ZERO;
110        }
111        self.e.cross(self.b) / b2
112    }
113}
114
115/// Gradient-B drift: v_grad = (m * v_perp^2) / (2*q*B^3) * (B × ∇B)
116#[derive(Clone, Debug)]
117pub struct GradBDrift;
118
119impl GradBDrift {
120    /// Compute gradient-B drift velocity.
121    /// `grad_b` is the gradient of |B| at the particle location.
122    pub fn drift_velocity(
123        particle: &ChargedParticle,
124        b: Vec3,
125        grad_b: Vec3,
126        v_perp: f32,
127    ) -> Vec3 {
128        let b_mag = b.length();
129        if b_mag < 1e-10 || particle.charge.abs() < 1e-10 {
130            return Vec3::ZERO;
131        }
132        let b_hat = b / b_mag;
133        let coeff = particle.mass * v_perp * v_perp / (2.0 * particle.charge * b_mag * b_mag);
134        coeff * b_hat.cross(grad_b)
135    }
136}
137
138/// Curvature drift: v_curv = (m * v_parallel^2) / (q * R_c^2 * B) * (R_c × B)
139#[derive(Clone, Debug)]
140pub struct CurvatureDrift;
141
142impl CurvatureDrift {
143    /// Compute curvature drift velocity.
144    /// `r_c` is the radius of curvature vector (pointing toward center of curvature).
145    pub fn drift_velocity(
146        particle: &ChargedParticle,
147        b: Vec3,
148        r_c: Vec3,
149        v_parallel: f32,
150    ) -> Vec3 {
151        let b_mag = b.length();
152        let rc_mag = r_c.length();
153        if b_mag < 1e-10 || rc_mag < 1e-10 || particle.charge.abs() < 1e-10 {
154            return Vec3::ZERO;
155        }
156        let coeff = particle.mass * v_parallel * v_parallel / (particle.charge * rc_mag * rc_mag * b_mag);
157        coeff * r_c.cross(b)
158    }
159}
160
161// ── Magnetic Mirror ───────────────────────────────────────────────────────
162
163/// Magnetic mirror: confinement between two high-field regions.
164#[derive(Clone, Debug)]
165pub struct MagneticMirror {
166    pub b_min: f32,  // minimum B (at center)
167    pub b_max: f32,  // maximum B (at mirror points)
168    pub length: f32,  // distance between mirrors
169}
170
171impl MagneticMirror {
172    pub fn new(b_min: f32, b_max: f32, length: f32) -> Self {
173        Self { b_min, b_max, length }
174    }
175
176    /// Mirror ratio: R = B_max / B_min
177    pub fn mirror_ratio(&self) -> f32 {
178        if self.b_min < 1e-10 {
179            return f32::INFINITY;
180        }
181        self.b_max / self.b_min
182    }
183
184    /// Loss cone angle: sin^2(alpha) = B_min / B_max = 1/R
185    pub fn loss_cone_angle(&self) -> f32 {
186        let r = self.mirror_ratio();
187        if r < 1.0 {
188            return PI / 2.0;
189        }
190        (1.0 / r).sqrt().asin()
191    }
192
193    /// Check if a particle with given pitch angle is confined.
194    /// pitch_angle is the angle between velocity and B field.
195    pub fn is_confined(&self, pitch_angle: f32) -> bool {
196        pitch_angle > self.loss_cone_angle()
197    }
198
199    /// Magnetic field magnitude as a function of position along the axis.
200    /// Simple model: B(z) = B_min * (1 + (R-1) * (2*z/L)^2) for z in [-L/2, L/2]
201    pub fn field_at_position(&self, z: f32) -> f32 {
202        let normalized_z = 2.0 * z / self.length;
203        let r = self.mirror_ratio();
204        self.b_min * (1.0 + (r - 1.0) * normalized_z * normalized_z)
205    }
206
207    /// Bounce period for a trapped particle.
208    /// Approximate: T ~ 2*L / v_parallel
209    pub fn bounce_period(&self, v_parallel: f32) -> f32 {
210        if v_parallel.abs() < 1e-10 {
211            return f32::INFINITY;
212        }
213        2.0 * self.length / v_parallel.abs()
214    }
215}
216
217// ── Particle Tracer ───────────────────────────────────────────────────────
218
219/// Traces particle trajectories and renders as glyph trails.
220pub struct ParticleTracer {
221    pub trail_length: usize,
222    pub color_by_velocity: bool,
223    pub min_speed_color: Vec4, // color at low speed
224    pub max_speed_color: Vec4, // color at high speed
225    pub max_speed: f32,
226}
227
228impl ParticleTracer {
229    pub fn new(trail_length: usize) -> Self {
230        Self {
231            trail_length,
232            color_by_velocity: true,
233            min_speed_color: Vec4::new(0.2, 0.3, 1.0, 1.0),
234            max_speed_color: Vec4::new(1.0, 0.2, 0.1, 1.0),
235            max_speed: 10.0,
236        }
237    }
238
239    /// Trace a particle for `steps` timesteps, recording positions.
240    pub fn trace(
241        &self,
242        particle: &mut ChargedParticle,
243        e_field: Vec3,
244        b_field: Vec3,
245        dt: f32,
246        steps: usize,
247    ) -> Vec<(Vec3, Vec4)> {
248        let mut trail = Vec::with_capacity(steps);
249        for _ in 0..steps {
250            let speed = particle.velocity.length();
251            let t = (speed / self.max_speed).clamp(0.0, 1.0);
252            let color = self.min_speed_color * (1.0 - t) + self.max_speed_color * t;
253            trail.push((particle.position, color));
254            boris_push(particle, e_field, b_field, dt);
255        }
256        // Keep only the most recent trail_length points
257        if trail.len() > self.trail_length {
258            trail.drain(0..trail.len() - self.trail_length);
259        }
260        trail
261    }
262
263    /// Get a glyph character for the particle based on charge sign.
264    pub fn particle_glyph(charge: f32) -> char {
265        if charge > 0.0 { '+' }
266        else if charge < 0.0 { '-' }
267        else { '○' }
268    }
269
270    /// Trail glyph based on direction of motion.
271    pub fn trail_glyph(direction: Vec3) -> char {
272        let angle = direction.y.atan2(direction.x);
273        let octant = ((angle / (PI / 4.0)).round() as i32).rem_euclid(8);
274        match octant {
275            0 => '→',
276            1 => '↗',
277            2 => '↑',
278            3 => '↖',
279            4 => '←',
280            5 => '↙',
281            6 => '↓',
282            7 => '↘',
283            _ => '·',
284        }
285    }
286}
287
288// ── Charged Particle System ───────────────────────────────────────────────
289
290/// Manages a collection of charged particles with field evaluation.
291pub struct ChargedParticleSystem {
292    pub particles: Vec<ChargedParticle>,
293    pub uniform_e: Vec3,
294    pub uniform_b: Vec3,
295}
296
297impl ChargedParticleSystem {
298    pub fn new() -> Self {
299        Self {
300            particles: Vec::new(),
301            uniform_e: Vec3::ZERO,
302            uniform_b: Vec3::ZERO,
303        }
304    }
305
306    pub fn add_particle(&mut self, particle: ChargedParticle) {
307        self.particles.push(particle);
308    }
309
310    /// Advance all particles one timestep.
311    pub fn step(&mut self, dt: f32) {
312        let e = self.uniform_e;
313        let b = self.uniform_b;
314        for p in &mut self.particles {
315            boris_push(p, e, b, dt);
316        }
317    }
318
319    /// Total kinetic energy of all particles.
320    pub fn total_kinetic_energy(&self) -> f32 {
321        self.particles.iter().map(|p| p.kinetic_energy()).sum()
322    }
323
324    /// Total momentum.
325    pub fn total_momentum(&self) -> Vec3 {
326        self.particles.iter().map(|p| p.momentum()).sum()
327    }
328
329    /// Evaluate the E field at a point including contributions from all particles
330    /// (Coulomb fields from each particle, treated as point charges).
331    pub fn e_field_at(&self, pos: Vec3) -> Vec3 {
332        let mut field = self.uniform_e;
333        for p in &self.particles {
334            let r_vec = pos - p.position;
335            let r2 = r_vec.length_squared();
336            if r2 < 1e-10 {
337                continue;
338            }
339            let r = r2.sqrt();
340            field += p.charge / r2 * (r_vec / r);
341        }
342        field
343    }
344
345    /// Step with self-consistent fields (N-body Coulomb + uniform B).
346    /// O(N^2) — fine for small N.
347    pub fn step_self_consistent(&mut self, dt: f32) {
348        let n = self.particles.len();
349        let mut forces = vec![Vec3::ZERO; n];
350
351        // Compute Coulomb forces between all pairs
352        for i in 0..n {
353            for j in (i + 1)..n {
354                let r_vec = self.particles[j].position - self.particles[i].position;
355                let r2 = r_vec.length_squared();
356                if r2 < 1e-6 {
357                    continue;
358                }
359                let r = r2.sqrt();
360                let f_mag = self.particles[i].charge * self.particles[j].charge / r2;
361                let f = f_mag * (r_vec / r);
362                forces[i] -= f; // repulsive for same sign
363                forces[j] += f;
364            }
365        }
366
367        // Apply forces + uniform fields using Boris push
368        let b = self.uniform_b;
369        for i in 0..n {
370            let e_total = self.uniform_e + forces[i] / self.particles[i].charge.abs().max(1e-10);
371            boris_push(&mut self.particles[i], e_total, b, dt);
372        }
373    }
374}
375
376impl Default for ChargedParticleSystem {
377    fn default() -> Self {
378        Self::new()
379    }
380}
381
382// ── Tests ─────────────────────────────────────────────────────────────────
383
384#[cfg(test)]
385mod tests {
386    use super::*;
387
388    #[test]
389    fn test_lorentz_force_electric_only() {
390        let f = lorentz_force(1.0, Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), Vec3::ZERO);
391        assert!((f - Vec3::new(1.0, 0.0, 0.0)).length() < 1e-6);
392    }
393
394    #[test]
395    fn test_lorentz_force_magnetic_only() {
396        // v along x, B along z => F along y (positive charge)
397        let f = lorentz_force(1.0, Vec3::X, Vec3::ZERO, Vec3::Z);
398        // v × B = (1,0,0) × (0,0,1) = (0*1 - 0*0, 0*0 - 1*1, 1*0 - 0*0) = (0,-1,0)
399        assert!((f.y - (-1.0)).abs() < 1e-6, "f={:?}", f);
400    }
401
402    #[test]
403    fn test_cyclotron_radius() {
404        let p = ChargedParticle::new(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), 1.0, 1.0);
405        let r = cyclotron_radius(&p, 1.0);
406        // r_L = m*v / (q*B) = 1*1 / (1*1) = 1
407        assert!((r - 1.0).abs() < 1e-6);
408    }
409
410    #[test]
411    fn test_cyclotron_frequency() {
412        let p = ChargedParticle::new(Vec3::ZERO, Vec3::ZERO, 2.0, 3.0);
413        let f = cyclotron_frequency(&p, 4.0);
414        // omega_c = |q|*B / m = 3*4/2 = 6
415        assert!((f - 6.0).abs() < 1e-6);
416    }
417
418    #[test]
419    fn test_boris_energy_conservation_pure_b() {
420        // In a pure magnetic field, the Boris pusher should conserve kinetic energy exactly
421        let mut p = ChargedParticle::new(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), 1.0, 1.0);
422        let b = Vec3::new(0.0, 0.0, 1.0);
423        let e = Vec3::ZERO;
424        let dt = 0.01;
425
426        let e0 = p.kinetic_energy();
427        for _ in 0..1000 {
428            boris_push(&mut p, e, b, dt);
429        }
430        let e1 = p.kinetic_energy();
431        assert!((e0 - e1).abs() < 1e-5, "Energy should be conserved: e0={}, e1={}", e0, e1);
432    }
433
434    #[test]
435    fn test_boris_circular_orbit() {
436        // A charged particle in uniform B should trace a circle
437        let mut p = ChargedParticle::new(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), 1.0, 1.0);
438        let b = Vec3::new(0.0, 0.0, 1.0);
439        let e = Vec3::ZERO;
440        let dt = 0.01;
441        let omega = cyclotron_frequency(&p, 1.0); // = 1
442        let period = 2.0 * PI / omega;
443        let steps = (period / dt) as usize;
444
445        for _ in 0..steps {
446            boris_push(&mut p, e, b, dt);
447        }
448
449        // After one full period, should return close to start
450        assert!(p.position.length() < 0.5, "Should return near origin after one period: {:?}", p.position);
451    }
452
453    #[test]
454    fn test_exb_drift() {
455        let drift = ExBDrift::new(Vec3::new(1.0, 0.0, 0.0), Vec3::new(0.0, 0.0, 1.0));
456        let v = drift.drift_velocity();
457        // E×B / B^2 = (1,0,0)×(0,0,1) / 1 = (0,-1,0)
458        // Actually: (1,0,0)×(0,0,1) = (0*1-0*0, 0*0-1*1, 1*0-0*0) = (0,-1,0)
459        assert!((v.y - (-1.0)).abs() < 1e-6, "drift={:?}", v);
460    }
461
462    #[test]
463    fn test_magnetic_mirror() {
464        let mirror = MagneticMirror::new(1.0, 4.0, 10.0);
465        assert!((mirror.mirror_ratio() - 4.0).abs() < 1e-6);
466
467        let loss_cone = mirror.loss_cone_angle();
468        // sin^2(alpha) = 1/R = 0.25, sin(alpha) = 0.5, alpha = pi/6
469        assert!((loss_cone - PI / 6.0).abs() < 0.01, "loss_cone={}", loss_cone);
470
471        // A particle with pitch angle > loss_cone should be confined
472        assert!(mirror.is_confined(PI / 4.0)); // 45° > 30°
473        assert!(!mirror.is_confined(PI / 12.0)); // 15° < 30°
474    }
475
476    #[test]
477    fn test_mirror_field_profile() {
478        let mirror = MagneticMirror::new(1.0, 4.0, 10.0);
479        // At center (z=0), B = B_min
480        assert!((mirror.field_at_position(0.0) - 1.0).abs() < 1e-6);
481        // At ends (z=±L/2), B = B_max
482        assert!((mirror.field_at_position(5.0) - 4.0).abs() < 1e-6);
483        assert!((mirror.field_at_position(-5.0) - 4.0).abs() < 1e-6);
484    }
485
486    #[test]
487    fn test_particle_system() {
488        let mut sys = ChargedParticleSystem::new();
489        sys.uniform_b = Vec3::new(0.0, 0.0, 1.0);
490        sys.add_particle(ChargedParticle::new(
491            Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), 1.0, 1.0,
492        ));
493        let e0 = sys.total_kinetic_energy();
494        for _ in 0..100 {
495            sys.step(0.01);
496        }
497        let e1 = sys.total_kinetic_energy();
498        assert!((e0 - e1).abs() < 1e-4, "System energy should be conserved in pure B");
499    }
500
501    #[test]
502    fn test_particle_tracer() {
503        let tracer = ParticleTracer::new(50);
504        let mut p = ChargedParticle::new(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), 1.0, 1.0);
505        let trail = tracer.trace(&mut p, Vec3::ZERO, Vec3::Z, 0.01, 100);
506        assert!(trail.len() <= 50, "Trail should be limited to trail_length");
507        assert!(trail.len() == 50);
508    }
509
510    #[test]
511    fn test_grad_b_drift() {
512        let p = ChargedParticle::new(Vec3::ZERO, Vec3::ZERO, 1.0, 1.0);
513        let b = Vec3::new(0.0, 0.0, 1.0);
514        let grad_b = Vec3::new(0.1, 0.0, 0.0); // B increases in x
515        let v_drift = GradBDrift::drift_velocity(&p, b, grad_b, 1.0);
516        // Should drift in y direction (B × ∇B for positive charge)
517        assert!(v_drift.y.abs() > 1e-6, "Should have y-drift: {:?}", v_drift);
518    }
519}
proof_engine/electromagnetic/charged_particles.rs

proof_engine/electromagnetic/
charged_particles.rs
