From 0913303f6e2c6130e9b60d7a350c0245704822da Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Stefanos=20Carlstr=C3=B6m?= <stefanos.carlstrom@gmail.com>
Date: Tue, 1 Aug 2023 08:04:10 +0200
Subject: [PATCH] Implemented rotation of fields, post facto

---
 src/field_types.jl  | 13 ++++++++++-
 src/rotations.jl    |  2 ++
 test/arithmetic.jl  |  1 +
 test/field_types.jl | 56 ++++++++++++++++++++++++++++++++++++---------
 4 files changed, 60 insertions(+), 12 deletions(-)

diff --git a/src/field_types.jl b/src/field_types.jl
index abee5fe..1d8a491 100644
--- a/src/field_types.jl
+++ b/src/field_types.jl
@@ -688,6 +688,11 @@ Base.convert(::Type{<:TransverseField}, f::LinearField) =
                     envelope(f), f.I₀, f.E₀, f.A₀,
                     I, f.params)
 
+rotate(f::TransverseField, R) =
+    TransverseField(f.carrier, f.env,
+                    f.I₀, f.E₀, f.A₀,
+                    compute_rotation(R*f.R), f.params)
+
 # * Linearly polarized transverse field
 
 @doc raw"""
@@ -748,6 +753,9 @@ phase_shift(f::LinearTransverseField, δϕ) =
 
 time_integral(f::LinearTransverseField) = time_integral(envelope(f))
 
+rotate(f::LinearTransverseField, R) =
+    LinearTransverseField(f.linear_field, compute_rotation(R*f.R))
+
 # * Constant field
 
 @doc raw"""
@@ -1006,6 +1014,8 @@ transverse_field(::LinearPolarization, f) = LinearTransverseField(f)
 transverse_field(::ArbitraryPolarization, f) = f
 transverse_field(f) = transverse_field(polarization(f), f)
 
+rotate(f, R) = rotate(transverse_field(f), R)
+
 # * Exports
 
 export LinearPolarization, ArbitraryPolarization, polarization,
@@ -1020,4 +1030,5 @@ export LinearPolarization, ArbitraryPolarization, polarization,
     params, dimensions,
     phase_shift, phase,
     field_amplitude_spectrum, vector_potential_spectrum,
-    transverse_field
+    transverse_field,
+    rotate, rotation_matrix
diff --git a/src/rotations.jl b/src/rotations.jl
index 5a13813..ba1892f 100644
--- a/src/rotations.jl
+++ b/src/rotations.jl
@@ -36,3 +36,5 @@ function rotation_axis(R::AbstractMatrix)
     i = argmin(abs.(ee.values .- 1))
     normalize(real(ee.vectors[:,i]))
 end
+
+export rotation_angle, rotation_axis, rotation_matrix
diff --git a/test/arithmetic.jl b/test/arithmetic.jl
index 6d32a7f..5045be6 100644
--- a/test/arithmetic.jl
+++ b/test/arithmetic.jl
@@ -334,6 +334,7 @@ test_addition(A, B) = test_addition(polarization(A+B), A, B)
         test_addition(C, D)
         test_addition(A, CpD)
         test_addition(ApB, CpD)
+        test_addition(ApB, rotate(CpD, [1 0 0; 0 1 1; 0 -1 1]))
 
         ApBpCpD = ApB + CpD
 
diff --git a/test/field_types.jl b/test/field_types.jl
index 462bed0..0bfd8c8 100644
--- a/test/field_types.jl
+++ b/test/field_types.jl
@@ -35,7 +35,7 @@
                                  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
         end
 
-        @test ElectricFields.rotation_matrix(F) == Matrix(I, 3, 3)
+        @test rotation_matrix(F) == Matrix(I, 3, 3)
 
         t = [0.4, 0.5]
         for fun in (vector_potential, field_amplitude, intensity)
@@ -119,9 +119,9 @@
                                  – Uₚ = 0.0507 Ha = 1.3785 eV => α = 0.1433 Bohr = 7.5826 pm"""
         end
 
-        @test ElectricFields.rotation_matrix(F) ≈ [1/√2 -1/√2 0
-                                                   1/√2 1/√2  0
-                                                   0    0     1]
+        @test rotation_matrix(F) ≈ [1/√2 -1/√2 0
+                                    1/√2 1/√2  0
+                                    0    0     1]
 
         t = [0.4, 0.5]
         for fun in (vector_potential, field_amplitude, intensity)
@@ -333,7 +333,7 @@
         end
     end
 
-    @testset "Conversion to transverse fields" begin
+    @testset "Various transverse fields" begin
         @field(A) do
             λ = 800u"nm"
             I₀ = 1e13u"W/cm^2"
@@ -370,12 +370,46 @@
         ApB = A+B
         CpD = C+D
 
-        tA = transverse_field(A)
-        @test tA isa ElectricFields.TransverseField
-        @test transverse_field(B) === B
-        @test transverse_field(ApB) === ApB
+        @testset "Conversion to transverse fields" begin
+            tA = transverse_field(A)
+            @test tA isa ElectricFields.TransverseField
+            @test transverse_field(B) === B
+            @test transverse_field(ApB) === ApB
 
-        tCpD = transverse_field(CpD)
-        @test tCpD isa ElectricFields.LinearTransverseField
+            tCpD = transverse_field(CpD)
+            @test tCpD isa ElectricFields.LinearTransverseField
+        end
+
+        @testset "Rotation of fields" begin
+            R = [1 0 0; 0 0 1; 0 1 0]
+
+            rA = rotate(A, R)
+            @test rA isa ElectricFields.TransverseField
+            @test rotation_matrix(rA) ≈ R
+            tA = timeaxis(A)
+            FA = field_amplitude(A, tA)
+            FrA = field_amplitude(rA, tA)
+            @test FrA[:,2] ≈ FA
+            @test FrA[:,3] ≈ zeros(length(tA)) atol=1e-14
+
+            rB = rotate(B, R)
+            @test rB isa ElectricFields.TransverseField
+            @test rotation_matrix(rB) ≈ R
+            tB = timeaxis(B)
+            FB = field_amplitude(B, tB)
+            FrB = field_amplitude(rB, tB)
+            @test FrB[:,1] ≈ FB[:,1]
+            @test FrB[:,2] ≈ FB[:,3]
+            @test FrB[:,3] ≈ zeros(length(tB)) atol=1e-14
+
+            rCpD = rotate(CpD, R)
+            @test rCpD isa ElectricFields.LinearTransverseField
+            @test rotation_matrix(rCpD) ≈ R
+            tCpD = timeaxis(CpD)
+            FCpD = field_amplitude(CpD, tCpD)
+            FrCpD = field_amplitude(rCpD, tCpD)
+            @test FrCpD[:,2] ≈ FCpD
+            @test FrCpD[:,3] ≈ zeros(length(tCpD)) atol=1e-14
+        end
     end
 end