More rotation fixes (#28)

* Allow SMatrix rotation matrices

* Rotate parent field(s) when possible; test NegatedField

* Bump version

* Fix ambiguity

Project.toml

@@ -1,7 +1,7 @@
 name = "ElectricFields"
 uuid = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca"
 authors = ["Stefanos Carlström <[email protected]>"]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/arithmetic.jl

@@ -146,6 +146,8 @@ dimensions(f::SumField) = dimensions(f.a)
 phase_shift(f::SumField, ϕ) =
     SumField(phase_shift(f.a, ϕ), phase_shift(f.b, ϕ))
 
+rotate(f::SumField, R) = SumField(rotate(f.a, R), rotate(f.b, R))
+
 # ** Wrapped fields
 
 """
@@ -159,7 +161,7 @@ for fun in [:params, :carrier, :envelope, :polarization,
             :wavelength, :period, :frequency, :max_frequency,
             :wavenumber, :fundamental, :photon_energy,
             :intensity, :amplitude, :duration, :continuity,
-            :span, :dimensions]
+            :span, :dimensions, :rotation_matrix]
     @eval $fun(f::WrappedField, args...) = $fun(parent(f), args...)
 end
 # [:vector_potential, :field_amplitude], should these be explicitly forwarded?
@@ -226,6 +228,8 @@ for fun in [:vector_potential, :vector_potential_spectrum]
     @eval $fun(f::NegatedField, t::Number) = -$fun(parent(f), t)
 end
 
+rotate(f::NegatedField, R) = NegatedField(rotate(f.a, R))
+
 # ** Delayed fields
 
 """
@@ -312,6 +316,8 @@ end
 
 Base.parent(f::DelayedField) = f.a
 
+rotate(f::DelayedField, R) = DelayedField(rotate(f.a, R), f.t₀)
+
 # *** DONE Delay operators
 #     Convention for delayed fields: a field delayed by a /positive/
 #     time, comes /later/, i.e. we write \(f(t-\delta t)\).
@@ -406,6 +412,8 @@ end
 time_integral(f::PaddedField) =
     time_integral(parent(f))
 
+rotate(f::PaddedField, R) = PaddedField(rotate(f.field, R), f.a, f.b)
+
 export PaddedField
 
 # ** Windowed fields
@@ -482,7 +490,7 @@ phase_shift(f::WindowedField, δϕ) =
     WindowedField(phase_shift(parent(f), δϕ), f.a, f.b)
 
 for fun in [:vector_potential, :field_amplitude, :intensity]
-    @eval function $fun(f::WindowedField{T}, t) where T
+    @eval function $fun(f::WindowedField{T}, t::Number) where T
         v = $fun(parent(f), t)
         t < f.a || t > f.b ? zero(v) : v
     end
@@ -497,9 +505,11 @@ function field_amplitude(f::WindowedField, a, b)
     end
 end
 
-function timeaxis(f::WindowedField, fs)
+function timeaxis(f::WindowedField, fs::Number)
     t = timeaxis(f.field, fs)
     t[findall(in(f.a..f.b), t)]
 end
 
+rotate(f::WindowedField, R) = WindowedField(rotate(f.field, R), f.a, f.b)
+
 export WindowedField

src/rotations.jl

@@ -19,7 +19,7 @@ function compute_rotation(R::AbstractMatrix{T}) where {T<:Real}
         throw(DimensionMismatch("Rotation matrix must have dimensions 3×3"))
     rank(R) < 3 &&
         throw(ArgumentError("Rotation matrix singular"))
-    R = copy(R)
+    R = Matrix{T}(R)
     for i = 2:3
         for j = 1:i-1
             R[:,i] -= dot(R[:,j],R[:,i])*R[:,j]

test/arithmetic.jl

@@ -1,10 +1,20 @@
+function get_me_three_dimensions(V)
+    if ndims(V) == 1
+        m = length(V)
+        hcat(zeros(m, 2), V)
+    elseif ndims(V) == 2
+        @assert size(V,2) == 3
+        V
+    else
+        error("This does not work")
+    end
+end
+
 function test_addition(::LinearPolarization, A, B)
     F = A + B
     F2 = B + A
     t = timeaxis(F)
 
-    tplot = 24.2e-3t
-
     Fv = field_amplitude(F, t)
     Fv2 = field_amplitude(F2, t)
     FvA = field_amplitude(A, t)
@@ -19,20 +29,6 @@ function test_addition(::ArbitraryPolarization, A, B)
     F2 = B + A
     t = timeaxis(F)
 
-    tplot = 24.2e-3t
-
-    function get_me_three_dimensions(V)
-        if ndims(V) == 1
-            m = length(V)
-            hcat(zeros(m, 2), V)
-        elseif ndims(V) == 2
-            @assert size(V,2) == 3
-            V
-        else
-            error("This does not work")
-        end
-    end
-
     Fv = get_me_three_dimensions(field_amplitude(F, t))
     Fv2 = get_me_three_dimensions(field_amplitude(F2, t))
     FvA = get_me_three_dimensions(field_amplitude(A, t))
@@ -44,6 +40,15 @@ end
 
 test_addition(A, B) = test_addition(polarization(A+B), A, B)
 
+function test_rotated_field(f, R)
+    R = ElectricFields.compute_rotation(R)
+    rf = rotate(f, R)
+    t = timeaxis(rf)
+    @test t ≈ timeaxis(f)
+    @test field_amplitude(rf, t) ≈ transpose(R*transpose(get_me_three_dimensions(field_amplitude(f, t))))
+    rf
+end
+
 @testset "Field arithmetic" begin
     @field(A) do
         I₀ = 1.0
@@ -133,6 +138,40 @@ test_addition(A, B) = test_addition(polarization(A+B), A, B)
                                │   – and a bandwidth of 0.3925 Ha = 10.6797 eV ⟺ 2.5823 PHz ⟺ 8.5598 Bohr = 452.9627 pm
                                └   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm"""
         end
+
+        @testset "Rotations" begin
+            rC = test_rotated_field(C, [1 0 0; 0 1 1; 0 -1 1])
+            @test rC isa ElectricFields.SumField
+        end
+    end
+
+    @testset "Negated fields" begin
+        nA = -A
+        @test nA isa ElectricFields.NegatedField
+        @test parent(nA) === A
+        t = timeaxis(nA)
+        @test t ≈ timeaxis(A)
+        @test field_amplitude(nA, t) ≈ -field_amplitude(A, t)
+
+        @field(B) do
+            I₀ = 0.5
+            T = 1.0
+            σ = 3.0
+            Tmax = 3.0
+        end
+
+        C = A - B
+        @test C.b isa ElectricFields.NegatedField
+        t = timeaxis(C)
+        @test field_amplitude(C, t) ≈ field_amplitude(A, t) - field_amplitude(B, t)
+
+        R = [1 0 0; 0 1 1; 0 -1 1]
+        rnA = test_rotated_field(nA, R)
+        @test rnA isa ElectricFields.NegatedField
+        @test rotation_matrix(rnA) ≈ ElectricFields.compute_rotation(R)
+        rC = test_rotated_field(C, R)
+        @test rC isa ElectricFields.SumField
+        @test rC.b isa ElectricFields.NegatedField
     end
 
     @testset "Delayed fields" begin
@@ -145,6 +184,11 @@ test_addition(A, B) = test_addition(polarization(A+B), A, B)
 
         @test span(B) == -5.6..6.4
 
+        R = [1 0 0; 0 1 1; 0 -1 1]
+        rB = test_rotated_field(B, R)
+        @test rB isa ElectricFields.DelayedField
+        @test rotation_matrix(rB) ≈ ElectricFields.compute_rotation(R)
+
         withenv("UNITFUL_FANCY_EXPONENTS" => true) do
             @test string(B) == """
                                Linearly polarized field with
@@ -215,6 +259,11 @@ test_addition(A, B) = test_addition(polarization(A+B), A, B)
 
         @test time_integral(B) == time_integral(A)
 
+        R = [1 0 0; 0 1 1; 0 -1 1]
+        rB = test_rotated_field(B, R)
+        @test rB isa ElectricFields.PaddedField
+        @test rotation_matrix(rB) ≈ ElectricFields.compute_rotation(R)
+
         withenv("UNITFUL_FANCY_EXPONENTS" => true) do
             @test string(B) == """
                                Padding before 124.0241 jiffies = 3.0000 fs and after 330.7310 jiffies = 8.0000 fs of
@@ -261,6 +310,11 @@ test_addition(A, B) = test_addition(polarization(A+B), A, B)
         @test field_amplitude(phase_shift(B, 2.0), 0.4) == field_amplitude(phase_shift(A, 2.0), 0.4)
         @test field_amplitude(phase_shift(B, 2.0), 3.0) == zero(field_amplitude(phase_shift(A, 2.0), 0.4))
 
+        R = [1 0 0; 0 1 1; 0 -1 1]
+        rB = test_rotated_field(B, R)
+        @test rB isa ElectricFields.WindowedField
+        @test rotation_matrix(rB) ≈ ElectricFields.compute_rotation(R)
+
         withenv("UNITFUL_FANCY_EXPONENTS" => true) do
             @test string(B) == """
                                Window from -4.1341 jiffies = -100.0000 as to 2.0671 jiffies = 50.0000 as of

test/field_types.jl

@@ -403,8 +403,9 @@
             @test FrB[:,3] ≈ zeros(length(tB)) atol=1e-14
 
             rCpD = rotate(CpD, R)
-            @test rCpD isa ElectricFields.LinearTransverseField
-            @test rotation_matrix(rCpD) ≈ R
+            @test rCpD isa ElectricFields.SumField
+            @test rotation_matrix(rCpD.a) ≈ R
+            @test rotation_matrix(rCpD.b) ≈ R
             tCpD = timeaxis(CpD)
             FCpD = field_amplitude(CpD, tCpD)
             FrCpD = field_amplitude(rCpD, tCpD)

test/rotations.jl

@@ -30,6 +30,9 @@ import ElectricFields: compute_rotation, rotation_angle, rotation_axis
                                 0.0 1.0 1.0
                                 0.0 0.0 1.0]) ≈ I rtol=1e-14
 
+        R = compute_rotation([1 0 0; 0 1 1; 0 -1 1])
+        @test compute_rotation(R) ≈ R rtol=1e-14
+
         @test_throws DimensionMismatch compute_rotation([1.0 0.0
                                                          0.0 1.0])