Add base_state from QEDbase (#28)

Copied the concrete implementations of `base_state` from `QEDbase.jl` (see https://github.com/QEDjl-project/QEDbase.jl/blob/dev/src/particles/particle_states.jl) This is part of the general restructuring of `QED.jl`, see QEDjl-project/QuantumElectrodynamics.jl#35 for details.
QEDjl-project · Jun 19, 2024 · 002d5fe · 002d5fe
2 parents a40938b + dc75423
commit 002d5fe
Show file tree

Hide file tree

Showing 4 changed files with 176 additions and 45 deletions.
diff --git a/src/QEDcore.jl b/src/QEDcore.jl
@@ -2,6 +2,7 @@ module QEDcore
 
 import Base: *
 import StaticArrays: similar_type
+import QEDbase: base_state
 
 # lorentz vectors
 export SLorentzVector, MLorentzVector
@@ -28,6 +29,9 @@ export IncomingFermionSpinor,
 export SpinorU, SpinorUbar, SpinorV, SpinorVbar
 export @valid_spinor_input
 
+# particle base states
+export base_state
+
 using QEDbase: QEDbase
 using DocStringExtensions
 using StaticArrays
@@ -42,5 +46,6 @@ include("algebraic_objects/gamma_matrices.jl")
 
 include("particles/particle_types.jl")
 include("particles/spinors.jl")
+include("particles/states.jl")
 
 end
diff --git a/src/particles/states.jl b/src/particles/states.jl
@@ -0,0 +1,136 @@
+function _booster_fermion(mom::QEDbase.AbstractFourMomentum, mass::Real)
+    return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(QEDbase.getT(mom)) + mass))
+end
+
+function _booster_antifermion(mom::QEDbase.AbstractFourMomentum, mass::Real)
+    return (mass * one(DiracMatrix) - slashed(mom)) / (sqrt(abs(QEDbase.getT(mom)) + mass))
+end
+
+function base_state(
+    particle::Fermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return BiSpinor(booster[:, QEDbase._spin_index(spin)])
+end
+
+function base_state(
+    particle::Fermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return SVector(BiSpinor(booster[:, 1]), BiSpinor(booster[:, 2]))
+end
+
+function base_state(
+    particle::AntiFermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return AdjointBiSpinor(BiSpinor(booster[:, QEDbase._spin_index(spin) + 2])) * GAMMA[1]
+end
+
+function base_state(
+    particle::AntiFermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return SVector(
+        AdjointBiSpinor(BiSpinor(booster[:, 3])) * GAMMA[1],
+        AdjointBiSpinor(BiSpinor(booster[:, 4])) * GAMMA[1],
+    )
+end
+
+function base_state(
+    particle::Fermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return AdjointBiSpinor(BiSpinor(booster[:, QEDbase._spin_index(spin)])) * GAMMA[1]
+end
+
+function base_state(
+    particle::Fermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return SVector(
+        AdjointBiSpinor(BiSpinor(booster[:, 1])) * GAMMA[1],
+        AdjointBiSpinor(BiSpinor(booster[:, 2])) * GAMMA[1],
+    )
+end
+
+function base_state(
+    particle::AntiFermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return BiSpinor(booster[:, QEDbase._spin_index(spin) + 2])
+end
+
+function base_state(
+    particle::AntiFermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return SVector(BiSpinor(booster[:, 3]), BiSpinor(booster[:, 4]))
+end
+
+function _photon_state(pol::QEDbase.AllPolarization, mom::QEDbase.AbstractFourMomentum)
+    cth = QEDbase.getCosTheta(mom)
+    sth = sqrt(1 - cth^2)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.getSinPhi(mom)
+    return SVector(
+        SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth),
+        SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0),
+    )
+end
+
+function _photon_state(pol::QEDbase.PolarizationX, mom::QEDbase.AbstractFourMomentum)
+    cth = QEDbase.getCosTheta(mom)
+    sth = sqrt(1 - cth^2)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.getSinPhi(mom)
+    return SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth)
+end
+
+function _photon_state(pol::QEDbase.PolarizationY, mom::QEDbase.AbstractFourMomentum)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.getSinPhi(mom)
+    return SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0)
+end
+
+@inline function base_state(
+    particle::Photon,
+    ::QEDbase.ParticleDirection,
+    mom::QEDbase.AbstractFourMomentum,
+    pol::QEDbase.AllPolarization,
+)
+    return _photon_state(pol, mom)
+end
+
+@inline function base_state(
+    particle::Photon,
+    ::QEDbase.ParticleDirection,
+    mom::QEDbase.AbstractFourMomentum,
+    pol::QEDbase.AbstractPolarization,
+)
+    return _photon_state(pol, mom)
+end
diff --git a/test/particles/base_states.jl → test/particles/states.jl b/test/particles/base_states.jl → test/particles/states.jl
@@ -1,5 +1,5 @@
-using QEDcore
 using QEDbase: QEDbase
+using QEDcore
 using StaticArrays
 using Random
 
@@ -52,22 +52,19 @@ test_broadcast(x::QEDbase.AbstractSpinOrPolarization) = x
         particle_mass = QEDbase.mass(p())
         groundtruth_states = FERMION_STATES_GROUNDTRUTH_FACTORY[(d, p)](mom, particle_mass)
         groundtruth_tuple = SVector(groundtruth_states(1), groundtruth_states(2))
-        @test QEDbase.base_state(p(), d(), mom, QEDbase.AllSpin()) == groundtruth_tuple
-        @test QEDbase.base_state(p(), d(), mom, QEDbase.SpinUp()) == groundtruth_tuple[1]
-        @test QEDbase.base_state(p(), d(), mom, QEDbase.SpinDown()) == groundtruth_tuple[2]
-
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.AllSpin())) isa
-            SVector
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.SpinUp())) isa
-            SVector
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.SpinDown())) isa
-            SVector
-
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.AllSpin())) ==
+        @test base_state(p(), d(), mom, QEDbase.AllSpin()) == groundtruth_tuple
+        @test base_state(p(), d(), mom, QEDbase.SpinUp()) == groundtruth_tuple[1]
+        @test base_state(p(), d(), mom, QEDbase.SpinDown()) == groundtruth_tuple[2]
+
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.AllSpin())) isa SVector
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.SpinUp())) isa SVector
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.SpinDown())) isa SVector
+
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.AllSpin())) ==
             groundtruth_tuple
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.SpinUp()))[1] ==
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.SpinUp()))[1] ==
             groundtruth_tuple[1]
-        @test QEDbase._as_svec(QEDbase.base_state(p(), d(), mom, QEDbase.SpinDown()))[1] ==
+        @test QEDbase._as_svec(base_state(p(), d(), mom, QEDbase.SpinDown()))[1] ==
             groundtruth_tuple[2]
     end
 end
@@ -87,9 +84,7 @@ end
             #@testset "$x $y $z" for (x,y,z) in Iterators.product(X_arr,Y_arr,Z_arr)
 
             mom = SFourMomentum(_cartesian_coordinates(om, om, cth, phi))
-            both_photon_states = QEDbase.base_state(
-                Photon(), D(), mom, QEDbase.AllPolarization()
-            )
+            both_photon_states = base_state(Photon(), D(), mom, QEDbase.AllPolarization())
 
             # property test the photon states
             @test isapprox((both_photon_states[1] * mom), 0.0, atol=ATOL, rtol=RTOL)
@@ -105,37 +100,28 @@ end
             )
 
             # test the single polarization states
-            @test QEDbase.base_state(Photon(), D(), mom, QEDbase.PolarizationX()) ==
+            @test base_state(Photon(), D(), mom, QEDbase.PolarizationX()) ==
                 both_photon_states[1]
-            @test QEDbase.base_state(Photon(), D(), mom, QEDbase.PolarizationY()) ==
+            @test base_state(Photon(), D(), mom, QEDbase.PolarizationY()) ==
                 both_photon_states[2]
-            @test QEDbase.base_state(Photon(), D(), mom, QEDbase.PolX()) ==
+            @test base_state(Photon(), D(), mom, QEDbase.PolX()) == both_photon_states[1]
+            @test base_state(Photon(), D(), mom, QEDbase.PolY()) == both_photon_states[2]
+
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.PolX())) isa
+                SVector
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.PolY())) isa
+                SVector
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.AllPol())) isa
+                SVector
+
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.PolX()))[1] ==
                 both_photon_states[1]
-            @test QEDbase.base_state(Photon(), D(), mom, QEDbase.PolY()) ==
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.PolY()))[1] ==
+                both_photon_states[2]
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.AllPol()))[1] ==
+                both_photon_states[1]
+            @test QEDbase._as_svec(base_state(Photon(), D(), mom, QEDbase.AllPol()))[2] ==
                 both_photon_states[2]
-
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.PolX())
-            ) isa SVector
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.PolY())
-            ) isa SVector
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.AllPol())
-            ) isa SVector
-
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.PolX())
-            )[1] == both_photon_states[1]
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.PolY())
-            )[1] == both_photon_states[2]
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.AllPol())
-            )[1] == both_photon_states[1]
-            @test QEDbase._as_svec(
-                QEDbase.base_state(Photon(), D(), mom, QEDbase.AllPol())
-            )[2] == both_photon_states[2]
         end
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -25,6 +25,10 @@ begin
         include("particles/types.jl")
     end
 
+    @time @safetestset "particle types" begin
+        include("particles/states.jl")
+    end
+
     @time @safetestset "particle spinors" begin
         include("particles/spinors.jl")
     end