Conversion to transverse fields, rotations (#27)

* Implemented addition of linearly and arbitrarily polarized fields

* Bumped version

* Convert any field to transverse

* Avoid extra newline when printin SumFields

* Fixed units of period & max_frequency for ConstantField & Ramp

* Implemented rotation of fields, post facto

* Fix doctests

Project.toml

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

docs/src/index.md

@@ -95,7 +95,7 @@ julia> instantaneous_intensity(IR, 4.0), intensity(IR, 4.0)
 (0.0026984944767521166, 0.002847529708372214)
 
 julia> span(IR)
--661.9198939608042..661.9198939608042
+-661.9198939608042 .. 661.9198939608042
 
 julia> timeaxis(IR)
 -661.9198939608042:1.1041199232040104:661.9198939608042

src/arithmetic.jl

@@ -98,18 +98,17 @@ function show(io::IO, f::SumField)
         write(io, "$s $l\n")
     end
 
-    write(io, "⊕\n")
+    write(io, "⊕")
 
     for (s,l) in zip(repeat("│", length(b_str)-1) * "└", b_str)
-        write(io, "$s $l\n")
+        write(io, "\n$s $l")
     end
 end
 
-function +(a::AbstractField, b::AbstractField)
-    polarization(a) == polarization(b) ||
-        throw(ArgumentError("Cannot add fields of different polarization"))
-    SumField(a, b)
-end
++(::Pol, a::AbstractField, ::Pol, b::AbstractField) where {Pol<:Polarization} = SumField(a, b)
++(::LinearPolarization, a::AbstractField, ::ArbitraryPolarization, b::AbstractField) = SumField(transverse_field(a), b)
++(::ArbitraryPolarization, a::AbstractField, ::LinearPolarization, b::AbstractField) = SumField(a, transverse_field(b))
++(a::AbstractField, b::AbstractField) = +(polarization(a), a, polarization(b), b)
 
 for fun in [:vector_potential, :field_amplitude, :vector_potential_spectrum]
     @eval $fun(f::SumField, t::Number) =
@@ -364,7 +363,7 @@ Linearly polarized field with
   – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
 
 julia> span(A), span(B)
-(-6.0..6.0, -16.0..36.0)
+(-6.0 .. 6.0, -16.0 .. 36.0)
 ```
 """
 struct PaddedField{Field<:AbstractField,T} <: WrappedField
@@ -447,7 +446,7 @@ Linearly polarized field with
   – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm
 
 julia> span(A), span(B)
-(-6.0..6.0, -3.0..5.0)
+(-6.0 .. 6.0, -3.0 .. 5.0)
 
 julia> field_amplitude(A, -4)
 -0.6395632362683398

src/field_types.jl

@@ -683,6 +683,79 @@ end
 
 time_integral(f::TransverseField) = time_integral(envelope(f))
 
+Base.convert(::Type{<:TransverseField}, f::LinearField) =
+    TransverseField(LinearTransverseCarrier(carrier(f)),
+                    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"""
+LinearTransverseField
+
+Linearly polarized field, but explicitly in 3d, i.e the polarization
+vector is confined to a line in the plane that is perpendicular to the
+direction of propagation. If the propagation vector is not rotated,
+this plane is the ``z-x``, with ``z`` by convention being the
+principal polarization axis, and ``y`` is the direction of
+propagation.
+
+This is mostly used as a proxy, if one wants to calculate using a
+manifestly linearly polarized field in 3d. See also
+[`LinearTransverseCarrier`](@ref).
+"""
+struct LinearTransverseField{LinearField<:AbstractField,Rotation} <: AbstractField
+    linear_field::LinearField
+    R::Rotation
+end
+
+LinearTransverseField(::LinearPolarization, f::AbstractField, R=I) =
+    LinearTransverseField(f, compute_rotation(R))
+
+LinearTransverseField(f::AbstractField, args...) =
+    LinearTransverseField(polarization(f), f, args...)
+
+vector_potential(f::LinearTransverseField, t::Number) =
+    f.R*SVector(0, 0, vector_potential(f.linear_field, t))
+
+vector_potential_spectrum(f::LinearTransverseField, ω::Number) =
+    f.R*SVector(0, 0, vector_potential_spectrum(f.linear_field, ω))
+
+for fun in (:vector_potential, :intensity, :amplitude, :carrier, :envelope, :params, :max_frequency, :span)
+    @eval $fun(f::LinearTransverseField) = $fun(f.linear_field)
+end
+
+rotation_matrix(f::LinearTransverseField) = f.R
+
+function show(io::IO, f::LinearTransverseField)
+    print(io, """
+Linearly polarized field in transverse plane constructed from
+""")
+    show(io, f.linear_field)
+    isrotated = f.R isa AbstractMatrix
+    if isrotated
+        printfmt(io, "\n  – and a rotation of {1:.2f}π about [{2:.3f}, {3:.3f}, {4:.3f}]",
+                 rotation_angle(f.R)/π, rotation_axis(f.R)...)
+    end
+end
+
+polarization(::LinearTransverseField) = ArbitraryPolarization()
+
+dimensions(::LinearTransverseField) = 3
+
+phase_shift(f::LinearTransverseField, δϕ) =
+    LinearTransverseField(phase_shift(f.linear_field, δϕ), f.R)
+
+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"""
@@ -777,8 +850,8 @@ span(f::ConstantField{T}) where T = zero(T)..f.tmax
 
 # These methods are just for convenience (to be able to establish a
 # time base for calculations), but they are not physical as such.
-period(f::ConstantField{T}) where T = one(T)
-max_frequency(f::ConstantField{T}) where T = one(T)
+period(f::ConstantField{T}) where T = one(T)*aunit(u"s")
+max_frequency(f::ConstantField{T}) where T = one(T)/aunit(u"s")
 photon_energy(f::ConstantField{T}) where T = 2T(π)
 
 dimensions(::ConstantField) = 1
@@ -928,12 +1001,21 @@ span(f::Ramp{T}) where T = zero(T)..f.tmax
 
 # These methods are just for convenience (to be able to establish a
 # time base for calculations), but they are not physical as such.
-period(f::Ramp{T}) where T = one(T)
-max_frequency(f::Ramp{T}) where T = one(T)
+period(f::Ramp{T}) where T = one(T)*aunit(u"s")
+max_frequency(f::Ramp{T}) where T = one(T)/aunit(u"s")
 photon_energy(f::Ramp{T}) where T = 2T(π)
 
 dimensions(::Ramp) = 1
 
+# * Conversion to transverse field
+
+transverse_field(::LinearPolarization, f::LinearField) = convert(TransverseField, f)
+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,
@@ -947,4 +1029,6 @@ export LinearPolarization, ArbitraryPolarization, polarization,
     field_envelope,
     params, dimensions,
     phase_shift, phase,
-    field_amplitude_spectrum, vector_potential_spectrum
+    field_amplitude_spectrum, vector_potential_spectrum,
+    transverse_field,
+    rotate, rotation_matrix

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

test/arithmetic.jl

@@ -1,3 +1,49 @@
+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)
+    FvB = field_amplitude(B, t)
+
+    @test Fv  ≈ FvA + FvB
+    @test Fv2 ≈ FvA + FvB
+end
+
+function test_addition(::ArbitraryPolarization, A, B)
+    F = 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))
+    FvB = get_me_three_dimensions(field_amplitude(B, t))
+
+    @test Fv  ≈ FvA + FvB
+    @test Fv2 ≈ FvA + FvB
+end
+
+test_addition(A, B) = test_addition(polarization(A+B), A, B)
+
 @testset "Field arithmetic" begin
     @field(A) do
         I₀ = 1.0
@@ -65,8 +111,7 @@
                                │   – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV, f = 41.3414 PHz)
                                │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
                                │   – 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
-                               """
+                               └   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm"""
         end
         withenv("UNITFUL_FANCY_EXPONENTS" => false) do
             @test string(C) == """
@@ -86,8 +131,7 @@
                                │   – a Fixed carrier @ λ = 7.2516 nm (T = 24.1888 as, ω = 6.2832 Ha = 170.9742 eV, f = 41.3414 PHz)
                                │   – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±1.00σ)
                                │   – 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
-                               """
+                               └   – Uₚ = 0.0032 Ha = 86.1591 meV => α = 0.0179 Bohr = 947.8211 fm"""
         end
     end
 
@@ -242,4 +286,117 @@
                                  – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
         end
     end
+
+    @testset "Add various kinds of fields" begin
+        @field(A) do
+            λ = 800u"nm"
+            I₀ = 1e13u"W/cm^2"
+            τ = 1.45u"fs"
+            σoff = 4.0
+            σmax = 6.0
+            env = :trunc_gauss
+        end
+
+        @field(B) do
+            λ = 100u"nm"
+            I₀ = 1e12u"W/cm^2"
+            τ = 1.45u"fs"
+            σoff = 4.0
+            σmax = 6.0
+            env = :trunc_gauss
+            ξ = 1.0
+        end
+
+        @field(C) do
+            tmax = 3.0u"fs"
+            E₀ = 0.1
+            kind = :constant
+        end
+        C = delay(C, -3.0u"fs")
+
+        @field(D) do
+            tmax = 3.0u"fs"
+            E₀ = 0.1
+            kind = :sin²_ramp
+            ramp = :down
+        end
+
+        ApB = A+B
+        CpD = C+D
+
+        F = A + B
+        F2 = B + A
+
+        @test polarization(F) == ArbitraryPolarization()
+        @test polarization(F2) == ArbitraryPolarization()
+
+        test_addition(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
+
+        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
+            @test string(ApBpCpD) == """
+                  ┌ ┌ Transversely polarized field with
+                  │ │   - I₀ = 2.8495e-04 au = 1.0e13 W cm⁻² =>
+                  │ │     - E₀ = 1.6880e-02 au = 8.6802 GV m⁻¹
+                  │ │     - A₀ = 0.2964 au
+                  │ │   – a LinearTransverseCarrier: Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz)
+                  │ │   – and a Truncated Gaussian envelope of duration 59.9450 jiffies = 1.4500 fs (intensity FWHM; turn-off from 2.4630 fs to 3.6945 fs)
+                  │ │   – and a bandwidth of 0.0463 Ha = 1.2586 eV ⟺ 304.3250 THz ⟺ 12277.0979 Bohr = 649.6760 nm
+                  │ │   – Uₚ = 0.0220 Ha = 597.5865 meV => α = 5.2039 Bohr = 275.3788 pm
+                  │ ⊕
+                  │ │ Transversely polarized field with
+                  │ │   - I₀ = 2.8495e-05 au = 1.0e12 W cm⁻² =>
+                  │ │     - E₀ = 5.3380e-03 au = 2.7449 GV m⁻¹
+                  │ │     - A₀ = 0.0117 au
+                  │ │   – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 100.0000 nm (T = 333.5641 as, ω = 0.4556 Ha = 12.3984 eV, f = 2.9979 PHz)
+                  │ │   – and a Truncated Gaussian envelope of duration 59.9450 jiffies = 1.4500 fs (intensity FWHM; turn-off from 2.4630 fs to 3.6945 fs)
+                  │ │   – and a bandwidth of 0.0463 Ha = 1.2586 eV ⟺ 304.3250 THz ⟺ 191.8297 Bohr = 10.1512 nm
+                  │ └   – Uₚ = 0.0000 Ha = 933.7290 μeV => α = 0.0257 Bohr = 1.3607 pm
+                  ⊕
+                  │ Linearly polarized field in transverse plane constructed from
+                  │ ┌ Constant field of
+                  │ │   - 124.0241 jiffies = 3.0000 fs duration, and
+                  │ │   - E₀ = 1.0000e-01 au = 51.4221 GV m⁻¹
+                  │ │   – delayed by -124.0241 jiffies = -3.0000 fs
+                  │ ⊕
+                  │ │ sin² down-ramp of
+                  │ │   - 124.0241 jiffies = 3.0000 fs duration, and
+                  └ └   - E₀ = 1.0000e-01 au = 51.4221 GV m⁻¹"""
+        end
+        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
+            @test string(ApBpCpD) == """
+                  ┌ ┌ Transversely polarized field with
+                  │ │   - I₀ = 2.8495e-04 au = 1.0e13 W cm^-2 =>
+                  │ │     - E₀ = 1.6880e-02 au = 8.6802 GV m^-1
+                  │ │     - A₀ = 0.2964 au
+                  │ │   – a LinearTransverseCarrier: Fixed carrier @ λ = 800.0000 nm (T = 2.6685 fs, ω = 0.0570 Ha = 1.5498 eV, f = 374.7406 THz)
+                  │ │   – and a Truncated Gaussian envelope of duration 59.9450 jiffies = 1.4500 fs (intensity FWHM; turn-off from 2.4630 fs to 3.6945 fs)
+                  │ │   – and a bandwidth of 0.0463 Ha = 1.2586 eV ⟺ 304.3250 THz ⟺ 12277.0979 Bohr = 649.6760 nm
+                  │ │   – Uₚ = 0.0220 Ha = 597.5865 meV => α = 5.2039 Bohr = 275.3788 pm
+                  │ ⊕
+                  │ │ Transversely polarized field with
+                  │ │   - I₀ = 2.8495e-05 au = 1.0e12 W cm^-2 =>
+                  │ │     - E₀ = 5.3380e-03 au = 2.7449 GV m^-1
+                  │ │     - A₀ = 0.0117 au
+                  │ │   – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 100.0000 nm (T = 333.5641 as, ω = 0.4556 Ha = 12.3984 eV, f = 2.9979 PHz)
+                  │ │   – and a Truncated Gaussian envelope of duration 59.9450 jiffies = 1.4500 fs (intensity FWHM; turn-off from 2.4630 fs to 3.6945 fs)
+                  │ │   – and a bandwidth of 0.0463 Ha = 1.2586 eV ⟺ 304.3250 THz ⟺ 191.8297 Bohr = 10.1512 nm
+                  │ └   – Uₚ = 0.0000 Ha = 933.7290 μeV => α = 0.0257 Bohr = 1.3607 pm
+                  ⊕
+                  │ Linearly polarized field in transverse plane constructed from
+                  │ ┌ Constant field of
+                  │ │   - 124.0241 jiffies = 3.0000 fs duration, and
+                  │ │   - E₀ = 1.0000e-01 au = 51.4221 GV m^-1
+                  │ │   – delayed by -124.0241 jiffies = -3.0000 fs
+                  │ ⊕
+                  │ │ sin² down-ramp of
+                  │ │   - 124.0241 jiffies = 3.0000 fs duration, and
+                  └ └   - E₀ = 1.0000e-01 au = 51.4221 GV m^-1"""
+        end
+    end
 end