versorc_math/sta/
mod.rs

1// ============================================================================
2// VERSORC-MATH: SPACETIME ALGEBRA (STA) Cl(1, 3) (MINKOWSKI SPACETIME)
3// License: Apache-2.0
4// ============================================================================
5
6use std::ops::{Add, Div, Mul, Neg, Sub};
7
8/// Minkowski 4-vector in Spacetime Algebra Cl(1, 3).
9/// Basis: \gamma_0^2 = +1 (timelike), \gamma_1^2 = \gamma_2^2 = \gamma_3^2 = -1 (spacelike).
10/// Metric interval: s^2 = (ct)^2 - (x^2 + y^2 + z^2).
11#[repr(C, align(16))]
12#[derive(Debug, Clone, Copy, PartialEq)]
13#[cfg_attr(feature = "vulkan", derive(bytemuck::Pod, bytemuck::Zeroable))]
14pub struct StaVector {
15    pub ct: f32, // \gamma_0: Time component (c * t)
16    pub x: f32,  // \gamma_1: Spatial x
17    pub y: f32,  // \gamma_2: Spatial y
18    pub z: f32,  // \gamma_3: Spatial z
19}
20
21impl Default for StaVector {
22    #[inline(always)]
23    fn default() -> Self {
24        Self::zero()
25    }
26}
27
28impl StaVector {
29    #[inline(always)]
30    pub const fn new(ct: f32, x: f32, y: f32, z: f32) -> Self {
31        Self { ct, x, y, z }
32    }
33
34    #[inline(always)]
35    pub const fn zero() -> Self {
36        Self { ct: 0.0, x: 0.0, y: 0.0, z: 0.0 }
37    }
38
39    /// Spacetime interval invariant s^2 = (ct)^2 - (x^2 + y^2 + z^2).
40    #[inline(always)]
41    pub fn interval_sq(&self) -> f32 {
42        self.ct * self.ct - (self.x * self.x + self.y * self.y + self.z * self.z)
43    }
44
45    /// Minkowski inner product: u . v = u_0 * v_0 - (u_x*v_x + u_y*v_y + u_z*v_z).
46    #[inline(always)]
47    pub fn dot(&self, other: &Self) -> f32 {
48        self.ct * other.ct - (self.x * other.x + self.y * other.y + self.z * other.z)
49    }
50
51    /// Checks if interval is timelike (s^2 > 0).
52    #[inline(always)]
53    pub fn is_timelike(&self) -> bool {
54        self.interval_sq() > 0.0
55    }
56
57    /// Checks if interval is spacelike (s^2 < 0).
58    #[inline(always)]
59    pub fn is_spacelike(&self) -> bool {
60        self.interval_sq() < 0.0
61    }
62
63    /// Checks if interval is lightlike / null (|s^2| < 1e-12).
64    #[inline(always)]
65    pub fn is_lightlike(&self) -> bool {
66        self.interval_sq().abs() < 1e-12
67    }
68
69    /// Proper time interval \Delta \tau = sqrt(s^2) / c (with c = 1).
70    #[inline]
71    pub fn proper_time(&self) -> Option<f32> {
72        let s2 = self.interval_sq();
73        if s2 >= 0.0 {
74            Some(s2.sqrt())
75        } else {
76            None
77        }
78    }
79}
80
81impl Add for StaVector {
82    type Output = Self;
83    #[inline(always)]
84    fn add(self, rhs: Self) -> Self {
85        Self {
86            ct: self.ct + rhs.ct,
87            x: self.x + rhs.x,
88            y: self.y + rhs.y,
89            z: self.z + rhs.z,
90        }
91    }
92}
93
94impl Sub for StaVector {
95    type Output = Self;
96    #[inline(always)]
97    fn sub(self, rhs: Self) -> Self {
98        Self {
99            ct: self.ct - rhs.ct,
100            x: self.x - rhs.x,
101            y: self.y - rhs.y,
102            z: self.z - rhs.z,
103        }
104    }
105}
106
107impl Neg for StaVector {
108    type Output = Self;
109    #[inline(always)]
110    fn neg(self) -> Self {
111        Self {
112            ct: -self.ct,
113            x: -self.x,
114            y: -self.y,
115            z: -self.z,
116        }
117    }
118}
119
120impl Mul<f32> for StaVector {
121    type Output = Self;
122    #[inline(always)]
123    fn mul(self, rhs: f32) -> Self {
124        Self {
125            ct: self.ct * rhs,
126            x: self.x * rhs,
127            y: self.y * rhs,
128            z: self.z * rhs,
129        }
130    }
131}
132
133impl Div<f32> for StaVector {
134    type Output = Self;
135    #[inline(always)]
136    fn div(self, rhs: f32) -> Self {
137        let inv = 1.0 / rhs;
138        Self {
139            ct: self.ct * inv,
140            x: self.x * inv,
141            y: self.y * inv,
142            z: self.z * inv,
143        }
144    }
145}
146
147/// Electromagnetic Faraday Bivector in Spacetime Algebra: F = E + I*c*B.
148/// Unifies electric and magnetic field vectors into a single grade-2 bivector.
149#[repr(C, align(32))]
150#[derive(Debug, Clone, Copy, PartialEq)]
151#[cfg_attr(feature = "vulkan", derive(bytemuck::Pod, bytemuck::Zeroable))]
152pub struct StaBivector {
153    pub ex: f32, // \gamma_1 \gamma_0: Electric field x
154    pub ey: f32, // \gamma_2 \gamma_0: Electric field y
155    pub ez: f32, // \gamma_3 \gamma_0: Electric field z
156    pub bx: f32, // \gamma_3 \gamma_2: Magnetic field x
157    pub by: f32, // \gamma_1 \gamma_3: Magnetic field y
158    pub bz: f32, // \gamma_2 \gamma_1: Magnetic field z
159    pub _pad0: f32,
160    pub _pad1: f32,
161}
162
163impl Default for StaBivector {
164    #[inline(always)]
165    fn default() -> Self {
166        Self::zero()
167    }
168}
169
170impl StaBivector {
171    #[inline(always)]
172    pub const fn new(ex: f32, ey: f32, ez: f32, bx: f32, by: f32, bz: f32) -> Self {
173        Self {
174            ex,
175            ey,
176            ez,
177            bx,
178            by,
179            bz,
180            _pad0: 0.0,
181            _pad1: 0.0,
182        }
183    }
184
185    #[inline(always)]
186    pub const fn zero() -> Self {
187        Self::new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
188    }
189
190    /// First relativistic invariant: E^2 - c^2*B^2 (with c = 1).
191    #[inline(always)]
192    pub fn invariant_1(&self) -> f32 {
193        (self.ex * self.ex + self.ey * self.ey + self.ez * self.ez)
194            - (self.bx * self.bx + self.by * self.by + self.bz * self.bz)
195    }
196
197    /// Second relativistic invariant: 2 * E . B.
198    #[inline(always)]
199    pub fn invariant_2(&self) -> f32 {
200        2.0 * (self.ex * self.bx + self.ey * self.by + self.ez * self.bz)
201    }
202
203    /// Poynting vector direction (energy flux): S = E x B.
204    #[inline]
205    pub fn poynting_vector(&self) -> [f32; 3] {
206        [
207            self.ey * self.bz - self.ez * self.by,
208            self.ez * self.bx - self.ex * self.bz,
209            self.ex * self.by - self.ey * self.bx,
210        ]
211    }
212}
213
214/// Lorentz Boost Versor (hyperbolic rotor) along velocity vector v = (vx, vy, vz).
215/// B = cosh(phi/2) - sinh(phi/2) * (v_hat . \sigma).
216#[repr(C, align(16))]
217#[derive(Debug, Clone, Copy, PartialEq)]
218#[cfg_attr(feature = "vulkan", derive(bytemuck::Pod, bytemuck::Zeroable))]
219pub struct StaLorentzBoost {
220    pub cosh_half: f32,
221    pub sinh_half_x: f32,
222    pub sinh_half_y: f32,
223    pub sinh_half_z: f32,
224}
225
226impl StaLorentzBoost {
227    /// Constructs a Lorentz boost from velocity vector v = (vx, vy, vz) with c = 1.
228    /// beta = ||v|| < 1.0.
229    #[inline]
230    pub fn from_velocity(vx: f32, vy: f32, vz: f32) -> Self {
231        let v2 = vx * vx + vy * vy + vz * vz;
232        let beta = v2.sqrt();
233        assert!(beta < 1.0, "Velocity must be strictly subluminal (|v| < c)");
234
235        if beta < 1e-12 {
236            return Self {
237                cosh_half: 1.0,
238                sinh_half_x: 0.0,
239                sinh_half_y: 0.0,
240                sinh_half_z: 0.0,
241            };
242        }
243
244        // Rapidity phi = atanh(beta)
245        let gamma = 1.0 / (1.0 - beta * beta).sqrt();
246        // Half-rapidity identities:
247        // cosh(phi/2) = sqrt((gamma + 1) / 2)
248        // sinh(phi/2) = sqrt((gamma - 1) / 2)
249        let cosh_half = ((gamma + 1.0) * 0.5).sqrt();
250        let sinh_half = ((gamma - 1.0) * 0.5).sqrt();
251
252        let inv_beta = 1.0 / beta;
253        Self {
254            cosh_half,
255            sinh_half_x: sinh_half * (vx * inv_beta),
256            sinh_half_y: sinh_half * (vy * inv_beta),
257            sinh_half_z: sinh_half * (vz * inv_beta),
258        }
259    }
260
261    /// Transforms a Minkowski 4-vector via 2-versor sandwich x' = B * x * B~.
262    #[inline]
263    pub fn transform_vector(&self, v: &StaVector) -> StaVector {
264        let ch = self.cosh_half;
265        let sx = self.sinh_half_x;
266        let sy = self.sinh_half_y;
267        let sz = self.sinh_half_z;
268
269        let s2 = sx * sx + sy * sy + sz * sz;
270        let gamma = ch * ch + s2;
271        let two_ch = 2.0 * ch;
272
273        // Boost transformation formulas:
274        let s_dot_r = sx * v.x + sy * v.y + sz * v.z;
275        let new_ct = gamma * v.ct - two_ch * s_dot_r;
276
277        let factor = if s2 > 1e-12 {
278            (gamma - 1.0) * s_dot_r / s2
279        } else {
280            0.0
281        };
282
283        let new_x = v.x - two_ch * v.ct * sx + factor * sx;
284        let new_y = v.y - two_ch * v.ct * sy + factor * sy;
285        let new_z = v.z - two_ch * v.ct * sz + factor * sz;
286
287        StaVector::new(new_ct, new_x, new_y, new_z)
288    }
289}
290
291/// The 6 Spinor Classes according to Pertti Lounesto's Classification Theorem:
292/// Lounesto (2001), "Clifford Algebras and Spinors", Chapter 12 (Flags, Poles and Dipoles).
293#[derive(Debug, Clone, Copy, PartialEq, Eq)]
294pub enum LounestoClass {
295    /// Class 1: Regular Dirac spinor with both non-zero scalar invariants (sigma != 0, omega != 0).
296    Class1DiracGeneral,
297    /// Class 2: Regular Dirac electron spinor with omega == 0 (Yvon-Takabayasi angle beta in {0, pi}).
298    Class2DiracElectron,
299    /// Class 3: Dirac spinor with sigma == 0 and omega != 0.
300    Class3DiracPseudoscalar,
301    /// Class 4: Singular Flag-Dipole spinor (sigma == omega == 0, S != 0, K != 0).
302    Class4FlagDipole,
303    /// Class 5: Singular Flag-Pole / Majorana spinor (sigma == omega == K == 0, S != 0).
304    Class5MajoranaFlagPole,
305    /// Class 6: Singular Dipole / Weyl neutrino spinor (sigma == omega == S == 0, K != 0).
306    Class6WeylDipole,
307}
308
309/// Bilinear Covariants of a Dirac Spinor in Cl(1, 3):
310/// sigma (scalar R1), omega (pseudoscalar R2), J (current), S (spin density bivector), K (chiral current).
311#[derive(Debug, Clone, Copy, PartialEq)]
312pub struct SpinorBilinearCovariants {
313    pub sigma: f32,       // Scalar invariant <psi psi~>_0
314    pub omega: f32,       // Pseudoscalar invariant <psi psi~>_4
315    pub j: StaVector,     // Dirac probability current 4-vector
316    pub s: StaBivector,   // Spin density bivector
317    pub k: StaVector,     // Chiral / pseudovector current
318}
319
320/// Dirac Spinor in Spacetime Algebra Cl(1, 3) as an Even Multivector Cl+(1, 3):
321/// psi = scalar + (ex \gamma_10 + ey \gamma_20 + ez \gamma_30 + bx \gamma_32 + by \gamma_13 + bz \gamma_21) + pseudoscalar \gamma_0123.
322/// Alignment: 32 bytes (cache line half, 8 floats).
323#[repr(C, align(32))]
324#[derive(Debug, Clone, Copy, PartialEq)]
325#[cfg_attr(feature = "vulkan", derive(bytemuck::Pod, bytemuck::Zeroable))]
326pub struct StaSpinor {
327    pub scalar: f32,       // Grade 0: alpha
328    pub ex: f32,           // Grade 2: boost x (\gamma_1 \gamma_0)
329    pub ey: f32,           // Grade 2: boost y (\gamma_2 \gamma_0)
330    pub ez: f32,           // Grade 2: boost z (\gamma_3 \gamma_0)
331    pub bx: f32,           // Grade 2: rotation x (\gamma_3 \gamma_2)
332    pub by: f32,           // Grade 2: rotation y (\gamma_1 \gamma_3)
333    pub bz: f32,           // Grade 2: rotation z (\gamma_2 \gamma_1)
334    pub pseudoscalar: f32, // Grade 4: beta (\gamma_0123)
335}
336
337impl Default for StaSpinor {
338    #[inline(always)]
339    fn default() -> Self {
340        Self::zero()
341    }
342}
343
344impl StaSpinor {
345    #[inline(always)]
346    pub const fn zero() -> Self {
347        Self {
348            scalar: 0.0,
349            ex: 0.0,
350            ey: 0.0,
351            ez: 0.0,
352            bx: 0.0,
353            by: 0.0,
354            bz: 0.0,
355            pseudoscalar: 0.0,
356        }
357    }
358
359    #[inline(always)]
360    pub const fn new(
361        scalar: f32,
362        ex: f32,
363        ey: f32,
364        ez: f32,
365        bx: f32,
366        by: f32,
367        bz: f32,
368        pseudoscalar: f32,
369    ) -> Self {
370        Self {
371            scalar,
372            ex,
373            ey,
374            ez,
375            bx,
376            by,
377            bz,
378            pseudoscalar,
379        }
380    }
381
382    /// Reversion of spinor: psi~ = scalar - B + pseudoscalar * I.
383    #[inline(always)]
384    pub fn reverse(&self) -> Self {
385        Self {
386            scalar: self.scalar,
387            ex: -self.ex,
388            ey: -self.ey,
389            ez: -self.ez,
390            bx: -self.bx,
391            by: -self.by,
392            bz: -self.bz,
393            pseudoscalar: self.pseudoscalar,
394        }
395    }
396
397    /// Computes Lounesto's 5 Bilinear Covariants (sigma, omega, J, S, K) (Lounesto 2001, Eq. 12.1).
398    #[inline]
399    pub fn bilinear_covariants(&self) -> SpinorBilinearCovariants {
400        let s = self.scalar;
401        let p = self.pseudoscalar;
402        let ex = self.ex;
403        let ey = self.ey;
404        let ez = self.ez;
405        let bx = self.bx;
406        let by = self.by;
407        let bz = self.bz;
408
409        // Scalar invariant: sigma = <psi psi~>_0 = s^2 - p^2 - (ex^2 + ey^2 + ez^2) + (bx^2 + by^2 + bz^2)
410        let sigma = (s * s - p * p) - (ex * ex + ey * ey + ez * ez) + (bx * bx + by * by + bz * bz);
411
412        // Pseudoscalar invariant: omega = <psi psi~>_4 = 2*s*p - 2*(ex*bx + ey*by + ez*bz)
413        let omega = 2.0 * (s * p - (ex * bx + ey * by + ez * bz));
414
415        // Dirac probability current J^mu = <psi \gamma^mu \psi~>_1
416        let j0 = s * s + p * p + (ex * ex + ey * ey + ez * ez) + (bx * bx + by * by + bz * bz);
417        let jx = 2.0 * (s * ex + p * bx + ey * bz - ez * by);
418        let jy = 2.0 * (s * ey + p * by + ez * bx - ex * bz);
419        let jz = 2.0 * (s * ez + p * bz + ex * by - ey * bx);
420        let j = StaVector::new(j0, jx, jy, jz);
421
422        // Spin density bivector S = <psi \gamma_{12} \psi~>_2
423        let sex = 2.0 * (s * bx - p * ex + ey * ez);
424        let sey = 2.0 * (s * by - p * ey + ez * ex);
425        let sez = 2.0 * (s * bz - p * ez + ex * ey);
426        let sbx = s * s - p * p + ex * ex - ey * ey - ez * ez - bx * bx + by * by + bz * bz;
427        let sby = 2.0 * (ex * ey - bx * by);
428        let sbz = 2.0 * (ex * ez - bx * bz);
429        let s_biv = StaBivector::new(sex, sey, sez, sbx, sby, sbz);
430
431        // Chiral current / pseudovector K^mu
432        let k0 = 2.0 * (s * p + ex * bx + ey * by + ez * bz);
433        let kx = 2.0 * (p * ex - s * bx + ey * bz - ez * by);
434        let ky = 2.0 * (p * ey - s * by + ez * bx - ex * bz);
435        let kz = 2.0 * (p * ez - s * bz + ex * by - ey * bx);
436        let k = StaVector::new(k0, kx, ky, kz);
437
438        SpinorBilinearCovariants {
439            sigma,
440            omega,
441            j,
442            s: s_biv,
443            k,
444        }
445    }
446
447    /// Classifies the spinor into one of the 6 Lounesto classes (Lounesto 2001, Ch. 12).
448    #[inline]
449    pub fn lounesto_class(&self, tol: f32) -> LounestoClass {
450        let cov = self.bilinear_covariants();
451        let sig_nonzero = cov.sigma.abs() > tol;
452        let ome_nonzero = cov.omega.abs() > tol;
453
454        if sig_nonzero && ome_nonzero {
455            LounestoClass::Class1DiracGeneral
456        } else if sig_nonzero && !ome_nonzero {
457            LounestoClass::Class2DiracElectron
458        } else if !sig_nonzero && ome_nonzero {
459            LounestoClass::Class3DiracPseudoscalar
460        } else {
461            // Singular spinors (sigma == 0 and omega == 0)
462            let s_norm_sq = cov.s.invariant_1().abs() + cov.s.invariant_2().abs();
463            let k_norm_sq = cov.k.interval_sq().abs();
464            let s_nonzero = s_norm_sq > tol * tol;
465            let k_nonzero = k_norm_sq > tol * tol;
466
467            if s_nonzero && k_nonzero {
468                LounestoClass::Class4FlagDipole
469            } else if s_nonzero && !k_nonzero {
470                LounestoClass::Class5MajoranaFlagPole
471            } else {
472                LounestoClass::Class6WeylDipole
473            }
474        }
475    }
476}
versorc_math/sta/mod.rs

versorc_math/sta/
mod.rs
