From 76f31c72b85e6fc430afec512c94ab430de5ea58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefanos=20Carlstr=C3=B6m?= Date: Tue, 1 Aug 2023 08:47:07 +0200 Subject: [PATCH] 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 | 2 +- docs/src/index.md | 2 +- src/arithmetic.jl | 17 +++-- src/field_types.jl | 94 +++++++++++++++++++++++-- src/rotations.jl | 2 + test/arithmetic.jl | 165 ++++++++++++++++++++++++++++++++++++++++++-- test/field_types.jl | 96 +++++++++++++++++++++++--- 7 files changed, 350 insertions(+), 28 deletions(-) diff --git a/Project.toml b/Project.toml index 35d24c5..9dd637c 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.0" +version = "0.2.1" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/docs/src/index.md b/docs/src/index.md index 3d08731..1643f8d 100644 --- a/docs/src/index.md +++ b/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 diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 5011a05..2f47e27 100644 --- a/src/arithmetic.jl +++ b/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 diff --git a/src/field_types.jl b/src/field_types.jl index ed99447..1d8a491 100644 --- a/src/field_types.jl +++ b/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 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 50c9b42..5045be6 100644 --- a/test/arithmetic.jl +++ b/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 diff --git a/test/field_types.jl b/test/field_types.jl index e5a07e6..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) @@ -200,8 +200,8 @@ @test duration(F) == 4 @test span(F) == 0..4 - @test isone(period(F)) - @test isone(max_frequency(F)) + @test isone(austrip(period(F))) + @test isone(austrip(max_frequency(F))) @test photon_energy(F) == 2π @test dimensions(F) == 1 @@ -283,8 +283,8 @@ # @test duration(F) == 4 @test span(F) == 0..4 - @test isone(period(F)) - @test isone(max_frequency(F)) + @test isone(austrip(period(F))) + @test isone(austrip(max_frequency(F))) @test photon_energy(F) == 2π @test dimensions(F) == 1 @@ -332,4 +332,84 @@ @test Fv ≈ Fv2 rtol=1e-14 end end + + @testset "Various transverse 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 + + @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 + 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