From cfb79574bb91229aeea601e827ab8eafe8b0b25a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanos=20Carlstr=C3=B6m?= Date: Tue, 1 Aug 2023 20:59:09 +0200 Subject: [PATCH] More rotation fixes (#28) * Allow SMatrix rotation matrices * Rotate parent field(s) when possible; test NegatedField * Bump version * Fix ambiguity --- Project.toml | 2 +- src/arithmetic.jl | 16 +++++++-- src/rotations.jl | 2 +- test/arithmetic.jl | 86 ++++++++++++++++++++++++++++++++++++--------- test/field_types.jl | 5 +-- test/rotations.jl | 3 ++ 6 files changed, 91 insertions(+), 23 deletions(-) diff --git a/Project.toml b/Project.toml index 9dd637c..3e50f37 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ElectricFields" uuid = "2f84ce32-9bb1-11e8-0d9f-3dce90a4beca" authors = ["Stefanos Carlström "] -version = "0.2.1" +version = "0.2.2" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 2f47e27..18d0515 100644 --- a/src/arithmetic.jl +++ b/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 diff --git a/src/rotations.jl b/src/rotations.jl index ba1892f..5bb383b 100644 --- a/src/rotations.jl +++ b/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] diff --git a/test/arithmetic.jl b/test/arithmetic.jl index 5045be6..78bbf48 100644 --- a/test/arithmetic.jl +++ b/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 diff --git a/test/field_types.jl b/test/field_types.jl index 0bfd8c8..518d6b1 100644 --- a/test/field_types.jl +++ b/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) diff --git a/test/rotations.jl b/test/rotations.jl index 09f5c16..9840aa0 100644 --- a/test/rotations.jl +++ b/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])